Kopia_notatnika_nw_skarbek_affine_fit_pl.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/TayTT/AI_exercises/blob/main/Kopia_notatnika_nw_skarbek_affine_fit_pl.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XlGnsLIBtaG0"
      },
      "source": [
        "copyright 2020: *Władysław Skarbek, Warsaw University of Technology*<br>\n",
        "(Jupyter notebook addressed to students of: EADIR, GSN, MATMU)\n",
        "$\\newcommand{\\xx}[2]{(#1)(#2)}$\n",
        "$\\newcommand{\\bb}[1]{\\mathbb{#1}}$\n",
        "$\\newcommand{\\cl}[1]{\\mathcal{#1}}$\n",
        "$\\newcommand{\\tp}[1]{{#1}^{\\intercal}}$\n",
        "$\\newcommand{\\tr}[1]{\\text{trace}\\left[#1\\right]}$\n",
        "$\\newcommand{\\inv}[1]{\\in\\bb{R}^{#1}}$\n",
        "$\\newcommand{\\inm}[2]{\\in\\bb{R}^{#1\\times#2}}$\n",
        "$\\newcommand{\\invc}[1]{\\in\\bb{C}^{#1}}$\n",
        "$\\newcommand{\\inmc}[2]{\\in\\bb{C}^{#1\\times#2}}$\n",
        "$\\newcommand{\\rbox}[2]{\\mathrel{\\raise{#1}{#2}}}$\n",
        "$\\newcommand{\\xconv}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{2.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-4pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\bb{C}}\\kern{-6.5pt}\\raise{4.5pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{0pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{2.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\newcommand{\\xdense}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{1.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-4pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\bb{F}}\\kern{-4.5pt}\\raise{1.5pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{2pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{1.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\newcommand{\\xpool}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{1.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-3pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\bb{P}}\\kern{-6.5pt}\\raise{0.8pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{0pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{1.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\newcommand{\\xinp}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{1.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-4pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\cl{I}}\\kern{-4.5pt}\\raise{1.5pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{2pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{1.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\newcommand{\\xdrop}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{1.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-4pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\bb{D}}\\kern{-4.5pt}\\raise{1.5pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{0.5pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{1.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\newcommand{\\xmerge}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{1.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-4pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\bb{M}}\\kern{-4.5pt}\\raise{1.5pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{2pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{1.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\newcommand{\\xgeneral}[5]{\n",
        "\\begin{array}[]{r}{\\tiny{#4\\ }}\\\\\\raise{1.5pt}{{\\tiny{#3}}}\\end{array}\n",
        "\\kern{-4pt}\n",
        "\\begin{array}[]{l}{\\LARGE\\bb{Q}}\\kern{-9.0pt}\\raise{4.5pt}{{\\tiny{#5}}}\\end{array}\n",
        "\\kern{2pt}\n",
        "\\begin{array}[]{l}{\\tiny {#2}}\\\\\\raise{1.5pt}{{\\tiny{#1}}}\\end{array}}$\n",
        "$\\def\\ds{\\displaystyle}$\n",
        "$\\def\\ass{\\leftarrow}$\n",
        "$\\def\\od#1#2{\\nabla_{#2}#1}$\n",
        "$\\def\\tod#1#2{\\tp{\\nabla}_{#2}{#1}}$\n",
        "$\\def\\cl#1{{\\cal#1}}$\n",
        "$\\def\\mbm#1{\\mathbf{#1}}$\n",
        "$\\def\\sp#1#2{\\frac{\\partial{#1}}{\\partial{#2}}}$\n",
        "$\\def\\spp#1#2#3{\\frac{\\partial{#1}}{\\partial{#2}\\partial{#3}}}$\n",
        "$\\def\\eqd{\\doteq}$\n",
        "$\\def\\ra{\\rightarrow}$\n",
        "$\\def\\lra{\\longrightarrow}$\n",
        "$\\def\\ovra#1{\\overset{#1}{\\lra}}$\n",
        "$\\def\\dra{\\overset{\\circ}{\\lra}}$\n",
        "$\\def\\xeq#1{\\overset{#1}{=}}$\n",
        "$\\def\\ov#1{\\overline{#1}}$\n",
        "$\\def\\dotp#1#2{\\left\\langle#1,#2\\right\\rangle}$\n",
        "$\\def\\th#1{\\ov{#1}^{\\intercal}}$\n",
        "$\\def\\rv#1{\\widetilde{#1}}$\n",
        "$\\def\\vars#1#2{\\mathtt{var}_{#1}{\\left[#2\\right]}}$\n",
        "$\\def\\pmodd#1{\\kern{-1mm}\\pmod{#1}}$\n",
        "$\\def\\pdt#1#2#3{\\frac{\\partial^2 #1}{\\partial #2\\partial #3}}$\n",
        "$\\def\\diag#1{\\mathtt{diag}\\left[#1\\right]}$\n",
        "$\\def\\rank#1{\\mathtt{rank}\\left[#1\\right]}$\n",
        "$\\def\\span#1{\\mathtt{span}\\left[#1\\right]}$\n",
        "$\\def\\spann#1{\\mathtt{span}\\left[#1\\right]}$\n",
        "$\\def\\tr#1{\\mathtt{tr}\\left[#1\\right]}$\n",
        "$\\def\\var#1{\\mathtt{var}\\left[#1\\right]}$\n",
        "$\\def\\kernn#1{\\kern{#1pt}}$\n",
        "$\\def\\raisee#1#2{\\raise{#1}{#2}}$\n",
        "$\\def\\colorr#1#2{\\color{#1}{#2}}$\n",
        "$\\def\\wtt{\\ \\text{ wtt }\\ }$\n",
        "$\\def\\mbm#1{\\mathbf{#1}}$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3DTnQZdyG0so"
      },
      "source": [
        "# Matematyka afinicznego dopasowania\n",
        "\n",
        "W szczególności narzędzia stosowane w dopasowaniu chmur punktów w $R^d$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GJ_Z7kADqJTS"
      },
      "source": [
        "###Dekompozycja EVD\n",
        "\n",
        "Powtórka z elementów klasycznej algebry liniowej $-$ wiek XIX."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YLanzZMMThpi"
      },
      "source": [
        "####Wartość własna $-$ definicja spektrum macierzy\n",
        "\n",
        "*Spektrum* macierzy zespolonej $C\\inmc{n}{n}$ stanowi zbiór liczb zespolonyc $\\lambda(C)\\subset\\bb{C}$ zdefiniowany następująco:\n",
        "$\n",
        "\\lambda(C) \\eqd \\{\\lambda\\in\\bb{C}: \\det(C-\\lambda I_n)=0\\}\n",
        "$\n",
        "- Funkcję  $p(\\lambda)\\eqd\\det(C-\\lambda I_n)$ nazywamy *wielomianem charakterystycznym macierzy* $C.$\n",
        "- *Spektrum jest zawsze zbiorem skończonym:*:<br>Jeśli $C\\inmc{n}{n}$, to $\\left|\\lambda(C)\\right|\\leq n.$\n",
        "- Elementy spektrum macierzy nazywamy jej *wartościami własnymi*: $\\lambda\\in\\lambda(C)$ wtedy i tylko wtedy gdy $\\exists u\\invc{n}, Cu=\\lambda u.$\n",
        "- Wektor $u$ skalowany w transformację $C$ przez czynnik $\\lambda$ nazywamu *wektorem własnym* macierzy $C$ określonym dla wartości własnej $\\lambda.$\n",
        "- Każda wartość własna określa *podprzestrzeń własną*:<br>\n",
        "Jesli $\\lambda\\in\\lambda(C)$, to zbiór $\\cl{E}_{\\lambda}\\eqd\\{u: Cu=\\lambda u\\}$ jest podprzestrzenią liniową w przestrzeni $\\bb{C}^n.$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OXe36kckgaIi"
      },
      "source": [
        "####Obliczanie wartości własnych w pakiecie `sympy`"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "po5AGYUnhYDK"
      },
      "source": [
        "#####Wielomian w reprezentacji symbolicznej: rozwijanie, wyznaczniki, pierwiastki"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "js3GA8kkg6kL",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6e838711-9243-48bc-95af-c2201691a737"
      },
      "source": [
        "import sympy as sym\n",
        "x = sym.Symbol('x'); y = sym.Symbol('y')\n",
        "print(sym.expand((x + y) ** 3))\n",
        "print(sym.expand((x-1)*(x+1)))"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "x**3 + 3*x**2*y + 3*x*y**2 + y**3\n",
            "x**2 - 1\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nP9unskkjHFf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c11ff633-ea28-4f25-bb7d-0ae75561a415"
      },
      "source": [
        "As = sym.Matrix([[1,2],[3,4]])\n",
        "xI = sym.Matrix([[x,0],[0,x]])\n",
        "ex = sym.expand(xI-As); print(ex)\n",
        "d = ex.det(); print(d)\n",
        "p = sym.expand(d); print(p)\n",
        "roots = sym.nroots(p)\n",
        "roots"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Matrix([[x - 1, -2], [-3, x - 4]])\n",
            "x**2 - 5*x - 2\n",
            "x**2 - 5*x - 2\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[-0.372281323269014, 5.37228132326901]"
            ]
          },
          "metadata": {},
          "execution_count": 2
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_JWVQCmiibAe"
      },
      "source": [
        "#####Wartość własna $-$ algorytm symboliczny versus algorytm numeryczny"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "D1H8ydlkuAug"
      },
      "source": [
        "import numpy as np\n",
        "from numpy.linalg import eig, svd, norm, cholesky"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AnAKbiWoi-GG",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e6ec6d21-5080-481d-f686-1befddda8f7c"
      },
      "source": [
        "A = np.array([[1,2],[3,4]])\n",
        "ls,vs = eig(A)\n",
        "A, ls, vs, np.dot(A,vs[:,0]), ls[0]*vs[:,0], ls, roots, ls-roots"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([[1, 2],\n",
              "        [3, 4]]),\n",
              " array([-0.37228132,  5.37228132]),\n",
              " array([[-0.82456484, -0.41597356],\n",
              "        [ 0.56576746, -0.90937671]]),\n",
              " array([ 0.30697009, -0.21062466]),\n",
              " array([ 0.30697009, -0.21062466]),\n",
              " array([-0.37228132,  5.37228132]),\n",
              " [-0.372281323269014, 5.37228132326901],\n",
              " array([0, 0], dtype=object))"
            ]
          },
          "metadata": {},
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hJ1ozMVkl18g"
      },
      "source": [
        "Wektory własne $-$ testujemy spełnienie własności definicyjnej"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3mrEjmC1i16Y",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "390345dc-4b22-4234-9b53-943faabc28c4"
      },
      "source": [
        "A = np.array([[1,-2,3],[-4,5,6],[7,8,-9]])\n",
        "ls, vs = eig(A)\n",
        "A, ls, vs, np.dot(A,vs[:,2]), ls[2]*vs[:,2]"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(array([[ 1, -2,  3],\n",
              "        [-4,  5,  6],\n",
              "        [ 7,  8, -9]]),\n",
              " array([-13.76691494,   2.78453226,   7.98238268]),\n",
              " array([[ 0.23127389,  0.84043331, -0.09343993],\n",
              "        [ 0.34064779,  0.05826955,  0.91498738],\n",
              "        [-0.91130208,  0.53877315,  0.3925138 ]]),\n",
              " array([-0.7458733 ,  7.30377939,  3.13319532]),\n",
              " array([-0.7458733 ,  7.30377939,  3.13319532]))"
            ]
          },
          "metadata": {},
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Zd05GZvIVKBZ"
      },
      "source": [
        "#### Problem własny w przypadku macierzy symetrycznych\n",
        "\n",
        "- *Transpozycja* $\\tp{C}$ macierzy zespolonej $C\\inmc{m}{n}$  w definicji po-elementowej:\n",
        "$\n",
        "(\\tp{C})_{ij}\\eqd C_{ji},\\ i=1,\\dots,m,\\ j=1,\\dots,n\n",
        "$\n",
        "- Definicja: Macierz $C\\inmc{m}{n}$ jest *symetryczna*  wtedy i tylko wtedy gdy $\\tp{C}=C.$\n",
        "- Własność: Jeśli $C\\inmc{m}{n}$ jest symetryczna, to $m=n$ oraz  $C\\inm{n}{n}$.\n",
        "- *Sprzężona transpozycja* $\\th{C}$ macierzy zespolonej $C$  w definicji po-elementowej:\n",
        "$\n",
        "(\\th{C})_{ij}\\eqd \\ov{C}_{ji},\\ i,j=1,\\dots,n\n",
        "$\n",
        "- $C$ jest *Hermitowska* wtedy i tylko wtedy gdy $\\th{C}=C.$\n",
        "- Własność: Spektra macierzy symetrycznej i macierzy Hermitowskiej są podzbiorami zbioru liczb rzeczywistych:<br>\n",
        "Jesli $C=\\tp{C}$ lub $C=\\th{C}$, to $\\lambda(C)\\subset\\bb{R}$\n",
        "- Własność: Dla różnych wartości własnych mamy *ortogonalność wektorów własnych*:<br>\n",
        "Jeśli $Cu_i=\\lambda_i u_i,\\ Cu_j=\\lambda_j u_j$ oraz $\\lambda_1\\neq\\lambda_2$, to $\\tp{u}_iu_j=0$ w przypadku symetrycznej macierzy $C$ oraz $\\th{u}_iu_j=0$ w przypadku macierzy Hermitowskiej.\n",
        "- Własność: Z każda macierzą symetryczną lub macierzą Hermitowską $C$ możemy związać *bazę ortogonalną* złożoną wyłącznie z wektorów własnych:<br>\n",
        "Niech $\\mathtt{obase}(\\cl{S})$ oznacza ortogonalną bazę w liniowej podprzestrzenie $\\cl{S}.$ Wtedy baza dla całej przestrzeni może być sumą mnogościową baz dla podprzestrzeni własnych wszystkich wartości własnych:\n",
        "$$\\ds\n",
        "\\mathtt{obase}\\left(\\bb{C}^n\\right) \\eqd \\bigcup_{\\lambda\\in\\lambda(C)}\\mathtt{obase}\\left(\\cl{E}_{\\lambda}\\right)\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jwC95vkFW6E7"
      },
      "source": [
        "####Dekompozycja EVD\n",
        "\n",
        "Niech ortogonalna baza wektorów własnych symetrycznej lub Hermitowskiej macierz $C\\inmc{n}{n}$ tworzy kolumny macierzy $U\\inmc{n}{n}$ zapisane według malejącego ciągu wartości własnych:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "Cu_i = \\lambda_iu_i,\\ i=1,\\dots,n\\\\[5pt]\n",
        "U \\eqd [u_1,\\dots,u_n],\\ \\lambda_1\\geq\\dots\\geq\\lambda_n,\\\n",
        "\\Lambda \\eqd \\diag{\\lambda_1,\\dots,\\lambda_n}\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "*Dekompozycja EVD* (EigenValue Decomposition) macierzy $C\\overset{evd}{=}U\\Lambda\\th{U}$ wynika z własności ortogonalności bazy własnej, a więc i macierzy własnej $U$:<br>\n",
        "Dekompozycja EVD jest treścią tzw. twierdzenia spektralnego:\n",
        "Dla macierzy symetrycznej symmetric $C\\inm{n}{n}$ oraz dla macierzy hermitowskiej $C\\inmc{n}{n}$\n",
        "$$\n",
        "C = U\\Lambda\\th{U} = \\sum_{i=1}^n\\lambda_iu_i\\th{u}_i\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ioWXQ096ugNW",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "cc821c0f-d08e-4171-f56b-c4b091adde60"
      },
      "source": [
        "B = np.dot(A,A.T)\n",
        "eigvals, eigvecs = eig(B)\n",
        "print(B, '\\n', np.dot(np.dot(eigvecs,np.diag(eigvals)),eigvecs.T))"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 14   4 -36]\n",
            " [  4  77 -42]\n",
            " [-36 -42 194]] \n",
            " [[ 14.   4. -36.]\n",
            " [  4.  77. -42.]\n",
            " [-36. -42. 194.]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m7XZ3GARCgIf"
      },
      "source": [
        "###Macierze nieujemnie określone\n",
        "\n",
        "*Definicje*:\n",
        "- Jeśli  $A\\inm{n}{n}$ oraz  $\\tp{A}=A$ to $A$ jest nieujemnie określona wtedy i tylko wtedy gdy $\\tp{x}Ax\\geq 0$ dla dowolnego $x\\inv{n}.$\n",
        "- Jeśli  $A\\inmc{n}{n}$ oraz  $\\th{A}=A$ to $A$ jest nieujemnie określona wtedy i tylko wtedy gdy $\\th{x}Ax\\geq 0$ dla dowolnego $x\\invc{n}.$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lFfBLWVopRYA"
      },
      "source": [
        "*Własności*\n",
        "\n",
        "- W symetrycznej i nieujemnie określonej macierzy $A\\inm{n}{n}$, wektory własne o różnych wartościach własnych są wzajemnie prostopadłe, tj. gdy $\\tp{A}=A, \\forall x\\inv{n}\\tp{x}Ax\\geq 0$:\n",
        "$$\n",
        "\\left[Au=\\lambda u,\\ Au'=\\lambda' u',\\ \\lambda'\\neq\\lambda\\right] \\lra \\tp{u}u' = 0\n",
        "$$\n",
        "- W przypadku nieujemnie okreslonej macierzy Hermitowskiej, tj. gdy $\\th{A}=A, \\forall x\\inv{n}\\th{x}Ax\\geq 0$:\n",
        "$$\n",
        "\\left[Au=\\lambda u, Au'=\\lambda' u',\\ \\lambda'\\neq\\lambda\\right] \\lra \\th{u}u' = 0\n",
        "$$\n",
        "- Jeśli $B\\inm{m}{n}$, to macierze $\\tp{B}B$ oraz $B\\tp{B}$ są symetryczne i nieujemnie określone.\n",
        "- Jeśli $C\\inmc{m}{n}$, to macierze $\\th{C}C$ oraz $C\\th{C}$ są Hermitowskie i nieujemnie określone.\n",
        "- Jeśli $A\\inm{m}{n}$, to spektra macierzy $A\\tp{A}$ oraz $\\tp{A}A$ są identyczne i co więcej wektory własne odpowiadają sobie wzajemnie w odwzorowaniach $\\tp{A}$ oraz $A$, odpowiednio.\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "A\\tp{A}u = \\lambda u \\lra \\tp{A}A\\underbrace{(\\tp{A}u)}_v = \\lambda\\underbrace{(\\tp{A}u)}_v\\\\\n",
        "\\tp{A}Av = \\lambda v \\lra A\\tp{A}\\underbrace{(Av)}_u = \\lambda\\underbrace{(Av)}_u\n",
        "\\end{array}\n",
        "$$\n",
        "- Jeśli $A\\inmc{m}{n}$, to spektra macierzy $A\\th{A}$ oraz $\\th{A}A$ są identyczne, co więcej wektory własne odpowiadają sobie wzajemnie w odwzorowaniach $\\th{A}$ oraz $A$, odpowiednio.\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "A\\th{A}u = \\lambda u \\lra \\th{A}A\\underbrace{(\\th{A}u)}_v = \\lambda\\underbrace{(\\th{A}u)}_v\\\\[5pt]\n",
        "\\th{A}Av = \\lambda v \\lra A\\th{A}\\underbrace{(Av)}_u = \\lambda\\underbrace{(Av)}_u\n",
        "\\end{array}\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dGEh7iy8tA8l"
      },
      "source": [
        "- Jeśli wektory własne $u,u'$ macierzy $A\\tp{A}$ są prostopadłe, to ich odpowiedniki $v\\eqd\\tp{A}u$ oraz $v'\\eqd\\tp{A}u'$  są też prostopadłe:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "\\left[A\\tp{A}u = \\lambda u,\\ v\\eqd\\tp{A}u,\\ A\\tp{A}u' = \\lambda' u',\\ v'\\eqd\\tp{A}u'\\ \\tp{u}u'=0\\right]\\\\\n",
        "\\lra \\tp{v}v'=\\tp{u}(A\\tp{A}u') = \\lambda'(\\tp{u}u') = 0\n",
        "\\end{array}\n",
        "$$\n",
        "- Jeśli wektory własne $u,u'$ macierzy $A\\th{A}$ są prostopadłe, to ich odpowiedniki $v\\eqd\\th{A}u$ oraz $v'\\eqd\\th{A}u'$  są też prostopadłe:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "\\left[A\\th{A}u = \\lambda u,\\ v\\eqd\\th{A}u,\\ A\\th{A}u' = \\lambda' u',\\ v'\\eqd\\th{A}u',\\ \\th{u}u'=0\\right]\\\\\n",
        "\\lra \\th{v}v'=\\th{u}(A\\th{A}u') = \\lambda'(\\th{u}u') = 0\n",
        "\\end{array}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### M23L: Zadanie EVD-S\n",
        "\n",
        "Na przykładzie losowej symetrycznej macierzy $A\\inm{4}{4}$ porównaj symboliczne rozwiązanie problemu własnego z rozwiązaniem numerycznym."
      ],
      "metadata": {
        "id": "6ztfKR4Ppg5h"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import sympy as sym\n",
        "import numpy as np\n",
        "\n",
        "x = sym.Symbol('x')\n",
        "\n",
        "# Losowa macierz 4x4 z ograniczonym zakresem\n",
        "A = np.random.randint(-10, 10, size=(4, 4))\n",
        "A = (A + A.T) / 2\n",
        "As = sym.Matrix(A)\n",
        "\n",
        "xI = sym.diag(x, x, x, x)\n",
        "ex = sym.expand(xI - As)\n",
        "d = ex.det()\n",
        "p = sym.expand(d)\n",
        "\n",
        "# Rozwiązanie symboliczne\n",
        "roots = sym.nroots(p, n=9)\n",
        "print(\"Rozwiązanie symboliczne:\", roots)\n",
        "\n",
        "# Rozwiązanie numeryczne\n",
        "ls, vs = np.linalg.eig(A)\n",
        "print(\"Rozwiązanie numeryczne:\", ls)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Tp3pyVRQeaoo",
        "outputId": "42395b5e-27af-4c3e-e67e-6afa060dbd5f"
      },
      "execution_count": 92,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Rozwiązanie symboliczne: [-14.8201021, -5.97058318, 3.32565586, 8.46502942]\n",
            "Rozwiązanie numeryczne: [-14.82010211   8.46502942   3.32565586  -5.97058318]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eT8QTDan_smR"
      },
      "source": [
        "### SVD (Singular Value Decomposition) $-$ dekompozycja na macierze singularne\n",
        "\n",
        "*Twierdzenie SVD*\n",
        "\n",
        "Jeśli $A\\inm{m}{n}$, to istnieją znormalizowane macierze ortogonalne $U\\inm{m}{m}$, $V\\inm{n}{n}$, oraz macierz diagonalna $\\Sigma\\inm{m}{n}$ takie, że:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "A \\overset{svd}{=} U\\Sigma\\tp{V} = U_r\\Sigma_r\\tp{V}_r = \\ds\\sum_{i=1}^r\\sigma_iu_i\\tp{v}_i\\\\\n",
        "\\tp{A} \\overset{svd}{=} V\\tp{\\Sigma}\\tp{U} = V_r\\Sigma_r\\tp{U}_r =  \\ds\\sum_{j=1}^r\\sigma_jv_j\\tp{u}_j\\\\\n",
        "U\\inm{m}{m},\\ \\tp{U}U=I_m,\\ V\\inm{n}{n},\\ \\tp{V}V=I_n\\ -\\ \\text{warunki znormalizowanej ortogonalności}\\\\\n",
        "U=[\\underbrace{u_1,\\dots,u_r}_{U_r},\\underbrace{u_{r+1},\\dots,u_m}_{U_{r:}}],\\ u_i\\inv{m},\\ i=1,\\dots,m\\\\\n",
        "V=[\\underbrace{v_1,\\dots,v_r}_{V_r},\\underbrace{v_{r+1},\\dots,v_n}_{V_{r:}}],\\ v_j\\inv{m},\\ j=1,\\dots,n\\\\\n",
        "\\lra U=[U_r,U_{r:}],\\ V=[V_r,V_{r:}]\\\\\n",
        "r=\\rank{A}, \\sigma_1\\geq\\sigma_2\\geq\\dots\\geq\\sigma_r>0\\ -\\ \\text{dodatnie wartości singularne}\\\\\n",
        "\\Sigma_r \\eqd \\diag{\\sigma_1,\\dots,\\sigma_r}\n",
        "\\lra\n",
        "\\Sigma = \\left[\n",
        "  \\begin{array}{cc}\n",
        "  \\Sigma_r & 0_{r\\times(n-r)}\\\\\n",
        "  0_{(m-r)\\times r} & 0_{(m-r)\\times(n-r)}\n",
        "  \\end{array}\n",
        "\\right]\n",
        "\\end{array}\n",
        "$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sODXbV0Iyggx"
      },
      "source": [
        "*Uwaga:* Jesli zastąpimy transpozycję macierz przez transpozycję ze sprzężeniem, to powyższe twierdzenie SVD jest prawdziwe również dla macierzy zespolonych $A\\inmc{m}{n}$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dACj1xdPBr30"
      },
      "source": [
        "*Dowód twierdzenia SVD* na podstawie dekompozycji na macierze wektorów własnych.\n",
        "\n",
        "1. Realizujemy dekompozycje EVD dla macierzy $A\\tp{A}$:\n",
        "$$\n",
        "r = \\rank{A},\\ A\\tp{A} \\overset{evd}{=} U\\Lambda\\tp{U} \\lra A\\tp{A}U = U\\Lambda\n",
        "\\lra A\\tp{A}U_r = U_r\\Lambda_r\n",
        "$$\n",
        "1. Na podstawie znormalizowanej ortogonalnej bazy własnej $U_r$, uzyskanej dla macierzy $A\\tp{A}$ dobieramy najpierw ortogonalną bazę własną dla macierzy $\\tp{A}A$, a następnie ją normalizujemy:\n",
        "$$\n",
        "V'_r\\eqd \\tp{A}U_r,\\ V_r\\eqd V_r'\\Lambda_r^{-\\frac{1}{2}}=\\tp{A}U_r\\Lambda_r^{-\\frac{1}{2}} \\lra \\tp{A}AV_r = V_r\\Lambda_r\n",
        "$$\n",
        "1. *Wniosek*: $\\tp{U}_rAV_r=\\Lambda_r^{\\frac{1}{2}}$\n",
        "$$\n",
        "V_r\\eqd \\tp{A}U_r\\Lambda_r^{-\\frac{1}{2}} \\lra  \\tp{U}_rAV_r = \\tp{U}_rA\\tp{A}U_r\\Lambda_r^{-\\frac{1}{2}} = \\tp{U}_rU_r\\Lambda_r\\Lambda_r^{-\\frac{1}{2}}=\\Lambda_r^{\\frac{1}{2}}\n",
        "$$\n",
        "1. Bazy $U_r$, $V_r$ dla podprzestrzeni uzupełniamy do pełnych znormalizowanych i ortogonalnych baz. W przypadku $U_{r:}$ jest już znane z EVD, a w przypadku $V_{r:}$ istnienie takiego rozszerzenia wynika np. z konstrukcji Grama-Schmidta:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "\\underbrace{\n",
        "\\left[\n",
        "\\begin{array}{c}\n",
        "\\tp{U}_r\\\\\n",
        "\\tp{U}_{r:}\n",
        "\\end{array}\n",
        "\\right]\n",
        "}_{\\tp{U}}\n",
        "A\n",
        "\\begin{array}{l}\n",
        "\\underbrace{\n",
        "\\left[\n",
        "V_r,\\ V_{r:}\n",
        "\\right]\n",
        "}_{V}\n",
        "\\end{array}\n",
        "\\end{array}\n",
        "=\n",
        "\\left[\n",
        "\\begin{array}{cc}\n",
        "\\tp{U}_rAV_r & \\tp{U}_rAV_{r:}\\\\\n",
        "\\tp{U}_{r:}AV_r & \\tp{U}_{r:}AV_{r:}\n",
        "\\end{array}\n",
        "\\right]\n",
        "=\n",
        "\\underbrace{\n",
        "\\left[\n",
        "\\begin{array}{cc}\n",
        "\\Lambda_r^{\\frac{1}{2}} & 0_{r\\times(n-r)}\\\\\n",
        "0_{(m-r)\\times r} & 0_{(m-r)\\times(n-r)}\n",
        "\\end{array}\n",
        "\\right]\n",
        "}_{\\Sigma}\n",
        "$$\n",
        "1. Otrzymaliśmy równanie macierzowe typu $\\tp{U}AV=\\Sigma$, które mnożymy lewostronnie i prawostronnie przez odwrotności macierzy ortogonalnych $\\left(\\tp{U}\\right)^{-1}=U$ oraz $V^{-1}=\\tp{V}:$\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "\\tp{U}AV=\\Sigma \\lra \\left(\\tp{U}\\right)^{-1}\\tp{U}AVV^{-1}=\\left(\\tp{U}\\right)^{-1}\\Sigma V^{-1}\\\\\n",
        " \\lra A \\overset{svd}{=} U\\Sigma\\tp{V},\\ \\text{ gdzie}\\\\\n",
        "U \\eqd \\left[\n",
        "U_r,\\ U_{r:}\n",
        "\\right]\n",
        ",\\\n",
        "\\Sigma \\eqd\n",
        "\\left[\n",
        "\\begin{array}{cc}\n",
        "\\Lambda_r^{\\frac{1}{2}} & 0_{r\\times(n-r)}\\\\\n",
        "0_{(m-r)\\times r} & 0_{(m-r)\\times(n-r)}\n",
        "\\end{array}\n",
        "\\right],\\\n",
        "V \\eqd \\left[\n",
        "V_r,\\ V_{r:}\n",
        "\\right]\n",
        "\\end{array}\n",
        "$$\n",
        "1. Gdy wykonamy mnożenie blokowe w wyrażeniu $U\\Sigma\\tp{V}$, podstawiając $\\Sigma_r\\eqd\\Lambda_r^{\\frac{1}{2}}$, otrzymamy zredukowaną postać SVD:\n",
        "$$\n",
        "A = U\\Sigma\\tp{V} = \\left[\n",
        "U_r,\\ U_{r:}\n",
        "\\right]\n",
        "\\left[\n",
        "\\begin{array}{cc}\n",
        "\\Sigma_r & 0_{r\\times(n-r)}\\\\\n",
        "0_{(m-r)\\times r} & 0_{(m-r)\\times(n-r)}\n",
        "\\end{array}\n",
        "\\right]\n",
        "\\left[\n",
        "\\begin{array}{c}\n",
        "\\tp{V}_r\\\\\n",
        "\\tp{V}_{r:}\n",
        "\\end{array}\n",
        "\\right]\n",
        "= U_r\\Sigma_r\\tp{V}_r\n",
        "$$\n",
        "1. Gdy na postaci zredukowanej dokonamy rozwinięcia (na podstawie ogólnego twierdzenia o iloczynach zewnętrznych macierzy i ich kolumn) oraz uwzględnimy oznaczenie $\\sigma_i\\eqd\\sqrt{\\lambda}_i$, to otrzymamy:\n",
        "$$\n",
        "A \\overset{svd}{=} U\\Sigma\\tp{V} = U_r\\Sigma_r\\tp{V}_r = \\sum_{i=1}^r\\sigma_iu_i\\tp{v}_i\\quad\\Box\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7_oDSGvKrn6e"
      },
      "source": [
        "*Wnioski* (SVD macierzy danych a charakterystyczne liniowe bazy)\n",
        "1. $U=[u_1,\\dots,u_m]$ jest znormalizowaną ortogonalną bazą w przestrzeni $\\bb{R}^m$, tj.\n",
        "$$ x\\inv{m}\\lra x = U(\\tp{U}x) = \\ds\\sum_{i=1}^m(\\tp{u_i}x)u_i$$\n",
        "1. $U_r=[u_1,\\dots,u_r]$ jest znormalizowaną ortogonalną bazą w podprzestrzeni $\\span{A}\\subseteq\\bb{R}^m$.\n",
        "1. $U_{r:}=[u_{r+1},\\dots,u_m]$ jest znormalizowaną ortogonalną bazą w podprzestrzeni $\\ker{[\\tp{A}]}\\subset\\bb{R}^{m}$.\n",
        "1. $V=[v_1,\\dots,v_n]$ jest znormalizowaną ortogonalną bazą w przestrzeni $\\bb{R}^n$, tj.\n",
        "$$ x\\inv{n}\\lra x = V(\\tp{V}x) = \\ds\\sum_{j=1}^n(\\tp{v_j}x)v_j$$\n",
        "1. $V_r=[v_1,\\dots,v_r]$ jest znormalizowaną ortogonalną bazą w podprzestrzeni $\\span{\\tp{A}}\\subseteq\\bb{R}^n$.\n",
        "1. $V_{r:}=[v_{r+1},\\dots,v_n]$ jest znormalizowaną ortogonalną bazą w podprzestrzeni $\\ker{[A]}\\subset\\bb{R}^n$.\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WHro_kVZLa4K",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "33d70c5c-801c-4430-8ee1-eccbdb94f872"
      },
      "source": [
        "s2,U = eig(np.dot(A,A.T))\n",
        "s3,V = eig(np.dot(A.T,A))\n",
        "print(s2,s3,'\\n',U,'\\n',V)\n",
        "U2,s,V2t = svd(A)\n",
        "print(s*s,'\\n',U2,'\\n',V2t.T)"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[213.8312464   6.8032865  64.3654671] [213.8312464   6.8032865  64.3654671] \n",
            " [[ 0.17517411  0.977064   -0.12107839]\n",
            " [ 0.29357441  0.06554757  0.95368631]\n",
            " [-0.93974895  0.20260667  0.27535875]] \n",
            " [[-0.51818244 -0.81781672 -0.25032533]\n",
            " [-0.43769958  0.00212276  0.89911877]\n",
            " [ 0.73478298 -0.57547486  0.35905802]]\n",
            "[213.8312464  64.3654671   6.8032865] \n",
            " [[ 0.17517411  0.12107839 -0.977064  ]\n",
            " [ 0.29357441 -0.95368631 -0.06554757]\n",
            " [-0.93974895 -0.27535875 -0.20260667]] \n",
            " [[-0.51818244  0.25032533 -0.81781672]\n",
            " [-0.43769958 -0.89911877  0.00212276]\n",
            " [ 0.73478298 -0.35905802 -0.57547486]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yQjmkXrwqXxH"
      },
      "source": [
        "###Dekompozycja Choleskiego\n",
        "\n",
        "Jeśli $B\\inm{n}{n}$ jest macierzą symetryczną, nieujemnie określoną, to istnieje dolno-trókątna macierz $X\\inm{n}{n}$ taka, że $B=X\\tp{X}.$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bxWLY0puqqjy",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "37bccf20-ba6d-49f9-fbba-2cd318fb5a3a"
      },
      "source": [
        "B = np.dot(A,A.T)\n",
        "X = cholesky(B)\n",
        "print(B,'\\n',X,'\\n',np.dot(X,X.T))"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 14   4 -36]\n",
            " [  4  77 -42]\n",
            " [-36 -42 194]] \n",
            " [[ 3.74165739  0.          0.        ]\n",
            " [ 1.06904497  8.70960061  0.        ]\n",
            " [-9.62140471 -3.64130195  9.3898611 ]] \n",
            " [[ 14.   4. -36.]\n",
            " [  4.  77. -42.]\n",
            " [-36. -42. 194.]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rASiB5c5lgXp"
      },
      "source": [
        "###Uogólniony problem Procrustesa\n",
        "\n",
        "Uogólniony problem Procrustesa dotyczy optymalnego dopasowania źródłowej bryły 3D do bryły docelowej, gdy w procesie dopasowania dysponujemy dowolnymi transformacjami afinicznymi.\n",
        "\n",
        "W zastosowaniach związanych z identyfikacją modelu 3D, najczęściej zamiast bryły 3D, mamy zbiór punktów 3D, tzw. chmurę punktów.\n",
        "\n",
        "W terminach matematycznych możemy ten problem wyrazić jako klasyczny problem najmniejszych kwadratów:\n",
        "- mając źródłową chmurę punktów $C\\inm{3}{N}$ i docelową chmurę punktów $C'\\inm{3}{N}$,\n",
        "- mamy ustalić macierz obrotu $R\\inm{3}{3}$, współczynnik skalowania $s\\in\\bb{R}_+$ oraz wektor przesunięcia $b\\inv{3}$, tak by\n",
        "- zminimalizować miarę niedopasowania:\n",
        "$$\\boxed{\\cl{E}(R,s,b) \\eqd \\left\\|sRC+b\\tp{1}_N-C'\\right\\|^2}$$.\n",
        "\n",
        "Powyższy problem optymlizacyjny nazywamy *uogólnionym problemem Procustesa*. Problem ten ma ogólne rozwiązanie oparte na dekompozycji SVD dla macierzy $\\ov{C}'\\tp{\\ov{C}}\\overset{svd}{=}U\\Sigma\\tp{V}$, gdzie $\\ov{C}\\eqd C-\\ov{c}\\tp{1}_N$ jest operatorem centrowania chmury punktów, tj. takiej relokacji chmury punktów, by jej centroid $\\ov{c}\\eqd\\ds\\frac{1}{N}\\sum_{i\\in[N]}c_i,$ znalazł się w początku układu współrzędnych.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e9suVZXCuMKe"
      },
      "source": [
        "###Rozwiązanie ogólnego problemu Procrustesa\n",
        "\n",
        "Następujące formuły podaja jedyne optymalne rozwiązanie problemu Procrustesa, gdy chmury punktów $C, C'\\inm{M}{N}$ leżą w przestrzeni dowolnego wymiaru $M$:\n",
        "$$\n",
        "\\ov{C}'\\tp{\\ov{C}}\\overset{svd}{=}U\\Sigma\\tp{V}\n",
        "\\ds \\longrightarrow \\left[\\ds\n",
        "R^{\\ast} \\ass U\\tp{V},\\  s^{\\ast} \\ass \\ds\n",
        "\\frac{\\tr{\\Sigma}}{\\tr{\\ov{C}\\tp{\\ov{C}}}}=\n",
        "\\frac{\\ds\\sum_{i\\in[M]}\\sigma_i}{\\|\\ov{C}\\|^2},\\ b^{\\ast} \\ass \\ov{c}'-s^{\\ast}R^{\\ast}\\ov{c}\\right]\n",
        "$$\n",
        "\n",
        "*Uwaga:* Wzory na optymalne parametry się nie zmienią, gdy rozważymy macierz kowariancji $Cov(C',C)\\eqd\\frac{1}{N}\\ov{C}'\\tp{\\ov{C}}$ między zbiorami punktów $C'$ oraz $C$, bo w efekcie $\\sigma_i'=\\frac{\\sigma_i}{N}$, a jednocześnie stała działa normalizacyjnie na  obie chmury jako mnożenie przez czynnik $\\frac{1}{\\sqrt{N}}$:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "Cov(C',C)\\eqd\\frac{1}{N}\\ov{C}'\\tp{\\ov{C}} = \\left(\\frac{\\ov{C}'}{\\sqrt{N}}\\right)\\left(\\frac{\\tp{\\ov{C}}}{\\sqrt{N}}\\right)\n",
        "\\overset{svd}{=}U\\underbrace{\\left(\\frac{\\Sigma}{N}\\right)}_{\\Sigma'}\\tp{V}\\\\\n",
        "\\ds \\longrightarrow \\left[\\ds\n",
        "R^{\\ast} \\ass U\\tp{V},\\  s^{\\ast} \\ass\n",
        "\\frac{\\tr{\\Sigma'}}{\\tr{Cov(C,C)}}=\n",
        "\\ds\\frac{\\ds\\sum_{i\\in[M]}\\frac{1}{N}\\sigma_i}{\\left\\|\\frac{1}{\\sqrt{N}}\\ov{C}\\right\\|^2}=\\frac{\\ds\\sum_{i\\in[M]}\\sigma_i}{\\|\\ov{C}\\|^2}\\right]\n",
        "\\end{array}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KTkNnvqc_MLq"
      },
      "source": [
        "####Interpretacja statystyczna rozwiązania Procrustes\n",
        "1. Znaleziony optymalny ruch afiniczny transformuje punkt średni (centroid) źródłowego zbioru punktów na punkt średni (centroid) docelowego  zbioru punktów.\n",
        "1. Optymalne skalowanie $s^{\\ast}$ jest ilorazem kowariancji wzajemnej między zbiorem źródłowym a zbiorem docelowym a wariancją zbioru źródłowego.\n",
        "1. Optymalna macierz obrotu $R^{\\ast}$ mapuje wierszową singularną bazę $V$ macierzy kowariancji $cov(C',C)$ na jej kolumnową singularną bazę $U$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g1m21LcG1Ip_"
      },
      "source": [
        "#### M23L: Pytanie (o Procrustesa)\n",
        "\n",
        "Kim był Procrustes? Czy transformacje, których się imał były afiniczne?"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "Prokrustes, znany również jako Prokrustes z Attiki, był postacią z mitologii greckiej. Był on właścicielem gospody, w której gościł podróżnych. Jednak jego gościnność była podstępna. Prokrustes miał zwyczaj poddawać swoich gości brutalnym przekształceniom, aby dopasować ich do swoich łóżek. Miał dwie łóżka: jedno zbyt krótkie dla wysokich gości i drugie zbyt długie dla niskich gości. Aby dostosować długość ciała podróżnych do rozmiaru łóżka, Prokrustes stosował przemoc, cięcie lub rozciąganie ciał.\n",
        "\n",
        "Transformacje, których Prokrustes używał, nie były afiniczne. Transformacje afiniczne to takie przekształcenia, które zachowują proporcje i linie proste. Natomiast Prokrustes dokonywał przemocy i brutalnych manipulacji ciałem gości, co nie jest charakterystyczne dla przekształceń afinicznych. Były to przekształcenia nieliniowe, które miały na celu dopasowanie ciała gości do konkretnych rozmiarów jego łóżek, niezależnie od tego, czy byli zbyt wysocy czy zbyt niscy."
      ],
      "metadata": {
        "id": "hv-osUbQPKFX"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WMpFVUWOWAr9"
      },
      "source": [
        "#### Ćwiczenie (Funkcja kosztu - od normy do sumy kwadratów błędów dopasowania)\n",
        "\n",
        "Niech macierze punktów $C,C'\\inm{3}{N}$ będą w analizie funkcji strat zastąpione przez sekwencje punktów $c_i,c'_i$:\n",
        "$$C=[c_0,c_1,\\dots,c_{N-1}],\\quad C'=[c'_0,c'_1,\\dots,c'_{N-1}]$$\n",
        "Przedstaw funkcję strat $\\cl{E}(R,s,b) \\eqd \\left\\|sRC+b\\tp{1}_N-C'\\right\\|^2$ w postaci sumy:\n",
        "$$\n",
        "\\cl{E}(R,s,b) = \\sum_{i=0}^{N-1}\\|???-c'_i\\|^2\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8tyL1_F8uU8t"
      },
      "source": [
        "####Funkcja `generalProcrustes`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8zApHiZjyeRt"
      },
      "source": [
        "def cov_stats(X,Y): # X,Y shape: M x N, here M = 3\n",
        "  M,N = X.shape\n",
        "  xb = np.mean(X,axis=1); yb = np.mean(Y,axis=1)\n",
        "  Xb = (X.T-xb).T; Yb = (Y.T-yb).T\n",
        "  Cov_XX = np.dot(Xb,Xb.T)/N\n",
        "  Cov_XY = np.dot(Xb,Yb.T)/N\n",
        "  Cov_YX = np.dot(Yb,Xb.T)/N\n",
        "  Cov_YY = np.dot(Yb,Yb.T)/N\n",
        "  return xb, yb, Cov_XX, Cov_XY, Cov_YX, Cov_YY"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IOVuPCxHyeRu"
      },
      "source": [
        "def generalProcrustes(X,Y):\n",
        "  xb, yb, Cov_XX, Cov_XY, Cov_YX, Cov_YY = cov_stats(X,Y)\n",
        "  U,sigmas,Vt = svd(Cov_YX)\n",
        "  R_ast = np.dot(U,Vt)\n",
        "  s_ast = np.sum(sigmas)/np.trace(Cov_XX)\n",
        "  b_ast = yb-s_ast*np.dot(R_ast,xb)\n",
        "  return R_ast, s_ast, b_ast"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "afW3t4GIwMrB"
      },
      "source": [
        "#### M23L: Zadanie (testowanie funkcji `generalProcrustes`)\n",
        "\n",
        "Niech $C\\inm{3}{1000}$ będzie chmurą punktów o $n=1000$ punktów w przestrzeni 3D generowaną losowo według rozkładu Gaussa. Obracamy, skalujemy i przesuwamy tę chmurę zgodnie z formułami:\n",
        "$$\n",
        "c'_i =\n",
        "3\\cdot\n",
        "\\left[\n",
        "\\begin{array}{ccc}\n",
        "\\frac{\\sqrt{3}}{2} & -\\frac{1}{2} & 0\\\\\n",
        "\\frac{1}{2} & \\frac{\\sqrt{3}}{2} & 0\\\\\n",
        "0 & 0 & 1\n",
        "\\end{array}\n",
        "  \\right]\n",
        "  c_i +\n",
        "\\left[\n",
        "\\begin{array}{c}\n",
        "1\\\\\n",
        "0\\\\\n",
        "-1\n",
        "\\end{array}\n",
        "  \\right],\\ i=0,\\dots,999\n",
        "$$\n",
        "Na tych danych sprawdzono poprawność twierdzenia Procrustesa (patrz rozwiązanie poniżej).\n",
        "\n",
        "Twoim zadaniem jest sprawdzenie twierdzenia Procrustesa na następujących, innych danych:\n",
        "\n",
        "$$\n",
        "c'_i =\n",
        "3\\cdot\n",
        "\\left[\n",
        "\\begin{array}{ccc}\n",
        "\\frac{1}{2} & \\frac{\\sqrt{3}}{2} & 0\\\\\n",
        "-\\frac{\\sqrt{3}}{2} & \\frac{1}{2} & 0\\\\\n",
        "0 & 0 & 1\n",
        "\\end{array}\n",
        "  \\right]\n",
        "  c_i +\n",
        "\\left[\n",
        "\\begin{array}{c}\n",
        "-1\\\\\n",
        "0\\\\\n",
        "1\n",
        "\\end{array}\n",
        "  \\right],\\ i=0,\\dots, 999\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d8sWiY3qCRe9"
      },
      "source": [
        "**Rozwiązanie**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0cWy9m3fsjkR"
      },
      "source": [
        "*Definiujemy transformację afiniczną*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EsuKxSYF74s2"
      },
      "source": [
        "s = 3\n",
        "b = np.array([1,0,-1.])\n",
        "sq3 = np.sqrt(3)\n",
        "R = np.array([[sq3/2,-1/2,0],[1/2,sq3/2,0],[0,0,1]])"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "print(\"Sprawdzenie dla innych danych.\")\n",
        "s = 3\n",
        "b = np.array([-1, 0, 1])\n",
        "sq3 = np.sqrt(3)\n",
        "R = np.array([[1/2, sq3/2,0],\n",
        "              [-sq3/2, 1/2 ,0],\n",
        "              [0, 0, 1]])"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7ClP-uiJP4nb",
        "outputId": "f497184f-f31c-4780-efab-69d3efec63fe"
      },
      "execution_count": 32,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Sprawdzenie dla innych danych.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7kU5Kojds5wJ"
      },
      "source": [
        "*Sprawdzamy ortogonalność*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GpvNbSR6s-mL",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "9820935f-5766-418c-b683-b3b71d40a1c9"
      },
      "source": [
        "print('Orthogonal?:\\n',np.dot(R.T,R))"
      ],
      "execution_count": 34,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Orthogonal?:\n",
            " [[1. 0. 0.]\n",
            " [0. 1. 0.]\n",
            " [0. 0. 1.]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AHYqbiectHb9"
      },
      "source": [
        "*Generujemy losową chmurę i ją centrujemy*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3iuw6Jmz8nEB"
      },
      "source": [
        "from numpy.linalg import det,svd,norm\n",
        "from numpy.random import randn"
      ],
      "execution_count": 35,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8IX93UB589n9",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "07671c04-20cb-4aba-ce69-f8798f55b2dc"
      },
      "source": [
        "C = randn(3,100); c_avg = np.average(C,axis=1)\n",
        "C = C-np.outer(c_avg,np.ones(100)); c_avg = np.average(C,axis=1)\n",
        "print('Is centered?:',c_avg)"
      ],
      "execution_count": 36,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Is centered?: [-7.21644966e-18  4.32986980e-17  1.11022302e-17]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g-D4_ZiCt6q_"
      },
      "source": [
        "*Stosujemy transformację afiniczną*\n",
        "\n",
        "W efekcie mamy dwie chmury punktów: wejściową $C\\inm{3}{100}$ i wyjściową $C'\\inm{3}{100}$."
      ]
    },
    {
      "cell_type": "markdown",
      "source": [],
      "metadata": {
        "id": "uBqC4qSyzqct"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Lhc5forBtsKR"
      },
      "source": [
        "Cp = s*np.dot(R,C)+np.outer(b,np.ones(100))"
      ],
      "execution_count": 37,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iEf1cATjuaLj"
      },
      "source": [
        "*Implementujemy w sposób jawny rozwiązanie problemu Procrustesa*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "53vtU8xc9ZCI"
      },
      "source": [
        "cb = np.mean(C,axis=1); cpb = np.mean(Cp,axis=1)\n",
        "Cb = (C.T-cb).T; Cpb = (Cp.T-cpb).T\n",
        "U,Sigma,Vt = svd(np.dot(Cpb,Cb.T))\n",
        "R_star = np.dot(U,Vt)\n",
        "s_star = np.sum(Sigma)/np.sum(Cb*Cb)\n",
        "b_star = np.average(Cp,axis=1)"
      ],
      "execution_count": 38,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zWXOqn0Ju-7z"
      },
      "source": [
        "*Sprawdzamy zidentyfikowane parametry transformacji*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9zuRWhi-vGzT",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e169a91f-0743-4cb3-de86-fcc161c0b53c"
      },
      "source": [
        "print('R-R_star:',R-R_star)\n",
        "print('error for rotation matrix ||R-R_star||:',norm(R-R_star))\n",
        "print('error for scale = |s-s_star|:',abs(s-s_star))\n",
        "print('error for shift = ||b-b_star||:',norm(b-b_star))"
      ],
      "execution_count": 39,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "R-R_star: [[-4.44089210e-16  0.00000000e+00  2.71616187e-17]\n",
            " [ 5.55111512e-16  5.55111512e-17 -1.44384151e-16]\n",
            " [-4.81548234e-17  1.09298465e-17 -4.44089210e-16]]\n",
            "error for rotation matrix ||R-R_star||: 8.542152051096907e-16\n",
            "error for scale = |s-s_star|: 4.440892098500626e-16\n",
            "error for shift = ||b-b_star||: 1.666962068746921e-16\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XKZmpotxuqVa"
      },
      "source": [
        "*A teraz stosujemy funkcję `generalProcrustes`*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OAr7-qH8-PqJ"
      },
      "source": [
        "R_star,s_star,b_star = generalProcrustes(C,Cp)"
      ],
      "execution_count": 40,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YpnsyDo_w5F-"
      },
      "source": [
        "*i również sprawdzamy ...*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z1c3RWgowuN2",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "423db7d4-eace-4bfa-d41b-3d20ff4e966e"
      },
      "source": [
        "print('R-R_star:',R-R_star)\n",
        "print('error for rotation matrix ||R-R_star||:',norm(R-R_star))\n",
        "print('error for scale = |s-s_star|:',abs(s-s_star))\n",
        "print('error for shift = ||b-b_star||:',norm(b-b_star))"
      ],
      "execution_count": 41,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "R-R_star: [[-2.22044605e-16 -1.11022302e-16  2.35699940e-16]\n",
            " [ 0.00000000e+00 -1.11022302e-16 -6.78892305e-17]\n",
            " [-8.50796550e-17 -2.33745655e-17 -4.44089210e-16]]\n",
            "error for rotation matrix ||R-R_star||: 5.82339474391554e-16\n",
            "error for scale = |s-s_star|: 4.440892098500626e-16\n",
            "error for shift = ||b-b_star||: 1.1822950248236328e-16\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D7sSF_5CCYSW"
      },
      "source": [
        "**Koniec Rozwiązania**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9zPPbT1WXO3F"
      },
      "source": [
        "#### M23L: Pytanie (P.1)\n",
        "Tworzymy chmurę $C'=RC$ przez rotację o macierzy $R=[r_x,r_y,r_z]$, tj. $\\tp{R}R=I_3,$ $\\det{R}=1$. Modyfikujemy transformację zamieniając miejscami kolumny macierzy $R$ i  otrzymujemy macierz $R'\\eqd[r_y,r_z,r_x].$ Czy macierz $R'$ jest macierzą obrotu? <br>Uzasadnij swoją odpowiedź."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import math\n",
        "import numpy as np\n",
        "\n",
        "C = randn(3,100); c_avg = np.average(C,axis=1)\n",
        "\n",
        "theta = math.pi\n",
        "R = [[math.cos(theta), -math.sin(theta), 0],\n",
        "      [math.sin(theta),  math.cos(theta), 0],\n",
        "      [0, 0, 1]]\n",
        "\n",
        "R_prim = [[-math.sin(theta), 0,math.cos(theta)],\n",
        "    [math.cos(theta), 0, math.sin(theta)],\n",
        "    [0, 1, 0]]\n",
        "\n",
        "R = np.array(R)\n",
        "R_prim = np.array(R_prim)\n",
        "\n",
        "print(\"Czy R_prim to macierz jednostkowa? \\n\", np.allclose(np.dot(R_prim, R_prim.transpose()), np.eye(3)))\n",
        "print(\"Wyznacznik R_prim to: \", det(R_prim))\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0xaSB27JRc7O",
        "outputId": "07c735e8-788c-4b21-bba3-e86c075880af"
      },
      "execution_count": 73,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Czy R_prim to macierz jednostkowa? \n",
            " True\n",
            "Wyznacznik R_prim to:  1.0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Tak, ponieważ kolumny są zamieniane cyklicznie, co powoduje, że wyznacznik macierzy będzie równy 1 nawet po przesunięciu kolumn. Jest to również macierz ortogonalna, czyli R_prim spełnia obydwa zadane warunki na bycie macierzą obrotu."
      ],
      "metadata": {
        "id": "fVA3J4X2ZXiz"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NkXavR0-aoGr"
      },
      "source": [
        "#### M23L: Pytanie (P.2)\n",
        "\n",
        "Niech $C\\inm{3}{1000}$ będzie chmurą o $n=1000$ punktach w 3D losowanych z rozkładu Gaussowskiego. Wykonujemy obrót, skalowanie i przesunięcie według formuły:\n",
        "$$\n",
        "c'_i =\n",
        "3\\cdot\n",
        "\\left[\n",
        "\\begin{array}{ccc}\n",
        "\\frac{1}{2} & \\frac{\\sqrt{3}}{2} &  0\\\\\n",
        "\\frac{\\sqrt{3}}{2} & -\\frac{1}{2} & 0\\\\\n",
        "0 & 0 & 1\n",
        "\\end{array}\n",
        "  \\right]\n",
        "  c_i +\n",
        "\\left[\n",
        "\\begin{array}{c}\n",
        "1\\\\\n",
        "1\\\\\n",
        "1\n",
        "\\end{array}\n",
        "  \\right],\\ i=0,\\dots,999\n",
        "$$\n",
        "Zauważmy, że powyższa transformacja nie jest rotacją. Dlaczego? Zastosuj formuły Procrustesa dla chmur $C,C'$. Czy powyższe transformacje zostały odzyskane?"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "R = np.array([[1/2, sq3/2, 0],[sq3/2,-1/2,0],[0,0,1]])\n",
        "\n",
        "\n",
        "print(\"Czy podana macierz to macierz jednostkowa? \\n\", np.allclose(np.dot(R, R.transpose()), np.eye(3)))\n",
        "print(\"Wyznacznik macierzy to: \", det(R))\n",
        "print(\"\\nGdyby macierz była macierzą obrotu to obróciłaby chmurę punktów o kąty:\\n\", np.arccos(1/2), np.arccos(-1/2), np.arcsin(sq3/2), \"\\n Kąty obrotu nie są sobie równe, więc nie może być to macierz obrotu.\")\n",
        "print(\"\\nPonad to wyznacznik macierzy jest zbliżony do -1 , co oznacza, że nie jest to macierz obrotu, a macierz odbicia lustrzanego.\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nx1WgsKUaK3s",
        "outputId": "39658db3-c6ad-45fe-9293-d43d374adcf1"
      },
      "execution_count": 81,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Czy podana macierz to macierz jednostkowa? \n",
            " True\n",
            "Wyznacznik macierzy to:  -0.9999999999999999\n",
            "\n",
            "Gdyby macierz była macierzą obrotu to obróciłaby chmurę punktów o kąty:\n",
            " 1.0471975511965979 2.0943951023931957 1.0471975511965976 \n",
            " Kąty obrotu nie są sobie równe, więc nie może być to macierz obrotu.\n",
            "\n",
            "Ponad to wyznacznik macierzy jest zbliżony do -1 , co oznacza, że nie jest to macierz obrotu, a macierz odbicia lustrzanego.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from scipy.linalg import svd\n",
        "\n",
        "C = np.random.normal(size=(1000, 3))\n",
        "\n",
        "transform_matrix = np.array([\n",
        "    [0.5, np.sqrt(3)/2, 0],\n",
        "    [np.sqrt(3)/2, -0.5, 0],\n",
        "    [0, 0, 1]\n",
        "])\n",
        "\n",
        "C_prime = 3 * np.dot(C, transform_matrix) + np.array([1, 1, 1])\n",
        "\n",
        "centroid_C = np.mean(C, axis=0)\n",
        "centroid_C_prime = np.mean(C_prime, axis=0)\n",
        "\n",
        "C_centered = C - centroid_C\n",
        "C_prime_centered = C_prime - centroid_C_prime\n",
        "\n",
        "A = np.dot(C_centered.T, C_prime_centered)\n",
        "\n",
        "U, S, Vt = svd(A)\n",
        "\n",
        "R = np.dot(Vt.T, U.T)\n",
        "\n",
        "scale = np.sum(S) / np.sum(C_centered ** 2)\n",
        "\n",
        "t = centroid_C_prime - scale * np.dot(centroid_C, R)\n",
        "\n",
        "C_transformed = scale * np.dot(C, R) + t\n",
        "\n",
        "print(\"Macierz rotacji R:\")\n",
        "print(R)\n",
        "print(\"\\nSkalowanie:\")\n",
        "print(scale)\n",
        "print(\"\\nWektor przesunięcia t:\")\n",
        "print(t)\n",
        "\n",
        "difference = np.linalg.norm(C_transformed - C_prime)\n",
        "print(\"\\nBłąd dopasowania między C_transformed a C_prime:\")\n",
        "print(difference)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cKvhbDfSdCpJ",
        "outputId": "a900f7fb-5c9a-4419-c2e6-91f35a5a0328"
      },
      "execution_count": 91,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Macierz rotacji R:\n",
            "[[ 5.00000000e-01  8.66025404e-01  5.55192720e-17]\n",
            " [ 8.66025404e-01 -5.00000000e-01 -1.66131059e-16]\n",
            " [-3.33357193e-16  4.06294122e-17  1.00000000e+00]]\n",
            "\n",
            "Skalowanie:\n",
            "2.9999999999999996\n",
            "\n",
            "Wektor przesunięcia t:\n",
            "[1. 1. 1.]\n",
            "\n",
            "Błąd dopasowania między C_transformed a C_prime:\n",
            "1.1649377152629387e-13\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Metoda Procrustesa nie była w stanie dokładnie odzyskać początkowej transformacji, zgodnie z oczekwaniami, ponieważ początkowe przekształcenie zawierało nie tylko rotację i przesunięcie, ale także odbicie oraz nierównomierne skalowanie. Metoda Procrustesa ma na celu dopasowanie dwóch chmur punktów za pomocą jedynie rotacji, przesunięcia i jednolitego skalowania. Nie jest jednak w stanie odzwierciedlić odbicia ani uwzględnić nierównomiernego skalowania. Dlatego też odzyskana transformacja może różnić się od początkowej."
      ],
      "metadata": {
        "id": "21JQXKBYhAKu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "####*Dowód uogólnionego twierdzenia Procrustesa*\n",
        "\n",
        "Ten nietrywialny dowód opiera sie na kilku własnościach, z których dwie nie są oczywiste:\n",
        "1. *Własność centroidu* $c\\inv{M}$ zbioru $X=\\{x_1,\\dots,x_N\\}\\subset\\bb{R}^M:$<br>\n",
        "$$\n",
        "c = \\frac{1}{N}\\sum_{i=1}^Nx_i \\iff c = \\arg\\min_{a\\inv{M}}\\sum_{i=1}^N\\|x_i-a\\|^2_2\n",
        "$$<br>\n",
        "  *Dowód:* Niech $f_X(a)\\eqd\\ds\\sum_{i=1}^N\\|x_i-a\\|^2_2$. Wtedy\n",
        "$$\n",
        "f_X(a) = \\sum_{i=1}^N\\tp{x}_ix_i -2\\tp{a}\\sum_{i=1}^Nx_i+N\\cdot\\tp{a}a \\lra\n",
        "\\od{f_x(a)}{a} = 0_M -2\\sum_{i=1}^Nx_i+2N\\cdot a\n",
        "$$<br>\n",
        "  Minimalny punkt $c$ funkcji $f_X(a)$ musi być punktem stacjonarnym teh funkcji, tj. jej gradient musi się w nim zerować:\n",
        "$$\n",
        "\\left.\\od{f_x(a)}{a}\\right|_{a=c} = 0_M \\iff c = \\frac{1}{N}\\sum_{i=1}^Nx_i\\quad\n",
        "$$<br>\n",
        "  Co więcej $c$ jest jedynym punktem stacjonarnym, a więc nie tylko dostatecznym ale i warunkiem koniecznym na minimum funkcji $f_X(a).\\quad\\Box$\n",
        "1. *Symetria śladu złożenia macierzy $X,Y\\inm{M}{M}$:*\n",
        "$$\n",
        "X,Y\\inm{M}{M} \\lra \\tr{XY}=\\tr{YX}\n",
        "$$<br>\n",
        "  *Dowód:*\n",
        "  $$\n",
        "  \\begin{array}{l}\n",
        "  \\tr{XY} = \\ds\\sum_{i=1}^M(XY)_{ii}=\\sum_{i=1}^M\\sum_{k=1}^MX_{ik}Y_{ki} =\n",
        "  \\sum_{i=1}^M\\sum_{k=1}^MY_{ki}X_{ik}\n",
        "  \\\\[3pt]\\ds\\underbrace{=}_{i\\leftrightarrow k}\n",
        "  \\sum_{k=1}^M\\sum_{i=1}^MY_{ik}X_{ki} = \\sum_{i=1}^M\\sum_{k=1}^MY_{ik}X_{ki} =\n",
        "  \\sum_{i=1}^M(YX)_{ii} = \\tr{YX}\\quad\\Box\n",
        "  \\end{array}\n",
        "  $$\n",
        "1. *Liniowość:* $\\tr{(\\alpha X+\\beta Y)} =\n",
        "\\alpha\\tr{X}+\\beta\\tr{Y}$\n",
        "1. *Inwariantność względem transpozycji:* $X\\inm{M}{M} \\lra \\tr{\\tp{X}}=\\tr{X}$\n",
        "1. *Ślad iloczynu zewneętrznego jest iloczynem skalarnym*: $x,y\\inv{M}\\lra\\tr{x\\tp{y}} =\\tp{x}y$\n",
        "1. *Kwadrat normy macierzy to ślad iloczynu zewnętrznego tej macierzy:* $X\\inm{M}{N}\\lra\\|X\\|^2 = \\tr{X\\tp{X}}$:\n",
        "\n",
        "  *Dowód:*\n",
        "  $$\n",
        "  X=[x_1,\\dots,x_N], x_i\\inv{M} \\lra \\tr{X\\tp{X} }= \\tr{\\sum_{i=1}^Nx_i\\tp{x}_i} =\n",
        "  \\sum_{i=1}^N\\tr{x_i\\tp{x}_i} = \\sum_{i=1}^N\\tp{x}_ix_i = \\sum_{i=1}^N\\|x_i\\|^2 = \\|X\\|^2\n",
        "  $$"
      ],
      "metadata": {
        "id": "AQdjQ31N23A6"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H7xTDw4RlCdq"
      },
      "source": [
        "*Obserwacje wstępne:*<br>\n",
        "Problem minimalizacji funkcji błędu $\\cl{E}(R,s,b)\\eqd\\ds\\min_{s,R,b}\\left\\|sRC+b\\tp{1}_N-C'\\right\\|^2$ dla $C,C'\\inm{M}{N}$,\n",
        "możemy sprowadzić do problemu $\\min_{s,R}\\left\\|A-sRB\\right\\|^2$ w kilku krokach:\n",
        "1. Zauważmy, że przy ustalonym $s,R$ optymalny wektor $b$ minimalizuje błąd kwadratowy dla zbioru wektorów znajdujących się w kolumnach macierzy $C'-sRC$.\n",
        "$$\n",
        "b^{\\ast}=\\arg\\min_b\\|C'-sRC-b\\tp{1}_N\\|^2 = \\arg\\min_b\\sum_{i\\in[N]}\\|(c'_i-sRc_i)-b\\|^2 \\lra\n",
        "b^{\\ast}=\\frac{1}{N}\\sum_{i\\in[N]} (c'_i-sRc_i) = \\ov{c}'-sR\\ov{c}\n",
        "$$\n",
        "W takim razie optymalne $b$ musi być centroidem tego zbioru, tj. $b^{\\ast} = \\ov{c}'-sR\\ov{c}$.\n",
        "2. Zastępując w funkcji błędu $\\cl{E}$ parametr $b$ jego optymalną wartością względem $s,R$, definiujemy macierze $A,B$, jak następuje:\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "b = \\ov{c}'-sR\\ov{c},\\\n",
        "A \\eqd C'-\\ov{c}'\\tp{1}_N,\\\n",
        "B \\eqd C-\\ov{c}\\tp{1}_N\\\\[2pt]\n",
        "\\lra \\cl{E}(R,s,b) = \\left\\|sRC+b\\tp{1}_N-C'\\right\\|^2 = \\left\\|A-sRB\\right\\|^2\n",
        "\\end{array}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rCPyM7z6OpZz"
      },
      "source": [
        "*Zasadnicza część dowodu:*<br>\n",
        "1. Norma w postaci *trace*:\n",
        "$$\\|A-sRB\\|_F^2=\\tr{A\\tp{A}}+s^2\\cdot\\tr{B\\tp{B}}-2s\\cdot\\tr{A\\tp{B}\\tp{R}}$$\n",
        "1. Równoważny problem *trace* (przy ustalonym $s\\in\\bb{R}_+$):\n",
        "$$R^*=\\arg\\max_{R:\\tp{R}R=I_m}\\tr{A\\tp{B}\\tp{R}}$$\n",
        "1. Dekompozycja SVD: $$A\\tp{B}\\xeq{svd}{U\\Sigma \\tp{V}}$$\n",
        "1. Macierz ortogonalna $Z\\eqd \\tp{V}\\tp{R}U:$\n",
        "$$\\tr{A\\tp{B}\\tp{R}}=\\tr{U\\Sigma \\tp{V}\\tp{R}}=\\tr{\\Sigma \\tp{V}\\tp{R}U} =\\tr{\\Sigma Z}$$\n",
        "1. Górne ograniczenie:\n",
        "$$|Z_{ij}|\\leq 1 \\lra \\tr{\\Sigma Z}=\\sum_i\\sigma_i Z_{ii}\\leq \\sum_i\\sigma_i.$$\n",
        "1. *Konkluzja 1:*\n",
        "$R=U\\tp{V} \\lra Z=I_m$ maksymalizuje $\\tr{\\Sigma Z}= \\tr{A\\tp{B}\\tp{R}}$ i wtedy przy dowolnym $s\\in\\bb{R}_+$ wyrażenie  $\\|A-sRB\\|_F^2$ przyjmuje wartość minimalną równą:\n",
        "$$\\tr{A\\tp{A}}+s^2\\cdot\\tr{B\\tp{B}}-2s\\cdot\\tr{\\Sigma}$$\n",
        "1. *Konkluzja 2:* Minimalna wartość błędu modelowania $\\|A-sR^{\\ast}B\\|_F^2$ przyjmowana jest dla współczynnika skalowania $$s^{\\ast}=\\frac{\\tr{\\Sigma}}{\\tr{B\\tp{B}}}$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CPPAxNIS2uPj"
      },
      "source": [
        "### Geometryczna interpretacja dekompozycji SVD\n",
        "\n",
        "Wartości singularne $\\Sigma$ i macierze ortogonalne $U,V$ w dekompozycji SVD dla macierzy $A\\inm{m}{n}$ mają ciekawa interpretację geometryczną. Mianowicie rozważmy obraz $S'$ sfery jednostkowej $S(0_n;1)$ przy odwzorowaniu nieosobliwym $A:\\bb{R}^n\\ra\\bb{R}^n:$\n",
        "$$\n",
        "S'\\eqd\\{Ax|\\ x\\inv{n},\\ \\tp{x}x=1\\}\n",
        "$$\n",
        "Można łatwo sprawdzić, że zbiór $S'$ ma symetrię środkową, tj. $x\\in S'$ wtt $-x\\in S'.$\n",
        "Można też znaleźć punkty w $S'$ najbliższe i najdalsze względem punktu zerowego. Po prostu znajdujemy stacjonarne punkty funkcj Lagrange'a  $L(x,\\lambda)=\\|Ax\\|^2-\\lambda\\|x\\|^2=\\tp{x}(\\tp{A}A)x-\\lambda\\tp{x}x$:\n",
        "$$\n",
        "  \\od{L(x,\\lambda)}{x}=2\\tp{A}{A}x-2\\lambda x=0_n\\lra \\tp{A}Ax=\\lambda x \\text{ oraz } \\|Ax\\|^2=\\lambda\\|x\\|^2=\\lambda\n",
        "$$\n",
        "Opis wszystkich punktów stacjonarnych tej funkcji faktycznie prowadzi do EVD macierzy $\\tp{A}A$:\n",
        "$\\tp{A}A\\overset{evd}{=}V\\Sigma^2\\tp{V}$.\n",
        "\n",
        "Mamy więc:\n",
        "- punkt najdalszy\n",
        "$y_{max}=Av_1=U\\Sigma\\tp{V}v_1=U\\Sigma e_1 = U(\\sigma_1e_1)=\\sigma_1u_1$,\n",
        "- podobnie punkt najbliższy $y_{min}=Av_n=U\\Sigma\\tp{V}v_n=\\sigma_nu_n$\n",
        "- oraz punkty siodłowe $y_i=Av_i=U\\Sigma\\tp{V}v_i=\\sigma_iu_i$ dla $1<i<n$.\n",
        "\n",
        "W takim razie możemy scharakteryzowac kształt $S'$ będący obrazem sfery jednostkowej przy transformacji $A$ - jest hiper-elipsoidą, w której:\n",
        "- półosie mają kierunki zgodne (lub przeciwne) z lewostronnymi wektorami singularnymi macierzy $A$, tj. z wektorami  $u_i$, $i=1,\\dots,n$,\n",
        "- długości tych półosi sa równe  $\\sigma_i$, $i=1,\\dots,n$.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### M23L: Pytanie SVD-AX\n",
        "\n",
        "Czym są półosie hiperelipsoidy w terminach algebraicznych?"
      ],
      "metadata": {
        "id": "fxLD2Gc3rLUC"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Półosie hiperelipsoidy w terminach algebraicznych są opisane przy użyciu macierzy kowariancji. Hiperelipsoida jest wielowymiarowym odpowiednikiem elipsoidy w przestrzeni, gdzie każdy wymiar jest reprezentowany przez półosie elipsoidy.\n",
        "\n",
        "Jeśli mamy hiperelipsoidę w n-wymiarowej przestrzeni, to półosie hiperelipsoidy są związane z macierzą kowariancji (lub macierzą kowariancji) danej jako Σ (Sigma). Macierz kowariancji to symetryczna macierz kwadratowa o rozmiarze n x n, gdzie n jest liczbą wymiarów przestrzeni.\n",
        "\n",
        "Półosie hiperelipsoidy odpowiadają pierwiastkom (wartościom własnym) tej macierzy kowariancji. Pierwiastki te można przedstawić jako d1, d2, ..., dn, gdzie d1 >= d2 >= ... >= dn. Są to długości półosi hiperelipsoidy wzdłuż odpowiednich osi współrzędnych.\n",
        "\n",
        "Ogólnie, jeśli Σ to macierz kowariancji, to półosie hiperelipsoidy wzdłuż poszczególnych osi współrzędnych są d1, d2, ..., dn, gdzie d1, d2, ..., dn to pierwiastki (wartości własne) macierzy Σ. Ich wartości i wzajemne stosunki opisują, jak hiperelipsoida jest rozciągana i obracana w przestrzeni."
      ],
      "metadata": {
        "id": "_0O0lRJXRFoW"
      }
    }
  ]
}