numerics-2022/Task 3/Lecture 3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 3. Projectors. Least squares problem. QR factorization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Projectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* A projector is a square matrix $P$ which satisfies\n",
    "$$\n",
    "P^2=P\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Lets inspect the geometric meaning of a projection\n",
    "![](proj.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* For any vector $v$ in $\\mathrm{range}(P)$ one has $Pv=v$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* For $S_1 = \\mathrm{range}(P)$ and $S_2=\\mathrm{null}(P)$ it is said that $P$ projects onto $S_1$ along $S_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Orthogonal projectors\n",
    "* For orthogonal projector: $S_1\\perp S_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* One can show (homework) that this is equivalent to $P = P^T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Lets inspect the geometric meaning of an orthogonal projection\n",
    "![](projorth.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Projection with orthonormal basis\n",
    "* Let us consider orthogonal projection onto certain $S_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Let dimension of $S_1$ be $n$; $q_1,q_2,...,q_n$ form an orthonormal basis for $S_1$ and $q_{n+1},...,q_m$ form an orthonormal basis for $S_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Consider \n",
    "$$\n",
    "Q = \\bigg[q_1\\bigg|q_2\\bigg|...\\bigg|q_n\\bigg|q_{n+1}\\bigg|...\\bigg|q_m\\bigg]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Then:\n",
    "$$\n",
    "PQ = \\bigg[q_1\\bigg|q_2\\bigg|...\\bigg|q_n\\bigg|0\\bigg|...\\bigg|0\\bigg]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " * We thus have $Q^TPQ =\\mathrm{diag}(1,1,...,1,0,0,0)=\\Sigma$\n",
    "which provides SVD decomposition ot $P$: \n",
    "\n",
    "$$P = Q\\Sigma Q^T$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Reduced SVD of a projector $P$ reads:\n",
    "$$\n",
    "P = \\hat Q \\hat Q^T\n",
    "$$\n",
    "where $\\hat Q$ is not orthogonal matrix (as it is not square), but it has orthonormal columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Projection with an arbitrary basis \n",
    "* Suppose that we want to project on a space spanned by columns of arbitrary $A$ (they are not necessarily orthonormal, but assumed to be linearly independent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* In passing from $v$ to its orthogonal projection $y \\in \\mathrm{range}(A)$, the difference $y-v$ must be orthogonal to $\\mathrm{range}(A)$:\n",
    "$A^T (Ax-v) = 0$, indeed:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](projorth2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* If the columns of $A$ are linearly independent:\n",
    "$$\n",
    "x = (A^T A)^{-1}A^T v\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* The projection is thus given by the operator\n",
    "$$\n",
    "P = A(A^T A)^{-1}A^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Linear systems\n",
    "\n",
    "- Linear systems appear as:\n",
    "\n",
    "    - Linear regression problems\n",
    "    - Discretization of partial differential/integral equations\n",
    "    - At elementary steps of iterative optimization techniques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear equations and matrices\n",
    "A linear system of equations can be written in the form\n",
    "\n",
    "\\begin{align*}\n",
    "    &2 x + 3 y = 5\\quad &\\longrightarrow \\quad &2x + 3 y + 0 z = 5\\\\\n",
    "    &3 x + 2z = 4\\quad &\\longrightarrow\\quad &3 x + 0 y + 2 z = 4\\\\\n",
    "    &x + y = 2\\quad &\\longrightarrow\\quad  & 1 x + 1 y + 0 z = 2\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix form\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "2 & 3 & 0 \\\\\n",
    "3 & 0 & 2 \\\\\n",
    "1 & 1 & 0 \\\\\n",
    "\\end{pmatrix}\\begin{pmatrix}\n",
    "x \\\\\n",
    "y \\\\\n",
    "z \n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "5 \\\\\n",
    "4 \\\\\n",
    "2\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "or simply\n",
    "\n",
    "$$ A u = b  $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Over/under determined linear systems\n",
    "\n",
    "If the system $Ax = b$ has\n",
    "\n",
    "-  more equations than unknowns, it is called overdetermined system (typically no solution)\n",
    "\n",
    "- less equations than unknowns, it is called underdetermined system (typically many solutions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overdetermined linear systems\n",
    "\n",
    "- The number of equations is > than the number of unknowns. \n",
    "\n",
    "- The simplest example is linear fitting, fitting a set of 2D points by a line.\n",
    "\n",
    "- Then, a typical way is to minimize the residual (least squares)\n",
    "\n",
    "$$\\Vert A x - b \\Vert_2 \\rightarrow \\min$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Underdetermined linear systems\n",
    "\n",
    "- The number of equations is > than the number of unknowns. \n",
    "\n",
    "- In this case, there could be infinitely many solutions $x$\n",
    "\n",
    "- Typical problem in Machine Learning: employ some kind of regularization. For example, one may find the solution to $Ax=b$ with a minimal norm \n",
    "$$\\Vert x \\Vert_2\\to\\min$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example of least squares\n",
    "Consider a two-dimensional example. Suppose we have a linear model \n",
    "\n",
    "$$y = ax + b$$\n",
    "\n",
    "and noisy data $(x_1, y_1), \\dots (x_n, y_n)$. Then the linear system on coefficients will look as follows\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "a x_1 &+ b &= y_1 \\\\\n",
    "&\\vdots \\\\\n",
    "a x_n &+ b &= y_n \\\\\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "or in a matrix form\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "x_1 & 1 \\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "x_n & 1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "y_1 \\\\\n",
    "\\vdots  \\\\\n",
    "y_n \\\\\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "\n",
    "which represents overdetermined system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solving the least squares\n",
    "- Suppose we need to solve\n",
    "\n",
    "$$ \\Vert A x - b \\Vert_2 \\rightarrow \\min_x, $$\n",
    "\n",
    "where $A$ is $m \\times n$, $m \\geq n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](projorth3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- The residual $r = A x - b$ should be orthogonal to $\\mathrm{range} A$, and\n",
    "or equivalently:\n",
    "$$\n",
    "P b = A x,\n",
    "$$\n",
    "where $P$ is a projector onto the range of $A$, that is $P = A(A^T A)^{-1}A^T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- If $A$ is of a full rank, then \n",
    "$$\n",
    "x = (A^T A)^{-1} A^T b\n",
    "$$\n",
    "- The matrix $(A^T A)^{-1} A^T$ is known as a pseudoinverse, denoted as $A^+$. Can be easily computed in numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1fdc8b5e340>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD6CAYAAAC4RRw1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeNUlEQVR4nO3deXRV5b3/8feTEEhEMCjRC4kY2qvMQjBoLIJCvCAiELE4FH6KWnAoUkVTiBN4V2/hGqBqFf1ZxwqoVREQkaACalVQQpBUIQoImkEJ0SBDICE8948TQgKZzpQz7M9rra5lNufs8+1RPjw8+7u/21hrERGR0BMR6AJERMQzCnARkRClABcRCVEKcBGREKUAFxEJUQpwEZEQ1WiAG2OeM8bsMsb8u8axTGPMFmPMJmPMm8aYWL9WKSIiJzCN9YEbYwYC+4B/WGt7Vh0bAqyy1h42xvwvgLV2amMf1r59e5uYmOh10SIiTpKdnb3bWht3/PEWjb3RWvuhMSbxuGMra/y4FvhtU4pITExk/fr1TXmpiIhUMcbsrOu4L/bAbwLe8cF5RETEDV4FuDHmPuAwsKCB10w0xqw3xqwvLi725uNERKQGjwPcGHMDcAUw1jawkW6tfdpam2ytTY6LO2ELR0REPNToHnhdjDGXAVOBi621B7wpoKKigvz8fA4ePOjNacSPoqOjSUhIICoqKtCliEgNjQa4MeZl4BKgvTEmH5gOZACtgHeNMQBrrbW3elJAfn4+bdq0ITExkapzSRCx1lJSUkJ+fj6dO3cOdDkiUkNTulCuq+Pws74q4ODBgwrvIGaM4bTTTkPXL0Q8szingMysPApLy+gYG0P60C6kJcX75NwebaH4msI7uOnfj4hnFucUkLEol7KKSgAKSsvIWJQL4JMQ1630IiJ+kpmVVx3eR5VVVJKZleeT8yvARUT8pLC0zK3j7lKAi4j4ScfYGLeOu0sBXqWsrIyLL76YysrKxl/sgfLycgYOHMjhw4f9cn4RCT7pQ7sQExVZ61hMVCTpQ7v45PwhF+CLcwroP2sVnae9Tf9Zq1icU+CT8z733HOMHj2ayMjIxl/sgZYtW5Kamsqrr77ql/OLSPBJS4pn9hVn8nDrhbThAPGxMcwc3ctnXSghFeBHr+gWlJZhOXZF1xchvmDBAkaNGgXAoEGDePfddwG4//77mTx5cqPvz83NpX///tU/b9iwgcGDB9d6TVpaGgsW1Dt1QETCibWw6Z8M/2AEV9sV5N54Mh9PG+yz8IYgaSNsqoau6HrzpZSXl7N9+3aOjrp96KGHePDBB9m1axc5OTksXbq00XP06NGDbdu2UVlZSWRkJHfffTdz5syp9ZqePXvy+eefe1yniISIn3fC21Ng63sQnwwjH4Mzevj8Y0JqBe6vK7q7d+8mNja2+ueBAwdirWXu3Lm88sor1dsqDzzwQL3niIiIoEePHnz55Ze88cYbdOrUib59+9Z6X2RkJC1btmTv3r1e1SsiQaryMHzyOMxLge/WwrCH4eaVfglvCLEVeMfYGArqCGtvr+jGxMTUmsWSm5tLUVER7du3p02bNgD88MMPjV6ATElJ4eOPP2bevHmsWLGizvcdOnSI6Ohor+oVkSBUtAmW3gFFG+Gcy2D4HDglwa8fGVIrcH9d0W3Xrh2VlZUcPHiQoqIixo4dy5IlS2jdujVZWVkA5OTk0KdPn+r3pKamUlBQe+89JSWF+++/nyuvvJL4+PgT3ldSUkJcXJyGQomEk/ID8O6D8PQl8Esh/PZ5uO4Vv4c3hFiApyXFM3N0L+JjYzDg0yu6Q4YMYeXKlYwePZo5c+bQrVs3HnjgAWbMmAHAxo0bq4P4yJEjbN26lVNPPbXWObp27UqrVq2YOvXY0+Vqvm/16tVcfvnlXtcqIkFi22p48kL4+FFIGguTPoOeo6GZxk+E1BYKuELcl1dxj5o0aRJz587l008/rT42cODA6p+3bt3K2WefDcBXX33FVVddRUxM7a2bRx99lJkzZ9K6devqYzXft3DhQmbOnOnz2kWkmR34CbLuhS9ehlN/DTcsg84Dmr2MkFqB+1NSUhKDBg2q90aeZ599logI19fVs2dP5s6dW/1r27Zto2vXrpSVlXHDDTfU+b7y8nLS0tLo0sU3DfwiEgBVrYE8ngy5r8GAe+C2TwIS3hCCK3B/uummmzx6369//Wu2bNnS4GtatmzJ9ddf79H5RSQINFNroDsU4CIiDak8DOuegtX/AybC1RrY7/cQ4Z+7tt2hABeRsOSTBykEoDXQHQpwEQk7Xj9IofwAfDDLdVPOSae5WgN7XNls3SVNpQAXkbDj1diNbath2Z3w8w7oez38139DTDu/1eoNBbiIhB2Pxm4ESWugOxTgIhJ23Bq7Ya2rJXDFNDi4x9UaODAdooJ/5IX6wIGTTz7ZL+fdsWMHCxcurPPXjhw5wuTJk+nZsye9evWiX79+fPvtt36pQ8Rpmjx24+edsOC3sGgCtOsMt3wIqQ+ERHiDVuB+dTTAf/e7353wa6+++iqFhYVs2rSJiIgI8vPza93B6WvWWqy11TcjiYSzo/vc9XahBHFroDv0u7ke27Zt47LLLuO8885jwIAB1TfqvPXWW1xwwQUkJSVx6aWX8uOPPwLwwQcf0KdPH/r06UNSUhJ79+5l2rRpfPTRR/Tp04e//vWvtc5fVFREhw4dqgM1ISGBdu1cF0qef/55zjnnHC6++GImTJjApEmTABg/fjyvv/569TmO/s1h3759pKam0rdvX3r16sWSJUsA1x8g3bp14/bbb6dv3758//33ZGZm0q9fP84991ymT58OwP79+xk+fDi9e/emZ8+eemqQhIW0pHg+njaYb2cNr/0ghaJN8EwqrLwPOg+EP6yDC24JufCGYFuBvzMNfsj17Tn/oxcMm+X22yZOnMhTTz3F2Wefzbp167j99ttZtWoVF110EWvXrsUYwzPPPMPDDz/MnDlzmD17Nk888QT9+/dn3759REdHM2vWLGbPns2yZctOOP/VV1/NRRddxEcffURqairjxo0jKSmJoqIipk+fTnZ2NqeccgqDBg0iKSmpwVqjo6N58803adu2Lbt37yYlJYWRI0cCkJeXx/PPP8+8efNYuXIl33zzDZ999hnWWkaOHMmHH35IcXExHTt25O233wZgz549bn9fIkEvRFoD3RFcAR4k9u3bxyeffMKYMWOqjx06dAiA/Px8rrnmGoqKiigvL6dz584A9O/fnylTpjB27FhGjx5NQkLDzf4JCQnk5eWxatUqVq1aRWpqKq+99hp79+7lkksuIS4uDoBrrrmGr7/+usFzWWu59957+fDDD4mIiKCgoKD6bwZnnXUWKSkpAKxcuZKVK1dW/4Gwb98+vvnmGwYMGMA999zD1KlTueKKKxgwILivvIu4LYRaA90RXAHuwUrZH44cOUJsbCwbN2484dfuuOMOpkyZwsiRI1mzZk31uNlp06YxfPhwli9fTkpKCu+9916jn9OqVSuGDRvGsGHDOOOMM1i8eDGpqamYelYELVq04MiRI4ArtMvLywHX8zyLi4vJzs4mKiqKxMTE6gdU1NxXt9aSkZHBLbfccsK5s7OzWb58ORkZGQwZMoQHH3yw0fpFgl4Itga6Q3vgdWjbti2dO3fmtddeA1zB98UXXwCu7YWjD2t48cUXq9+zbds2evXqxdSpU0lOTmbLli20adOm3senbdiwgcLCQsD1B8amTZs466yzuOCCC1izZg0lJSVUVFRU1wCQmJhIdnY2AEuWLKGioqK6ptNPP52oqChWr17Nzp076/zMoUOH8txzz7Fv3z4ACgoK2LVrF4WFhZx00kmMGzeOe+65hw0bNnj83YkEhSCbGugvwbUCD5ADBw7U2vKYMmUKCxYs4LbbbuPPf/4zFRUVXHvttfTu3ZsZM2YwZswY4uPjSUlJqW79e+SRR1i9ejWRkZF0796dYcOGERERQYsWLejduzfjx4/nrrvuqv6MXbt2MWHChOqtmfPPP59JkyYRHR3NjBkzuPDCC+nQoQN9+/atHnE7YcIERo0axfnnn09
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "a0 = 1; b0 = 2; n = 10\n",
    "\n",
    "xi = np.arange(n)\n",
    "yi = a0 * xi + b0 + np.random.random(n)\n",
    "A = np.hstack([xi[:,None],  np.ones(n)[:,None]])\n",
    "coef = np.linalg.pinv(A) @ yi\n",
    "\n",
    "plt.plot(xi, yi, 'o', label='$(x_i, y_i)$')\n",
    "plt.plot(xi, coef[0]*xi + coef[1], label='Least Squares')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What happens inside pinv?\n",
    "We dont know yet -- computing pseudoinverse might be tough. Let us inspect a computationally cheaper way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solving the linear least squares via $QR$ decomposition.\n",
    "\n",
    "- Any matrix can be factored into a product \n",
    "$$ A = Q R$$\n",
    "where $Q$ is unitary, and $R$ is upper triangular.\n",
    "- Finding optimal $x$ is equivalent to solving\n",
    "$$\n",
    "A^T A x = A^T b\n",
    "$$\n",
    "If $A$ is of full rank, it is the same as\n",
    "$$ Rx = Q^T b. $$\n",
    "- Since $R$ is upper triangular, the solving of this linear system costs $\\mathcal{O}(n^2)$. \n",
    "- It is more stable than solving the normal equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$R^T Q^T Q R x = R^T Q^T b$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$R x = Q^T b$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# QR decomposition\n",
    "\n",
    "- The name implies that a matrix is represented as a product  \n",
    "$$\n",
    "    A = Q R, \n",
    "$$\n",
    "where $Q$ is a column orthonormal matrix and $R$ is upper triangular.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- The matrix sizes: $Q$ is $m \\times n$, $R$ is $n \\times n$ if $m\\geq n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- QR decomposition is defined for any rectangular matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- QR decomposition can be computed in a $O(mn^2)$ (for $m>n$) (via so called direct methods)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Full QR vs reduced QR\n",
    "![](full.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition: applications\n",
    "\n",
    "This decomposition plays a crucial role in many problems:\n",
    "- Computing orthonormal bases in a linear space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Solving overdetermined systems of linear equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- QR-algorithm for the computation of eigenvectors and eigenvalues ([one of the 10 most important algorithms of the 20th century](https://archive.siam.org/pdf/news/637.pdf)) is based on the QR decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Existence of QR decomposition\n",
    "Every rectangular $m \\times n$  matrix has a QR decomposition. \n",
    "\n",
    "\n",
    "There are several ways to prove it and compute it:\n",
    "\n",
    "- Geometrical: using the Gram-Schmidt orthogonalization\n",
    "- Practical: using Householder/Givens transformations "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gram-Schmidt orthogonalization\n",
    "\n",
    "- QR decomposition is a mathematical way of writing down the Gram-Schmidt orthogonalization process.  \n",
    "- Given a sequence of vectors $a_1, \\ldots, a_m$ we want to find orthogonal basis $q_1, \\ldots, q_m$ such that every $a_i$ is a linear combination of such vectors.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Gram-Schmidt:\n",
    "- $q_1 := a_1/\\Vert a_1\\Vert$\n",
    "- $q_2 := a_2 - (a_2, q_1) q_1, \\quad q_2 := q_2/\\Vert q_2 \\Vert$\n",
    "- $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \\quad q_3 := q_3/\\Vert q_3 \\Vert$\n",
    "\n",
    "And so on.\n",
    "\n",
    "Note that the transformation from $Q$ to $A$ has triangular structure, since from the $k$-th vector we subtract only the previous ones. It follows from the fact that the product of triangular matrices is a triangular matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR decomposition: the practical way\n",
    "\n",
    "* Gram-Schmidt orthogonalization is not very practical as just described for numerical issues (will see later). One typically does it differently: Householder rotations (implemented in numpy)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More on Pseudoinverse\n",
    "Recall that we introduced pseudoinverse as\n",
    "$$\n",
    "A^+ = (A^T A)^{-1} A^T\n",
    "$$\n",
    "- Matrix $A^T A$ can be singular in general case.\n",
    "- Therefore, we need to introduce the concept of pseudoinverse matrix  $A^{\\dagger}$.\n",
    "- The matrix $$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^T A)^{-1} A^T$$ is called Moore-Penrose pseudoinverse of the matrix $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Pseudoinverse and linear systems\n",
    "Consider a system of equations $Ax\\approx b$\n",
    "- $\\forall x$, we have $\\left\\|Ax-b\\right\\|_{2}\\geq \\left\\|Az-b\\right\\|_{2}$, where \n",
    "$z=A^{+}b$\n",
    "- If $A x = b$ is satisfiable, the vector $z=A^{+}b$ is a solution, and satisfies $\\|z\\|_{2}\\leq \\|x\\|_{2}$ for all other solutions\n",
    "\n",
    "As a result, the solution to the linear least squares problem (including least norm solutions) can formally be written as\n",
    "$$x = A^{\\dagger} b.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simple cases\n",
    "\n",
    "- If matrix $A$ is squared and non-singular we get standard inverse of $A$:\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^T A)^{-1} A^T = (A^T A)^{-1} A^T = A^{-1} (A^{T})^{-1} A^T = A^{-1}$$\n",
    "- If matrix $A$ has full column rank, then $A^T A$ is non-singular and we get \n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^T A)^{-1} A^T = (A^T A)^{-1} A^T. $$\n",
    "- Pseudoinverse is uniquely defined in all other cases, too! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Least squares via SVD \n",
    "Lets inspect SVD approach to over- and under- determined problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider solving $Ax \\approx b$ for $A$ of arbitrary shape and possibly not of a full rank (the system cen be under- or overdetermined).\n",
    "\n",
    "Let $A = U \\Sigma V^T$ be SVD decomposition of $A$, then we need to solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "U\\Sigma V^T x \\approx b\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $U$ and $V$ are orthogonal, we can just switch to new variables, $\\bar x = V^T x$ and $\\bar b = U^T b$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Sigma \\bar x \\approx \\bar b\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to now, its not important if the system is under- or overdetermined. Here it becomes important. Lets consider underdetermined case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\sigma_1 & 0 & 0\\\\\n",
    "0 & \\sigma_2 & 0 \\\\\n",
    "0 & 0 & 0 \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\bar x_1 \\\\\n",
    "\\bar x_2 \\\\\n",
    "\\bar x_3\n",
    "\\end{pmatrix} \\approx\n",
    "\\begin{pmatrix}\n",
    "\\bar b_1 \\\\\n",
    "\\bar b_2 \\\\\n",
    "\\bar b_3\n",
    "\\end{pmatrix},\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\bar x_1 = \\bar b_1 / \\sigma_1, \\bar x_2 = \\bar b_2 / \\sigma_2, \\bar x_3 = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Example: Signal recovery "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "with np.load('data.npz') as data:\n",
    "    A, C = data['A'], data[\"C\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(25, 60)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACsCAYAAABikvffAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAe4ElEQVR4nO2dW4xl11nn/9/e+9zq0l3d7W637TYxiSxBQIOBVsgo82Aug8JF4/AQiUggIyE1DyAFCTR4eMmAhJQHbg8gJM/EihEQFAkysUYIYgVGgReIEzwTByexsTt22+2+uLuqq+rUOWefvT8e6hiKXv+v65yu6tOzj/4/yequ1avWXmvttb+zff7fxdwdQgghmkd2tycghBDi9pABF0KIhiIDLoQQDUUGXAghGooMuBBCNBQZcCGEaCgHMuBm9kEz+7qZvWxmTxzWpIQQQuyP3a4fuJnlAL4B4D8DuADgiwA+4u7/FP1O2zrexfJ042fBZ4vRzlONCQDI2AAAouvdKaJ9Z/Ooa97XgrWwy5F1W3XwccN1sOZo76Pr1WSQGaYWwq4363PA5hatjxFd7zD2Yqb7x34/6jvDnGfZz1nmOytsHuHUyD9Ez160SXQvgjHYPILr3fBrV9395M3tRTC7aXgfgJfd/RUAMLM/BfAYgNCAd7GM78t+KP0HYoCz5SU6hrEN6nT4BZnR6vK+HrTP9FDO8iEwHE09DxsM+Rit6W9f3W0nbfnGdjCHtC8AgBn8qqJd2YeD94K9D/bNRiVpPPjD7p1036yKPoh4u+2k9yRaH923yDAE99T6g/R60f1n+xmdYza3Iud9yzFvZ/3HM5yLTouPyz60ACAn6wvOBbtP4d6zufX7vG8R3KdWuhYfpPcOADBO99ODZ/1zO3/0TdZ+kNfOBwC8vufnC5M2IYQQc+Agb+DsIy/5yDSzcwDOAUAX/K1aCCHE7BzkDfwCgAf3/HwGwJs3d3L3J939rLufbSH4X0whhBAzcxAD/kUAD5vZt5pZG8BPAnjmcKYlhBBiP277KxR3H5vZLwD4KwA5gKfc/au3+h3LMmS9Xtq+upK2haJi+plTH53OswUAxstcoIs8MkZH0/5Vl3/u5QMyRiCu1EUguhAxzfMjU/etenxuWUkEoTN834ptLkCxsbMh37e6nfZlnjAAYIFY5Xnan64DgBMdLdrjbDS9h0QWiJtsbIucDWbRwYPruaVngO3xbue0yQIxlu1nuG/j6fctWnM+IkJhJEzPsm/jYPPr9GvbrIwE1nR9NkxtEwBkW4EwSZwTmLAJAM5E0xm9oQ7yHTjc/S8A/MVBxhBCCHF7KBJTCCEaigy4EEI0FBlwIYRoKAf6DnxmzGAsapJFJHVW6RD1cjdpq1a4MMlEnnGPR5lVHa6YDI+QMZamFyAjISYjQYYAF+MinAXcBcFyNRk3mkPd5pOgYweaS9WeIcw/WHNOxEa2ZoALYZEgmJEg2HGPzzcfBqIieXJqIroCXKSN7lN0DpmAGO4Fi0COohoJkRibBwHBRjTBPBCb2Tqq1mwiJlt3NOfOOpkcu3kAWpvpA2EtvslR9KgNyRgD/qCx5WVBZDMCzVRv4EII0VBkwIUQoqHIgAshREORARdCiIYiAy6EEA1lzl4ooLl87UjqcVKtpt4mAPc4KVf4Msrl1L2hDDxIRqtB+1raNl4OckQzsTlMus+bszL9h7oIrsfSTAd3lHkKWFRTYsQnR71FovoD+fQeOZEHgdXpL3jGJ111SRh0IOjT680YwmwVmRtbM2JvEUYUrs72fpb1ReeCjhvcj2Jres8iC2LpqedTlIo98rJhHlWBVw97rvPAo6M4mg4ceSEVO3yTWjfSSUeOZTZKbVmU8gAbvFlv4EII0VBkwIUQoqHIgAshREORARdCiIYyXxHTHWCFagl1l4eqViS0decEXwYTQcrl6cVKgItjkdxVc92VE4hx46UZqmiT9UUiHxPd6m5UcTu43A7J8R0qNCR8fBiEJQfiH1t4Ng7C1YmIVfWivWBtUdx20Ez2jomuAFAtMVWRjztT8vBINGXtwbkAWXe0FyXPbkHF6ZyclYhoyVkg0rK9i84he95bW7xvQeoXt7ZnzVVO0noEW28kh3lYwDxAb+BCCNFQZMCFEKKhyIALIURDkQEXQoiGciAR08zOA9gEUAEYu/vZw5iUEEKI/TkML5Tvd/erU/XMclqBvjpBQuk7XFbeOZl6p7BiBQAwOpJKxcMTM3h/AKhWUqW4fYzH4o5H6UQ6Pe51MxpyL5ssim8nsJ558PutduqmUY747a9r/j9mNUlDYIF3gxPXgshDJsKJV8d4FPxPI2uOLsc8CIJ9ywruLsLmVrQDt4kh2efgehZ4lrQ76f0bl/zg50U6jyxKQVClG1cOAq8u0hcAslZ6vfoY7QofE0+mko9Lirbvwu5f5AFECnV4EXlDsfQIwRyCLy/yQTrpaonvZ75J0orkM1R0CWchhBDi/3sOasAdwOfM7Etmdu4wJiSEEGI6DvoVygfc/U0zOwXgWTP7mrt/YW+HiWE/BwDdPIoEEEIIMSsHegN39zcnf14G8BkA7yN9nnT3s+5+tp31DnI5IYQQe7htA25my2a2+s7fAfwwgBcOa2JCCCFuzUG+QrkXwGfM7J1x/sTd//LWv+JARZT6mii3vUCZ7qRKcRkUY2C5GyJvk/oe7i3SW0lzE9y3doP2ZR4kdZDoYTAOPEBI/xM9kqQhYGc8vXfLOPA2iRhVqUIepYQosvSeRteLPG9Y+/U+/7841jcjcwCAnHhkLLdHM82NUQXr6+SpB8n1AV/H0S73cMpJlYV+SfJuADjSSccoyb0D+HmpgjMbeTix/lHfrWE65/6Ar4N5MgHccya612U/HZvfacCNJRcK9iIoCjFeIl49y/xctDZSG+DjoDJFwG0bcHd/BcB33e7vCyGEOBhyIxRCiIYiAy6EEA1FBlwIIRrKnKvSZ0AvrXowXu0kbYNjXHQpV4iImUbn745xKhVMvcfFjuUjXJU4s7aetH3PsddpXyZ47VRcVBzWvJ2Py+d8pEjn3ApKlffrVMxZyricUwfS5IWdND466jus0qO11t6hfY+3tmn7G4O1tG+XV81YaU2fCL+Xp4J1HZRAj0To4+10zlHfK6P0gN6/zMuM39Pm1QbY2Kfam7RvRdbydrlM+14cHE3algt+LjZKvvesfzsoE39lkO7Fxe0jtO9KICyzvYjE5mtE9L6ec4MxrtJnJB8F59t4u+fp3nc2gvQdq+n12qv8PuEt3qw3cCGEaCgy4EII0VBkwIUQoqHIgAshREORARdCiIYyZy8Ug7fSS457qceJZ1zlHRMhnDhY7EIcWfJlHjJ/fJmHqz+4vJ60vad7mfbdrNPJrXX5uOvVEm2/Pk5V6Pta6RwArrx/S+tt2vfyOM0rcKrgXgzbdeoVBAAPtK8nbcy7BQDW8nTdpwvueTFw7pHznb20/cXB/bTv0Tz1cBkEnj5Leeqx8mDrGu1bBV42zGuFrRnge3+x5BUPvr37Bm3frFNvirWce+/kpJJFdE+/3rovaetm/BlZieLHCVsV91g5Sbxslgp+vSMtfj3m4VRk3PtqtZV6nAxLfi62t4ltWg6SRQQFJJgzU0kKoQCAB54ss6A3cCGEaCgy4EII0VBkwIUQoqHIgAshREOZs4gJgIiYXqRf5o9IyDwAjI6mAk2gl8DbaQj6PWs8VPnMyjpt/97V80nb+3uv0r4lUTAGPtsWvzY+nrRF4t+JLBXuukEo/X/sridteSDQrWR8Q/9hmAqka0E4fpcIrN1AtFkxLipt1OnY0V4sZdOH0pek1PjDLR7mH1F5ur71IB/4WpaKm9/WuUj7PhSEsXdIruro/nXIfl6tuLh9ulhP2srgzLaMh8cz0fTtmoeEszD/4wUXY6O0EBtVKuieaXMR+hut00nbdskF3VdH6bkog2chKlc/3iGV7YPX5KpD/mFGYVNv4EII0VBkwIUQoqHIgAshREORARdCiIayrwE3s6fM7LKZvbCn7biZPWtmL03+5GFlQggh7hjTuEh8EsDvAfjDPW1PAPi8u3/czJ6
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(A)\n",
    "A.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The image, stored in a matrix $A$ is obtained from a certain original $A_0$ via convolution with a filter $C$ addition of a noise. The filter $C$ 'blurs' the image, simultaneously changing the image size from $16\\times 51$ to $25\\times 60$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mat2vec(A):\n",
    "    return np.reshape(np.flipud(A), np.prod(A.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def vec2mat(a, shape):\n",
    "    return np.flipud(np.reshape(a, shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Representing the images as vectors, we can write down the filtering as follows:\n",
    "$$\n",
    "a_0\\to a = C a_0 + \\epsilon,\n",
    "$$\n",
    "where $\\epsilon$ is a noise vector (consisting from iid normal variables)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Lets inspect how the filter $C$ acts on the images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACJCAYAAAA8GMMgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAIW0lEQVR4nO3dT4xdZRnH8e/PtrRSSQShBCiKmsZIiNRkUklwwR+BqsTiwgQSky6M40ISTDSmsvFPQuLGPxtjUrWhiQIhaqUxRKhVgwujFK1SBKQhRZppOhA04gYFHhf3TJxAp9O5986c8t7vJ5mce945d94nz0x/PXnnnDOpKiRJb3xv6rsASdJ4GOiS1AgDXZIaYaBLUiMMdElqhIEuSY0YKdCTbE3yZJLDSXaMqyhJ0tJl2OvQk6wC/gZcBxwFHgZuqaq/LvSeM7K21rF+qPkkaVK9yD+er6rzFjtu9QhzbAEOV9XTAEnuAbYBCwb6OtbzgVw7wpSSNHl+WT9+5lSOG2XJ5SLg2Xn7R7sxSVIPRjlDzwnGXrd+k2QamAZYx5kjTCdJOplRztCPAhfP298IzLz2oKraWVVTVTW1hrUjTCdJOplRAv1hYFOSdyY5A7gZ2DuesiRJSzX0kktVvZzkVuABYBWwq6oeG1tlkqQlGWUNnaq6H7h/TLVIkkbgnaKS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGrB7lzUmOAC8CrwAvV9XUOIpq3QMzB/suoXc3XLi57xJOC/4snJw/J0szUqB3rq6q58fwdSRJI3DJRZIaMWqgF/BgkkeSTI+jIEnScEZdcrmyqmaSbAD2JXmiqh6af0AX9NMA6zhzxOkkSQsZ6Qy9qma67SywB9hygmN2VtVUVU2tYe0o00mSTmLoQE+yPslZc6+B64FD4ypMkrQ0oyy5nA/sSTL3de6qql+MparGnQ6XYi10udzpUJsGJul74eWb4zF0oFfV08DlY6xFkjQCL1uUpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1YtFAT7IryWySQ/PGzkmyL8lT3fbs5S1TkrSYUzlDvxPY+pqxHcD+qtoE7O/2JUk9WjTQq+oh4IXXDG8DdnevdwM3jbcsSdJSDbuGfn5VHQPothvGV5IkaRirl3uCJNPANMA6zlzu6SRpYg17hn48yQUA3XZ2oQOramdVTVXV1BrWDjmdJGkxwwb6XmB793o7cN94ypEkDetULlu8G/gd8J4kR5N8Cvg6cF2Sp4Drun1JUo8WXUOvqlsW+NS1Y65FkjQC7xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1YvViByTZBdwIzFbVZd3YV4BPA891h91eVfcvV5Eavxsu3Nx3CVrEAzMH+y5BbzCncoZ+J7D1BOPfqqrN3YdhLkk9WzTQq+oh4IUVqEWSNIJR1tBvTfKXJLuSnD22iiRJQxk20L8LvBvYDBwDvrHQgUmmkxxIcuC/vDTkdJKkxQwV6FV1vKpeqapXge8BW05y7M6qmqqqqTWsHbZOSdIihgr0JBfM2/04cGg85UiShpWqOvkByd3AVcC5wHHgy93+ZqCAI8BnqurYopMlzwHPdLvnAs8PVXVb7MOAfbAHc+zD63vwjqo6b7E3LRroyyXJgaqa6mXy04h9GLAP9mCOfRi+B94pKkmNMNAlqRF9BvrOHuc+ndiHAftgD+bYhyF70NsauiRpvFxykaRG9BLoSbYmeTLJ4SQ7+qihD91jEmaTHJo3dk6SfUme6rZNP0YhycVJfp3k8SSPJbmtG5+0PqxL8ockf+768NVufKL6AJBkVZI/Jfl5tz+JPTiS5NEkB5Mc6MaW3IcVD/Qkq4DvAB8GLgVuSXLpStfRkzt5/ZMrdwD7q2oTsL/bb9nLwOer6r3AFcBnu+//pPXhJeCaqrqcwT0dW5NcweT1AeA24PF5+5PYA4Cru6fXzl2uuOQ+9HGGvgU4XFVPV9V/gHuAbT3UseIWeHLlNmB393o3cNNK1rTSqupYVf2xe/0ig3/IFzF5faiq+ne3u6b7KCasD0k2Ah8Fvj9veKJ6cBJL7kMfgX4R8Oy8/aPd2KQ6f+4u2267oed6VkySS4D3A79nAvvQLTUcBGaBfVU1iX34NvBF4NV5Y5PWAxj8Z/5gkkeSTHdjS+7Don+xaBnkBGNeajNhkrwF+Anwuar6V3KiH4u2VdUrwOYkbwX2JLms55JWVJK5v4T2SJKrei6nb1dW1UySDcC+JE8M80X6OEM/Clw8b38jMNNDHaeL43MPO+u2sz3Xs+ySrGEQ5j+qqp92wxPXhzlV9U/gNwx+vzJJfbgS+FiSIwyWXq9J8kMmqwcAVNVMt50F9jBYml5yH/oI9IeBTUnemeQM4GZgbw91nC72Atu719uB+3qsZdllcCr+A+DxqvrmvE9NWh/O687MSfJm4EPAE0xQH6rqS1W1saouYZADv6qqTzJBPQBIsj7JWXOvgesZPMF2yX3o5caiJB9hsHa2CthVVXeseBE9WODJlT8D7gXeDvwd+ERVNfsn/5J8EPgt8Cj/Xze9ncE6+iT14X0MftG1isGJ1b1V9bUkb2OC+jCnW3L5QlXdOGk9SPIuBmflMFgGv6uq7himD94pKkmN8E5RSWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiP+B5iOd9fD+omsAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.zeros((16,51))\n",
    "X[5:15,15] = 1; X[5:15,30] = 1; X[5:15,40] = 1\n",
    "X[5,10:20] = 1; X[5,30:41] = 1; X[15,30:41] = 1\n",
    "plt.imshow(X);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACsCAYAAABikvffAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAN1klEQVR4nO3dbYwdV33H8d8vy+66fiDYSRysJCVJlRcgBI60SkHpC0MoSinC4UWkRGplJKTlBUhBCiIubwKVkPKipbwAVXIbK0ZAUFSSxqpQwQRQWgmFOGnaODiQKOTB2Hjz4NgbIvlh8++LO4svu2d87+zMnesz+/1I1r337PjOOR77l5M558xxRAgAkJ8Lxl0BAMDKEOAAkCkCHAAyRYADQKYIcADIFAEOAJmqFeC2b7T9K9vP2t7ZVKUAAIN5pfPAbU9I+rWkv5R0SNKjkm6NiF+W/Z4pT8carVvR+QBgtZrXsVci4pKl5W+r8Z3XSXo2Ip6TJNvfk7RdUmmAr9E6/blvqHFKAFh9fhz/9kKqvM4tlMskvdT3+VBRBgBoQZ0euBNly+7H2J6VNCtJa7S2xukAAP3q9MAPSbqi7/Plkg4vPSgidkXETETMTGq6xukAAP3qBPijkq6xfZXtKUm3SNrbTLUAAIOs+BZKRJyx/TlJP5Q0IWl3RDzVWM0AAOdU5x64IuIHkn7QUF0AABWwEhMAMkWAA0CmCHAAyBQBDgCZIsABIFMEOABkigAHgEwR4ACQKQIcADJFgANApghwAMgUAQ4AmSLAASBTBDgAZIoAB4BMEeAAkCkCHAAyRYADQKYIcADIVK09MW0/L2le0oKkMxEx00SlAACD1Qrwwoci4pUGvid7ExdtWl648cLksW9tXJ8sP3Xh1PKyd6Qv08kNy/8H6tTbnTz2dOJ0pzdE8tiFDW8lyyfml59vcj59vsk3lpdNnUifb3o+fb6p188sLzt+KnnsBccSJzx2PHnswquvLStLXjup0vVLXTspff1S104a7fVLqXtNpfR1Hdk1lZLXNXVNVwNuoQBApuoGeEj6ke3HbM82USEAwHDq3kK5PiIO294saZ/tpyPi4f4DimCflaQ1WlvzdACARbV64BFxuHidk/SApOsSx+yKiJmImJnUdJ3TAQD6rDjAba+zvWHxvaSPSjrQVMUAAOdW5xbKpZIesL34Pd+NiP9spFYYSmrGQmq2gpSesVA6W2HD6WTxgiYTpWV9gNRMhvTshrLvmHq95PBROI9nC0nNXL+U+te0rLzsO4aPnCml/zCS37xKZ6GsOMAj4jlJ72+wLgCACphGCACZIsABIFMEOABkqoml9BiT2surSwa71q4/mSx/M/W9yUEwKd03qDII1q4qg5VS/eXxVQYrpWauX0r9ayrVH7CuFkNlg5urET1wAMgUAQ4AmSLAASBTBDgAZIoAB4BMMQslY3WXV5fNVrh4/e+T5aldO1KzGKSmlmi3p8psE6n9zTSauH4p9a+pVH/GUf1l92U16zp64ACQKQIcADJFgANApghwAMgUg5gZq7u8umyw653rTgxdh9QgmNTUEu32VBmslNp/Fvuorl9KlWsqtf+c+FRsMYgJAMgKAQ4AmSLAASBTBDgAZGpggNvebXvO9oG+sk2299l+pnjdONpqAgCWGmYWyj2SviHpW31lOyU9FBF32d5ZfL6j+erhnGoury6brfCutfV3+G5miXZ7qsw2kdrfTKPt65cyuhlHVTf6WP4d60qO7LqBPfCIeFjS0r8R2yXtKd7vkXRTs9UCAAyy0nvgl0bEEUkqXjc3VyUAwDBGvpDH9qykWUlao7WjPh0ArBor7YEftb1FkorXubIDI2JXRMxExMykpld4OgDAUisN8L2SdhTvd0h6sJnqAACGNfAWiu17JW2TdLHtQ5LulHSXpPtsf1rSi5JuHmUlkVb3+RhlsxWumn65XsVKVH3GRpsbPVSZbSK1v5nG+XD9yrS/0cf4NwA5XwwM8Ii4teRHNzRcFwBABazEBIBMEeAAkCkCHAAyxYYOGau7vLpssOua6d/Vq1hFZYOb0prW6lBlsFJqfzON8/n6pYx2ow8GMRfRAweATBHgAJApAhwAMkWAA0CmCHAAyBSzUDJWd3l12WyFq992vF7FGtPiLJQqs02k1jfTyPP6LdfMRh/0OxfxJwEAmSLAASBTBDgAZIoAB4BMMYiZsbrLq8sGu66aLFlXrnYHx+Z0UWvnqjJYKbX/LPZmrl/K+Ac8qz4nvnxwc/WhBw4AmSLAASBTBDgAZIoAB4BMDQxw27ttz9k+0Ff2Zdu/tf1E8etjo60mAGCpYWah3CPpG5K+taT8nyLiHxqvUcYWXk3MLEiVnUNqfL1szP3+mQ9W+u7cXP3jny8rSy94lxZGWxXgvDSwBx4RD0uqlkIAgJGrcw/8c7b/r7jFsrGxGgEAhrLSAP9nSX8maaukI5L+sexA27O299vef1rpp7UBAKpbUYBHxNGIWIiItyT9i6TrznHsroiYiYiZSU2vtJ4AgCVWFOC2t/R9/KSkA2XHAgBGY+AsFNv3Stom6WLbhyTdKWmb7a3qTQp4XtJnRldFlJmc97ir0BkT8+m+TNlzN1LP6Sh7pkf7hn++yXNnLlxW9szJdyaP/c3JS5LlL7y5aVnZ737/9uSxr7yxblnZm2+U/J/5fPrPvuxarUYDAzwibk0U3z2CugAAKuA/ZQCQKQIcADJFgANAptjQIWOTb4y7Bt1RPiA8/OBm2QYE58/g5nKpAcsqg5VSesAyNVgplQxYVhysZPD+LHrgAJApAhwAMkWAA0CmCHAAyBQBDgCZYhZKxqZOlG1vgKrKZ/QMPzslx2X3qRknVWabSPWXx1edbcLsq7PogQNApghwAMgUAQ4AmSLAASBTDGJmbHr+rXFXoTPKB4TLBjFT5fktu6/7LG+p/vL4qoOVDN6fRQ8cADJFgANApghwAMgUAQ4AmRoY4LavsP1T2wdtP2X7tqJ8k+19tp8pXjeOvroAgEXDzEI5I+n2iHjc9gZJj9neJ+lTkh6KiLts75S0U9Ido6sqlpp6/cy4q9AZ5TN6yvo4qZkTo1l2L41udkrtzRik2svjq842YfbVWQN74BFxJCIeL97PSzoo6TJJ2yXtKQ7bI+mmEdURAJBQ6R647SslXSvpEUmXRsQRqRfykjY3XjsAQKmhA9z2eknfl/T5iDhR4ffN2t5ve/9pnVxJHQEACUMFuO1J9cL7OxFxf1F81PaW4udbJM2lfm9E7IqImYiYmVTJPTQAQGXDzEKxpLslHYyIr/X9aK+kHcX7HZIebL56AIAyw8xCuV7S30p60vYTRdmXJN0l6T7bn5b0oqSbR1JDlJo6fmrcVRipNp94UT6jp+yfSKrvM5rnpkij2xSi7mYMUv3nm1SdbcLsq7MGBnhE/LfK/2be0Gx1AADDYiUmAGSKAAeATBHgAJApNnTI2AXHur0990KL56o+IJz6pzOaZffS6DaFqLsZg1R/eXzVwcrUtVqtWzzQAweATBHgAJApAhwAMkWAA0CmCHAAyBSzUHJ27Pi4a9AZZTN6prS+wreMatl9+juqbgqRVHMzBqn+8vgqs02k9LVqc8bS+YQeOABkigAHgEwR4ACQKQIcADLFIGbGFl59bdxV6I6SAeGyHk79wc0qy+7Lyqs9Uzyl7rO8pfrL46sMVkpi8L4PPXAAyBQBDgCZIsABIFMEOABkigAHgEw5or1Hodt+WdILxceL1czG2uerLrevy22TaF/uuti+d0XEJUsLWw3wPzqxvT8iZsZy8hZ0uX1dbptE+3LX9fb14xYKAGSKAAeATI0zwHeN8dxt6HL7utw2ifblruvt+4Ox3QMHANTDLRQAyFTrAW77Rtu/sv2s7Z1tn79ptnfbnrN9oK9sk+19tp8pXjeOs4512L7C9k9tH7T9lO3bivLs22h7je1f2P7fom1fKcqzb1s/2xO2/8f2fxSfO9M+28/bftL2E7b3F2Wdad8grQa47QlJ35T0V5LeI+lW2+9psw4jcI+kG5eU7ZT0UERcI+mh4nOuzki6PSLeLekDkj5bXLMutPG
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = mat2vec(X)\n",
    "plt.imshow(vec2mat(C @ x, (25, 60)));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACJCAYAAAA8GMMgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYyElEQVR4nO2deXRc9XXHv3dGI2m02LK8IS94wTYGG28oZjEBO5glgQJpSQollLZJSFtSSEoOJWl70gbS5vScZulpQuKkHGgoS0jCEuKwhJ1ggu1gB+MN413W5kW7NJJmfv1D49aF+322JSOlT9/POT7S3J/uvN+77707z+/3nXsthAAhhBD//0kM9QSEEEKcGJTQhRAiJiihCyFETFBCF0KImKCELoQQMUEJXQghYsKAErqZXWpmW8xsm5ndfqImJYQQ4vix/urQzSwJYCuAiwDsBbAawLUhhI3MJ1lWGgoqK/3BHNlOERkAELrJ51EB3yfrMf+9UtwnkeRzyPX6c0h0+tsBgFyabCsRcSx6yL6m+Nysi39eh0KyrcDnXVHa7tpbG0upT66IbCZJXTjG45OMiEM242/MsnxTgYQuXZqhPplmf2dzKb6dqDmg2N8net6Dz6+rlRwIAEZClyvq53dUcuQcijh+9Lwr4McV5NpDMmo7fCjRRXJDQcTbFZL5Zfl11J/b6O6dNftDCGOP9ncRUz0qiwFsCyFsBwAzexDAlQBoQi+orETVbbe4Y8lOfy9tip9EACBXU+Las5U91CdVW+jae6q6qU9ZRScda2vwk9mIjfwqbjndn1+itJf6WJ1/QYbxPMEUbU7Tsc7p/v5aJ8+0Hz1ntWt/4XtnUZ+2Kb69pyLiQmWfNREf7hXjWulY8+6Rrj3VxK+s3lJ/EvOr36E+bz8x07V3TODzTjVHXN2n+fvUu6uMupy+eIdr3/LyNOpT0OEnn46p/DpCRL4yci3TmwgAliE3RmP4+Z076F/LKOfXUYhItOVv+e/XNYbPO3tylz+3Vn79WzrqU9xn1w1f3HUsfzeQRy4TAew54vXevE0IIcQQMJCE7n3UveejzMxuNLM1ZrYm29Y2gM0JIYSIYiAJfS+AyUe8ngRg37v/KISwIoRQHUKoTpbx/yoKIYQYGANJ6KsBzDSzaWZWCOAaAI+fmGkJIYQ4Xvq9KBpC6DWzzwJ4CkASwN0hhLcinRIBIe0vEJVNaXbtHW+O4u9X7C9WpNJ8MScUkIUUpiIBMOKBEXSs/XJ/cXHRtZupz8uvzHHt2QiFSZIMFW/kC58XXf06HXv+vsWu/exr36A+m684ybWPuOc9/zH7Hw7u8hfmSyv5QnP7AX+xu7CBn67ZzUQ9BWDk+Qdde2umgvosOcdf2999x6nUp/dP/UXM5Db+P9PJ5+2hYzUvTHbt2Yl80e+dJ6e79jCPiws6uvy4lo/qoD7Z1/l12TnbXygs2smVNgVzW1y7veovaANA2ww/DsUR28lU8gXJ3Hl+DkKGn3fp9b4oon0mF1mkt/jz66w6/sXSdzMQlQtCCCsBrBzwLIQQQgwYfVNUCCFighK6EELEBCV0IYSICUroQggRE/pdy6U/FE2fGCbccZM7lt7gqzXap3PFyqwZta694ZGTqU/TAn/1OdkcsT5c5a/aA0C2zf+Kb6KNf4U+V+6vZhfV8jn0lPnHaeYCrpLY9wT53j2ArsX+l7zGPcxVM+Nu2e7a39gylfqkR/lqltQrXDnUMt//unfyAP869YxFPA67XvTj0DOLK20m/9A/Fp03H6I+9TVE+RFVwySiLkuC1OLJjeAql7FVvlLj4KbRfDusvlHE7V6WKMwAIJnx3y85lX+xMLvdVwJVLayjPs0rq/ztXHiA+tgTXA3VU+7Pu3Mc39feCv9YGKszA2DkJP8YtW3myqHtt926NoRQTf8gj+7QhRAiJiihCyFETFBCF0KImKCELoQQMUEJXQghYoISuhBCxIQB1XI5brJGO3m0TyPyH9LJBAB2vuYXL+o+k3c5KSjyJYPWyCVxvRESpIKDfgh7x/PiPKwQWJRMrGKmX2Cq5kkuTWybFdEB6YAvTyy7aS/12fGQ35EnvcwvrAQABa/68sTEMn9/AKBonS8tG1HdSH12P8vj0Fvuy85sTzH1+cy3HnTtd9xzLfVJzvVlkOn1XAradio/T5j01Ur4ce161i+GZtW80FYPkd5GtTccNZpLENs2+Mcv/UI59em+0Jfy1WwcT32Kz/d9in7CpYnl19fQsR2bfBlkahyXt4K0HSyq57LlrgZ/fkXzSXGw40B36EIIEROU0IUQIiYooQshRExQQhdCiJighC6EEDFhUFUuiW5DyR5S9OhUX5lipXylPVnvqwdKR/JiWolf+S2tMhW8AE+6NEI1U+8rJSpf4J+Vuy/xV8Arz2ygPnW1fuGexISI+LTxOWTLfL8D9/HCZnM+5bdl2/Dg6dSnc7wf15GP8EJEXZf5rdwa91ZQnwpeMwuZsUTlMpYf1xV/9QeuveOaCPVSxj+ubTN5gblpD/G3y33BLz63q4YX2jLSxayXqMsAoKDcn9/olVwF1DyDz2HB8i2ufX3NROpT/JJ/XSbJ+QMAnXV+Qa/ORdyn8W2/jSIApKv8Nn3pZ7g6p3Wab08t4idk5+YK195V77ezOx50hy6EEDFBCV0IIWKCEroQQsQEJXQhhIgJSuhCCBETlNCFECImDKinqJntBNAKIAug92g970pnVYU5//Yn7ljj5jGuvaDD7/MHAD2TfAnZ+Ke5RKv+Et+nMM2lZSN+7sujAKDkOl9aNqKISyfffm66a++q4kWXSnb5cs/MaC5bzEX0fSxq9CV2iXm8QFDpz/xCWx1X8OJc7XVEipXic0u0+3NLNfP7j+4pXIIYuvz3K67jqt2uyf55UrapkPpkRvv7lJ3Mz4XRo3iRq0Mb/GsiEaGcnLZkt2vv/DqXDDac6cfhlGU7qM9bOybQsYo1fsGqlrN5katchz+HUev4MTq0yL9epkzhRdyaH+XzbppHCgRGFEMreseXdkYV2iuc70sae1dzKe+WO//6mHqKnggd+rIQwv4T8D5CCCEGgB65CCFETBhoQg8AnjaztWZ244mYkBBCiP4x0EcuS0II+8xsHIBnzGxzCOGlI/8gn+hvBIDCcf4zWCGEEANnQHfoIYR9+Z8NAB4BsNj5mxUhhOoQQnXByJKBbE4IIUQE/U7oZlZqZuWHfwdwMYANJ2piQgghjo+BPHIZD+ARMzv8PveHEJ6Mcsh2FODgOr/nYXK6X+kstzOiAlmrP/2mGfxzKnQS+V9E39D9y7hOrOgVXwa16+QIbdlEIpGM6OGYzPjzjpImpg5FxGGOX9Gws4ZLNDMXk56Zr/qV8gBgzFJfANX1oi/JA4DMIv9cSO3j50Kmm+9rst0fKzmTi7PSP/erCYYP816o9pLfK7J8IZeC1m71rwcACGP886R8I5dO7n7O762a+DM+h9xb/qPQjRt45c10RM/Mpjm+zG/hFN6vdsuTfr/apjmkfCSA9E5fnrzvgN8bFACy53LpZNF2X4LYW8rPrdIaIlX9PV5tsbnJf1KRHNl/Cflh+p3QQwjbAcwf8AyEEEKcECRbFEKImKCELoQQMUEJXQghYoISuhBCxIRB7SkakgE9FWTVutkv6FNxGl8tzj3vqwpaZ/KV8REn+eqO9I+5UqP+g/z9WEGt5CFeIKzydF9d0baKKx5mXPW2a9/9wxnUp2kWHcLc8X7/0j2/9AuHAcChuf4xapvDC2O1bycFhyJ80r/11Syd87lCIVnnzw0AsqQvbeYVrrTJkZaZ3W1+H1sAyM72VSltdbzokvXy4nMFjf451DYvomdugx8H28q/1Jcd6cenZA9XsoTFXDWT6PBVOOvXnkJ9isimCg/ye87Emf4ckuv5tdzTzlNeT7mvMinbxefQtMw/J3O7+Bws6x/z3hE8zxwrukMXQoiYoIQuhBAxQQldCCFighK6EELEBCV0IYSICUroQggREwZVtgiAfoQkW/yp9OzwpYkAkF3iSxATu3gRp+wqX0K2fwEvjGURhZ9Kany9VfvJvA/hwY2+XC47mfc
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = mat2vec(A)\n",
    "a0 = np.linalg.pinv(C) @ a\n",
    "A0 = vec2mat(a0, (16, 51))\n",
    "plt.imshow(A0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "u, s, vh = np.linalg.svd(C, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1500, 816), (816,))"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.shape, s.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACJCAYAAAA8GMMgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYyElEQVR4nO2deXRc9XXHv3dGI2m02LK8IS94wTYGG28oZjEBO5glgQJpSQollLZJSFtSSEoOJWl70gbS5vScZulpQuKkHGgoS0jCEuKwhJ1ggu1gB+MN413W5kW7NJJmfv1D49aF+322JSOlT9/POT7S3J/uvN+77707z+/3nXsthAAhhBD//0kM9QSEEEKcGJTQhRAiJiihCyFETFBCF0KImKCELoQQMUEJXQghYsKAErqZXWpmW8xsm5ndfqImJYQQ4vix/urQzSwJYCuAiwDsBbAawLUhhI3MJ1lWGgoqK/3BHNlOERkAELrJ51EB3yfrMf+9UtwnkeRzyPX6c0h0+tsBgFyabCsRcSx6yL6m+Nysi39eh0KyrcDnXVHa7tpbG0upT66IbCZJXTjG45OMiEM242/MsnxTgYQuXZqhPplmf2dzKb6dqDmg2N8net6Dz6+rlRwIAEZClyvq53dUcuQcijh+9Lwr4McV5NpDMmo7fCjRRXJDQcTbFZL5Zfl11J/b6O6dNftDCGOP9ncRUz0qiwFsCyFsBwAzexDAlQBoQi+orETVbbe4Y8lOfy9tip9EACBXU+Las5U91CdVW+jae6q6qU9ZRScda2vwk9mIjfwqbjndn1+itJf6WJ1/QYbxPMEUbU7Tsc7p/v5aJ8+0Hz1ntWt/4XtnUZ+2Kb69pyLiQmWfNREf7hXjWulY8+6Rrj3VxK+s3lJ/EvOr36E+bz8x07V3TODzTjVHXN2n+fvUu6uMupy+eIdr3/LyNOpT0OEnn46p/DpCRL4yci3TmwgAliE3RmP4+Z076F/LKOfXUYhItOVv+e/XNYbPO3tylz+3Vn79WzrqU9xn1w1f3HUsfzeQRy4TAew54vXevE0IIcQQMJCE7n3UveejzMxuNLM1ZrYm29Y2gM0JIYSIYiAJfS+AyUe8ngRg37v/KISwIoRQHUKoTpbx/yoKIYQYGANJ6KsBzDSzaWZWCOAaAI+fmGkJIYQ4Xvq9KBpC6DWzzwJ4CkASwN0hhLcinRIBIe0vEJVNaXbtHW+O4u9X7C9WpNJ8MScUkIUUpiIBMOKBEXSs/XJ/cXHRtZupz8uvzHHt2QiFSZIMFW/kC58XXf06HXv+vsWu/exr36A+m684ybWPuOc9/zH7Hw7u8hfmSyv5QnP7AX+xu7CBn67ZzUQ9BWDk+Qdde2umgvosOcdf2999x6nUp/dP/UXM5Db+P9PJ5+2hYzUvTHbt2Yl80e+dJ6e79jCPiws6uvy4lo/qoD7Z1/l12TnbXygs2smVNgVzW1y7veovaANA2ww/DsUR28lU8gXJ3Hl+DkKGn3fp9b4oon0mF1mkt/jz66w6/sXSdzMQlQtCCCsBrBzwLIQQQgwYfVNUCCFighK6EELEBCV0IYSICUroQggRE/pdy6U/FE2fGCbccZM7lt7gqzXap3PFyqwZta694ZGTqU/TAn/1OdkcsT5c5a/aA0C2zf+Kb6KNf4U+V+6vZhfV8jn0lPnHaeYCrpLY9wT53j2ArsX+l7zGPcxVM+Nu2e7a39gylfqkR/lqltQrXDnUMt//unfyAP869YxFPA67XvTj0DOLK20m/9A/Fp03H6I+9TVE+RFVwySiLkuC1OLJjeAql7FVvlLj4KbRfDusvlHE7V6WKMwAIJnx3y85lX+xMLvdVwJVLayjPs0rq/ztXHiA+tgTXA3VU+7Pu3Mc39feCv9YGKszA2DkJP8YtW3myqHtt926NoRQTf8gj+7QhRAiJiihCyFETFBCF0KImKCELoQQMUEJXQghYoISuhBCxIQB1XI5brJGO3m0TyPyH9LJBAB2vuYXL+o+k3c5KSjyJYPWyCVxvRESpIKDfgh7x/PiPKwQWJRMrGKmX2Cq5kkuTWybFdEB6YAvTyy7aS/12fGQ35EnvcwvrAQABa/68sTEMn9/AKBonS8tG1HdSH12P8vj0Fvuy85sTzH1+cy3HnTtd9xzLfVJzvVlkOn1XAradio/T5j01Ur4ce161i+GZtW80FYPkd5GtTccNZpLENs2+Mcv/UI59em+0Jfy1WwcT32Kz/d9in7CpYnl19fQsR2bfBlkahyXt4K0HSyq57LlrgZ/fkXzSXGw40B36EIIEROU0IUQIiYooQshRExQQhdCiJighC6EEDFhUFUuiW5DyR5S9OhUX5lipXylPVnvqwdKR/JiWolf+S2tMhW8AE+6NEI1U+8rJSpf4J+Vuy/xV8Arz2ygPnW1fuGexISI+LTxOWTLfL8D9/HCZnM+5bdl2/Dg6dSnc7wf15GP8EJEXZf5rdwa91ZQnwpeMwuZsUTlMpYf1xV/9QeuveOaCPVSxj+ubTN5gblpD/G3y33BLz63q4YX2jLSxayXqMsAoKDcn9/olVwF1DyDz2HB8i2ufX3NROpT/JJ/XSbJ+QMAnXV+Qa/ORdyn8W2/jSIApKv8Nn3pZ7g6p3Wab08t4idk5+YK195V77ezOx50hy6EEDFBCV0IIWKCEroQQsQEJXQhhIgJSuhCCBETlNCFECImDKinqJntBNAKIAug92g970pnVYU5//Yn7ljj5jGuvaDD7/MHAD2TfAnZ+Ke5RKv+Et+nMM2lZSN+7sujAKDkOl9aNqKISyfffm66a++q4kWXSnb5cs/MaC5bzEX0fSxq9CV2iXm8QFDpz/xCWx1X8OJc7XVEipXic0u0+3NLNfP7j+4pXIIYuvz3K67jqt2uyf55UrapkPpkRvv7lJ3Mz4XRo3iRq0Mb/GsiEaGcnLZkt2vv/DqXDDac6cfhlGU7qM9bOybQsYo1fsGqlrN5katchz+HUev4MTq0yL9epkzhRdyaH+XzbppHCgRGFEMreseXdkYV2iuc70sae1dzKe+WO//6mHqKnggd+rIQwv4T8D5CCCEGgB65CCFETBhoQg8AnjaztWZ244mYkBBCiP4x0EcuS0II+8xsHIBnzGxzCOGlI/8gn+hvBIDCcf4zWCGEEANnQHfoIYR9+Z8NAB4BsNj5mxUhhOoQQnXByJKBbE4IIUQE/U7oZlZqZuWHfwdwMYANJ2piQgghjo+BPHIZD+ARMzv8PveHEJ6Mcsh2FODgOr/nYXK6X+kstzOiAlmrP/2mGfxzKnQS+V9E39D9y7hOrOgVXwa16+QIbdlEIpGM6OGYzPjzjpImpg5FxGGOX9Gws4ZLNDMXk56Zr/qV8gBgzFJfANX1oi/JA4DMIv9cSO3j50Kmm+9rst0fKzmTi7PSP/erCYYP816o9pLfK7J8IZeC1m71rwcACGP886R8I5dO7n7O762a+DM+h9xb/qPQjRt45c10RM/Mpjm+zG/hFN6vdsuTfr/apjmkfCSA9E5fnrzvgN8bFACy53LpZNF2X4LYW8rPrdIaIlX9PV5tsbnJf1KRHNl/Cflh+p3QQwjbAcwf8AyEEEKcECRbFEKImKCELoQQMUEJXQghYoISuhBCxIRB7SkakgE9FWTVutkv6FNxGl8tzj3vqwpaZ/KV8REn+eqO9I+5UqP+g/z9WEGt5CFeIKzydF9d0baKKx5mXPW2a9/9wxnUp2kWHcLc8X7/0j2/9AuHAcChuf4xapvDC2O1bycFhyJ80r/11Syd87lCIVnnzw0AsqQvbeYVrrTJkZaZ3W1+H1sAyM72VSltdbzokvXy4nMFjf451DYvomdugx8H28q/1Jcd6cenZA9XsoTFXDWT6PBVOOvXnkJ9isimCg/ye87Emf4ckuv5tdzTzlNeT7mvMinbxefQtMw/J3O7+Bws6x/z3hE8zxwrukMXQoiYoIQuhBAxQQldCCFighK6EELEBCV0IYSICUroQggREwZVtgiAfoQkW/yp9OzwpYkAkF3iSxATu3gRp+wqX0K2fwEvjGURhZ9Kany9VfvJvA/hwY2+XC47mfc
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# equivalent to least squares\n",
    "n = 816\n",
    "s0 = np.zeros(816)\n",
    "s0[:n] = s[:n]\n",
    "C0 = np.dot(u, np.dot(np.diag(s0), vh))\n",
    "A0_guess = vec2mat(np.linalg.pinv(C0) @ a, (16, 51))\n",
    "plt.imshow(A0_guess)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhRUlEQVR4nO3deXxV5YH/8c9zs+97QlaTELYAYRUEwb0KVtGp1sHWGdtBGbtN22lnqtNOW3+dTpexrW3HsTrW1rEuo47jgrS4gVSlSJA1QICEQEJCNrKRELLc5/dHLjYygIEs59x7v+/XK6/ce3Jzz5ebyzcnz3nOOcZai4iIBD6P0wFERGRsqPBFRIKECl9EJEio8EVEgoQKX0QkSIQ6HeBsUlNTbX5+vtMxRET8xubNm5ustWmn+5qrCz8/P5/S0lKnY4iI+A1jzMEzfU1DOiIiQUKFLyISJFT4IiJBQoUvIhIkVPgiIkHClYVvjLneGPNwW1ub01FERAKGKwvfWvuytXZlQkKC01FERAKGKwt/OKy1PL6hile21zkdRUTEVVx94NX5MMbwTGkNYSGGj5dkOh1HRMQ1Am4LH+CKyelsqW6lvbvX6SgiIq4RkIU/IzcBa2HvkQ6no4iIuEZAFv6kcfEA7FHhi4h8ICALPzM+kvAQD9UtXU5HERFxjYAsfI/HkJ0URU3LcaejiIi4hisLfyQOvMpR4YuIfIgrC38kDrzKSYrisIZ0REQ+4MrCHwk5SdE0Heuhq6fP6SgiIq4QsIWfnxIDwLZqnY9HRAQCuPAvmZhKdHgI975c5nQUERFXCNjCj4sM428uLmDPkQ7aunTErYhIwBY+wJz8JAD2HGl3OImIiPMCuvALUwfG8Q8e1WwdEZGALvysxCg8BiobO52OIiLiuIAu/LAQD5kJUTy8voK+fq/TcUREHBXQhQ9wXUkmXgtHO3ucjiIi4qiAL/xZeQM7bhs6TjicRETEWa4s/JG8iHlaXAQAjSp8EQlyriz8kbyIebqv8Ovbu4f9XCIi/syVhT+SxiVEEuoxHNLUTBEJcgFf+GEhHnKTo6lq1tRMEQluAV/4AAWpMZqLLyJBLygKPz8lhoPNXVhrnY4iIuKYoCj8grQYjvf2c7hVV8ASkeAVFIV/UUEyAG/uaXA4iYiIc4Ki8CdkxDEhPZZV2+ucjiIi4pigKHyAyyens7W6VefUEZGgFTSFPyUzjp4+L5VNmq0jIsEpaAp/Tt7AOP4buzWOLyLBKWgKPy8lmgvzk3h2c7WmZ4pIUAqawge4eU4OlY2dbKludTqKiMiYC6rCv2JyBgDbVPgiEoSCqvBTYsIJCzHUt+tUySISfIKq8D0eQ3pcJA06VbKIBKGgKnyAjPgI6jtU+CISfMas8I0xMcaYx4wx/2mM+fRYrfdUGfGRGtIRkaA0rMI3xjxqjGkwxuw8ZfkSY0y5MWa/MeZu3+JPAM9Za+8Elg1nvcMxUPjawheR4DPcLfzfAksGLzDGhAAPAEuBYuBWY0wxkANU+x7WP8z1nreM+Eg6uvvo6ulzKoKIiCOGVfjW2vXA0VMWzwP2W2srrbU9wNPADUANA6V/1vUaY1YaY0qNMaWNjY3DiXda+SnRAGysPDW2iEhgG40x/Gz+vCUPA0WfDTwP3GSMeRB4+UzfbK192Fo711o7Ny0tbcTDXTklg/S4CJ5879CIP7eIiJuFjsJzmtMss9baTuCzo7C+cxIe6uHSiWms2zvyfz2IiLjZaGzh1wC5g+7nALXn8gTGmOuNMQ+3tbWNaLCTUmIjaOvq1Tl1RCSojEbhbwImGGMKjDHhwHLgpXN5Amvty9balQkJCaMQDxKjw+jp99LV49i+YxGRMTfcaZlPARuAScaYGmPMCmttH/BFYA2wG3jGWls2/KgjJyk6DIDW470OJxERGTvDGsO31t56huWrgdXDee7RlBAVDkD5kXayE6McTiMiMjZceWqF0R7Dz08dmJr5izf2axxfRIKGKwt/tMfwJ4+LZ/mFuWytbuWFrYdHZR0iIm7jysIfC/csnUJSdBjffqGMdeW67KGIBL6gLfyE6DCevWshoSGGFY+V8pbm5YtIgHNl4Y/2GP5JRemxvPG1y5iYEcfnfreZLYdaRnV9IiJOcmXhj/YY/mDJMeE89jcXkhwTzj/9786P/gYRET/lysIfa+lxkdy5uJDdde3sOdLudBwRkVGhwve5riST8FAPn/3NJp0vX0QCkgrfJyU2gu/dMJW6tm7WlB3R/HwRCTiuLPyx2ml7qk/OySU5Jpxvv1jGx362ng0VzWO6fhGR0eTKwh/LnbaDeTyGp1dexLc+PoWO7l6+/uw2bemLSMBwZeE7aWJGHHcsLuSrV03kcOtxHnu3yulIIiIjQoV/BstmZjEuPpLH/3QQr1db+SLi/1T4ZxAdHsrdSydT0djJl57eQucJXfRcRPybKwvfqZ22p7phZhafnp/HK9vreGxDlaNZRESGy5WF79RO21MZY/j+X0wnPyWaP+5tcjSLiMhwubLw3eaSiWlsqGzmuy+V0dfvdTqOiMh5GdYVr4LFP19XjNdafvtuFfXt3fzklhlEh+ulExH/otYagrAQD/9y43TyU2L4l1d2s7uunR/eVMLsvCTCQ/VHkoj4B7XVObhjcSEP/dUcmjt7WP7wn7jyp+s4ptk7IuInVPjn6Jqp43j7G1fwzWunUH30OE9uPKijcUXEL7iy8N0yLfNMEqLCWLGogDkXJPGvq/dw04Pv0q+Ds0TE5VxZ+G6Zlnk2Ho/hiTvm83dXFPH+oVY+85v3WL2jjhN9/U5HExE5Le20HYbIsBC++rGJRIWH8rPX9/LHfU0snpDKd64vpig9zul4IiIfYtw8/jx37lxbWlrqdIwh6e7t57415Tzy9gEALi5K4dvXTWXSOBW/iIwdY8xma+3c035NhT+yqo928WxpNQ++VUFvv2VCeiyXT05nxaICMuIjnY4nIgFOhe+AqqZOnt5UzY7Drbyzv5mEqDAe/cxc5lyQ7HQ0EQlgKnyH7TnSzm2PvEfTsRPkJEXx8+UzVfwiMirOVviunKUTaCaPi+e1r17CPUsnYy3c/KsN/PMLO+no7nU6mogEERX+GEmKCedvLx3PM3ct4MrJ6Tz+p4Ms/vFaHli7X0frisiYcOWQjjHmeuD6oqKiO/ft2+d0nBFnraX0YAsPrqvgzT0NJEaHcefiQm5fmE9shGbKisj50xi+i22rbuX+1/eytryRqLAQvnzVBO66dLzTsUTET52t8LU56bAZuYn85rPz2Fbdyi/f3McPf7+HqqZOFhalMikjTvP4RWTEaAvfRfr6vXzjf3bw/JYaTv5YFk9I5dPz87i6eBwej3E2oIi4noZ0/Ex3bz8HmjpZU3aEx96toqWrl5KcBO5eMpmFRalOxxMRF1Ph+7ETff2s2lbHv60p50h7N99YMpm7Li3EGG3ti8j/pXn4fiwiNISb5uTw5tcvZXZeIj/6wx6u//e3OdLW7XQ0EfEzKnw/ER0eyhN3XMS9y6ay83A7N//qXd47cNTpWCLiR1T4fiQqPITbF+Zzx6ICOrr7uOWhDXziP97h3f1N9PZ7nY4nIi6nMXw/1dXTx2/eqeLBdRUcO9FHeIiHT83P46tXTSQhOszpeCLiEM3DD0DR4aF84fIibp6Tw7sVTTy3uYbfvlvFzsNt/PJTs8hMiHI6ooi4jIZ0/FxGfCR/MSuH362Yz1evmkjpwRYu+fFantlUrYuri8iHuHJIJ9DPpTOadtW2880XdrDlUCuLJ6Ry0+wcpuckMD4t1uloIjIGNA8/yHi9lgffquD+1/fS228J8Rhump3NN68t1vi+SIBT4Qepnj4vB5o6eWh9Bc+/f5gQj+Gh2+ZwVXGG09FEZJTowKsgFR7qYdK4OH7yyRk8vmIeEzPi+PyT77OxstnpaCLiABV+EDDGsHhCGj++qWRg+uYjG3ns3Sq6enThFZFgosIPItNzEnj7G5czPTuB77xUxpzvvc6zpdW061KLIkFBY/hByFrLO/ub+f7q3eyuayci1MM3lkzmpjk5JERpp66IP9NOWzmtvn4
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s0)\n",
    "plt.yscale('log');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACJCAYAAAA8GMMgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUyElEQVR4nO2dXYwk11XH/6eqq7vna2e/1/baiUNkiUQWMdLKCjIPToDIgIXDQ1AsIfkBsTwQKUggZPISQIrECx8vCMkQy5YgCZbAxEIWiWWCzBN4DYE42AHLcuzVzu7sx8zOV39Udx8epgcm9v2fmemenXZq/j9pNd339q26fU7Vqdq6/z7H3B1CCCF++MkmPQEhhBD7gwK6EEJUBAV0IYSoCAroQghRERTQhRCiIiigCyFERRgroJvZQ2b2PTN7w8we369JCSGE2Ds2qg7dzHIA/w3gZwBcBPAygEfd/b/YmHo+7VPFkWSfF3myfVDj1xxPD6HtAOBkc258zChYZNbB3sfYCGPYfjbHHczvD9yIYQN7U1+MMmaHcfsKMWm4+8gNbHsjjdlff8f23rvB993nbHP7bG/rpztZOwBYj5yYvR4ds9K/ds3dT/GZbFLb6QMB9wN4w93fBAAz+xqARwDQgD5VHMFPfPCxZF9523yyvX2qTifQnk9H53KOe7rfJO0FHRKfkcRvGfcN8na6vdbmBwEd0wnGdHlfVpIDcbD3E98zbqBBLd03KKIx6fb+CGMAwEnfSBeBwDxZn2yKtANAVgZ9vfTOcuK7aF9sWwD4d4qCaXDTNMiJz8mxAAADcv6FY4I5MNiNEcDtGp2XxVp6g8XNLh1Tu76W7li8Rsd8Y+nL36ed2xjnkctZAO9se39x2CaEEGICjBPQU5fO91zKzOy8mV0wswvdfmuM3QkhhIgYJ6BfBHDXtvd3Arj07g+5+xPufs7dz9XzqTF2J4QQImKcgP4ygHvM7ENmVgfwWQDP7c+0hBBC7JWRF0XdvWdmnwPwDQA5gCfd/bvhmCJHeUd68XPtjkayvXWKX3M6x9Lt5Rxf+ejPkL4iWC2JLntsWMkHZa10X97mC0C1VrovGpPxdRm6GDfSoihZBAP4YiVbBIv6Bnx9HIOCz5uqoUZYFI0VD6S9N6KPyLg8GGNkMT5amB1FAXNgPg/GeJ6ed+TXrB/5gp1jfHv11fS53Jjj5/9UM31A1vMg0Czxru2Mo3KBuz8P4PlxtiGEEGJ/0C9FhRCiIiigCyFERVBAF0KIiqCALoQQFWGsRdG9MigytMhP+TdOp68trTN8Bb48kV7SL+Y7dMz8THrJeqbOf4OdZ4FqZkDmXfLl+bVWWtHT2eBjylbaVdkGvybnZNUeAIyqXAKJAJF4MLUBECge6sGYBlMiBT95r3MZR1ZLj8sCvzIGgYRi0Ev7wgPFE7qBGqqT7mNqDICrXKJUFFQWEkh6Qp8TVVHkc2+QY6se+Ki2d/8hULkwZVq+zn1UEpVLb4rnJWDpDGwwy+dGE6r8ILpDF0KIiqCALoQQFUEBXQghKoICuhBCVAQFdCGEqAgK6EIIUREOVraYA+1j6WtI53h6THmS663mTqcrf9w5f5OOOTud7jteX6djmkFZmTbJHrRc8lTBV1rpMnxXWzN0zM2N9PbaLZ6xqmwHJV2YlC7K08QUX4GEjckJa4HMcLqZtvdMg2elmiq4j6Zq6b7aCLLFHpGpAkCrlz4WIgnreof7r9NOj+t1uV+d9UVyPWaGKHlZ4HOW6C5vcJ83G2kfNQM5cbNIx4Zo2pH/Nrppe2+skzJnANrT6TGDGveRkTnUgmNht+gOXQghKoICuhBCVAQFdCGEqAgK6EIIUREU0IUQoiIcqMrFc6A7l16DLufTK+PNY7z+093H0nWZPnLkMh3z4eZisv222jIdUw/qd3VJfbOrvbSSBQAuNtKSnrdJOwAsNNLbW5riappIQVGWadcPguRcWZZWNtRq3D5TRKUw1+AJ1I41NtLt9RYdc7RIjwGA2Ty9r0i9xGCqJgBY6aXVEMvlNB1zvcOVTcudtG9XA7+2iVKj1wuUMcznQXKuWpAYq0HUJ5FK6UgjfZ4fJ8cCAMzU0n4tjM+tw7LFAVjqpu19ucnP5au1dEKt0vl52emkfVGsjX9/rTt0IYSoCAroQghRERTQhRCiIiigCyFERVBAF0KIiqCALoQQFWEs2aKZvQVgFUAfQM/dz0Wf9wzoETXPYCYtfTs6y2VLH5gmssWpS3TMj9YXku2nci6JawTZfjpE2XW9xhOEzWR7l9E18rQUbJrIpgDgZp1Lp1jCqCg3V05kbCz5FQDMN9J2PdngydBO1VeT7aeLFT6mxvuOZGlJXJMVVgXQJyme1gfperAAsNxPyxOvBRLWhfo87bvWSPs2kjqudNLSyXaPn+qsLq4FssV6zqWqs/X08R1JEE820on2bqtzv87X0tsrAplxJ5CdXuul7T1d43JLxkLJZaLlWtpH5WyUVmx37IcO/RPufm0ftiOEEGIM9MhFCCEqwrgB3QF808xeMbPz+zEhIYQQozHuI5cH3P2SmZ0G8IKZve7uL23/wDDQnweA2pFjY+5OCCEEY6w7dHe/NPy7COBZAPcnPvOEu59z93P5DF/MEUIIMR4jB3QzmzGzua3XAD4F4NX9mpgQQoi9Mc4jlzMAnjWzre18xd3/IRxhwIAkijNSb/AIkUABXN52WyAZPEPkiWdynsGuYdxMHU/LCQvj824Xabll27mkaoMYru9BdsRAdjZVI7UQg+3Vs/R3nS34dz1eT0vLbq9zH91O7HOWtAPAqZzLIOcsPe9moBJjwrfVoCblzTwtaZzJuOwtktgx/2WBuLRGtrdecrllz9PfKTp+mIwWAI6SrJhMmggAZxtp395RLNMxJ/L09gribyA+x+ZIbOgT+wDAGrHr8jSXDG9Mp8/l3lRQA3iXjBzQ3f1NAB8bewZCCCH2BckWhRCiIiigCyFERVBAF0KIiqCALoQQFeFAa4oCoHUKWSKgPOP1AXNSOzADH1MQZUNhfIU5N37dK0DqA9IRQE7mFykeGmTlPqqfWM/49iI1C6NJlA2zQfKi+VpaOTCf80RNJ2pp9cKJjI85ThQ4ADCXMR9xnw+Ij6L6sgOk7TCX8cRv0xmvNzpN1DGRwoQlzeoO+LyZpCc695jiCQAapI99n6gvsh3ra4YqF56QrczT4ZAlAQO4wqte43NYL9KxLih3umt0hy6EEBVBAV0IISqCAroQQlQEBXQhhKgICuhCCFERFNCFEKIiHKxs0QHrpeVyg1762rJR8qRZq/10bb7lAZeC3RikE3rVLV13EgAagWyx42lp140Bl8Sx+bGalABws59O9tPqc4Fku8/d2x1BI8WSNZXBd2WJjfrBvQQfw6WWXGAHDJzIxALJZ0n82iXb2hyTnncZyCOjZFEd4qPQ5710Xyc4Fsp+en45scFOtPppyR5LMBf1RfaJ+hilB+eEp+0wCJJzsT6PZMFR4d4x0R26EEJUBAV0IYSoCAroQghRERTQhRCiIiigCyFERVBAF0KIinCgskVzIGflJ1tpydDSBq/N987GsWT7PKkNCAA50Qwt11bomKbxDG1MOnWjP0vHvNM9kWx/u3Ocjllozyfbl7vcPuuB5JNJ1eI6pOnDJcrcmBFpYJRZktmbZdcEgD7SclQAWKfb43NgtVpXndfmvNqfS7YzfwPAQvco7bvSOZLeT5sfW8vt9PHQ6fFTvddP39exDKgA0Ch4NsF+UHeVwTKQsvMVANqkLm4z4+drJEG82kvbe7Gb9ivAz792l0sqrUtquPJp7xrdoQshREVQQBdCiIqggC6EEBVBAV0IISqCAroQQlSEHVUuZvYkgIcBLLr7vcO24wD+GsDdAN4C8EvuvrTjtvpAQcQIxc206mJ1hqs43qyn1QODIInTUi+dAOtkka5jCcSKjJIk9LnRm6FjFjvpVfMrLb6afrOTTkS21uaqi7IX1MwcpG2UZVxVUNTSdlhp8Dmslem+9V4wpp/uu9HnNj2er9O+OaJ6ihQUzK8rA348XivT/rvcTasnAOBym/ddbaXVLMu
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 10\n",
    "s0 = np.zeros(816)\n",
    "s0[:n] = s[:n]\n",
    "C0 = np.dot(u, np.dot(np.diag(s0), vh))\n",
    "A0_guess = vec2mat(np.linalg.pinv(C0) @ a, (16, 51))\n",
    "plt.imshow(A0_guess)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3.10.8 ('nummethods')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "8d489487a004dbbc79692a52da1077223bea5d5bb7772308e4a4df2310821984"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}