diff --git a/Task 3/Lecture 3.ipynb b/Task 3/Lecture 3.ipynb new file mode 100644 index 0000000..aac17da --- /dev/null +++ b/Task 3/Lecture 3.ipynb @@ -0,0 +1,1244 @@ +{ + "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": [ + "" + ] + }, + "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+/atHnE7YcIERo0axfnnn09qamr16nrs2LGMGDGC5ORk+vTpQ9euXev8/zlkyBA2b97MhRdeCLgugs6fP5+tW7eSnp5OREQEUVFRPPnkk77/kkWaSxBODfQXY61ttg9LTk6269evr3Vs8+bNdOvWrdlqCDUvvPAC69ev5/HHHw9oHfr3JEHv+NbA1AdDsjWwLsaYbGtt8vHHtQIXkdAX5FMD/UUBHuTGjx/P+PHjA12GSHAKw9ZAdwRFgFtr6+28kMBrzm02kSYL09ZAdwQ8wKOjoykpKeG0005TiAchay0lJSVER4fGbAhxgDBvDXRHowFujHkOuALYZa3tWXXsVOBVIBHYAVxtrf3ZkwISEhLIz8+nuLjYk7dLM4iOjm70xiQRvwvhqYH+0pQV+AvA48A/ahybBrxvrZ1ljJlW9fNUTwqIioqqvptRRKROP++EZXfBtvfDvjXQHY0GuLX2Q2NM4nGHRwGXVP3zi8AaPAxwEZF6hcnUQH/xdA/8DGttEYC1tsgYc3p9LzTGTAQmAnTq1MnDjxMRxyn6ApZOdlxroDv8fhHTWvs08DS4buTx9+eJSIgrPwBrZsKnT7haA8e8AN3THNMa6A5PA/xHY0yHqtV3B2CXL4sSEYdSa6BbPB1mtRS4oeqfbwCW+KYcEXGk/SXw5q3wUhpEtHC1Bo78m8K7EU1pI3wZ1wXL9saYfGA6MAv4pzHmZuA7YEz9ZxARqYdaA73SlC6U6+r5pVQf1yIiTqLWQK8F/E5MEXGYysOw7klY/Zeq1sBM6HezWgM9oAAXkeZTqzVwGAyfrdZALyjARcT/1BroFwpwEfGvbatce91qDfQ5BbiI+Mf+Elh5n2tq4Gn/CePfhsSLAl1VWFGAi4hvqTWw2SjARcR3ft4By6aoNbCZKMBFxHtqDQwIBbiIeEetgQGjABcRz9TTGrh4YyGZWasoLC2jY2wM6UO7kJYUH+hqw5ICXETcV09r4OKcAjIW5VJWUQlAQWkZGYtyARTifuDpNEIRcaLqqYFXuqYGjn+71tTAzKy86vA+qqyiksysvEBUG/a0AheRxjWxNbCwtKzOt9d3XLyjABeRhrnRGtgxNoaCOsK6Y2yMn4t0Jm2hiEjdKg/DJ3+DeRfC9+tcrYE3r2ywrzt9aBdiomq3DsZERZI+tIu/q3UkrcBF5EQetgYevVCZmZWnLpRmoAAXkWN8MDUwLSlegd1MFOAi4qKpgSFHAS7idJoaGLIU4CJOpamBIU8BLuJEmhoYFhTgIk6iqYFhRQEu4hSaGhh2FOAi4U4PFA5bCnCRcKbWwLCmABcJR2oNdAQFuEg4UWugoyjARcKFWgMdRwEuEupOaA18GPr9Xq2BDqAAFwllag10NAW4SCgqPwBr/gKfzqtuDVx8qB+ZT35NYekXGuPqEApwkVCzbRW8dSeU7qxuDVy85QAZb+phwk6jJ/KIhIr9JbDoFtcDhSOjaj1QWA8TdiavVuDGmLuA3wMWyAVutNYe9EVhIlLFWtj0T8jKcLUGDkx3tQfWaA3Uw4SdyeMVuDEmHpgMJFtrewKRwLW+KkxEcLUGzr8K3pwIp/4KbvkIBt9/Ql93fQ8N1sOEw5u3WygtgBhjTAvgJKDQ+5JEhMrD8PFj8ETKsQcK35QFZ3Sv8+V6mLAzebyFYq0tMMbMBr4DyoCV1tqVPqtMxKkKN8Jbk10tgk1sDdTDhJ3J4wA3xrQDRgGdgVLgNWPMOGvt/ONeNxGYCNCpUyfPKxUJd3W0BrozNVAPE3Yeb7ZQLgW+tdYWW2srgEXAb45/kbX2aWttsrU2OS4uzouPEwlj21bBvBT45G+QNA4mfQY9rtTIV2mQN10o3wEpxpiTcG2hpALrfVKViFPsL4Gse2HTK5oaKG7zZg98nTHmdWADcBjIAZ72VWEiYc1a2PQqrMiAQ7/U2Roo0hiv+sCttdOB6T6qRcQZft7hesjCtlWQ0A9GPFZvd4lIQ3QrvUhzqTwMa+e5pgZGtIDLZ0PyTZoaKB5TgIs0h5qtgV0ud4X3KeoYEe8owEX8qXx/1QOF50Hr9jDmReg+St0l4hMKcBF/2fq+a6+7dCf0vQH+6yE9UFh8SgEu4mv7S1yDpza9WtUauBwS+we6KglDCnARX6nVGrgXBv4JBtyt1kDxGwW4iC/8vMP1kIXtqyHhfNcDhU/vFuiqJMwpwEW8UWdr4M0QoWeliP8pwEU8VbgRlt4BP2xSa6AEhAJcxF3l+10r7rXzoHUcXP0P6DZSrYHS7BTgIu6o2Rp43ni49CGIiQ10VeJQCnCRpti/u2pq4ImtgYtzCvQgBQkIBbhIQxppDVycU0DGotzqJ8IXlJaRsSgXQCEufqcAF6nPT9+6tksaaA3MzMqrDu+jyioqyczKU4CL3ynARY7nRmtgYWlZnaeo77iILynARWpyszWwY2wMBXWEdcfYGD8WKeKiuw1EwNUamHUf/H0Q7PvR1Rp47cJG+7rTh3YhJqr2PO+YqEjSh3bxZ7UigFbgIl61Bh7d51YXigSCAlycq4HWQHekJcUrsCUgFODiPJoaKGFCAS7O0oTWQJFQoQAXZ9DUQAlDCnAJf5oaKGFKAS7hS1MDJcwpwCU8aWqgOIACXMKLj1oDRUKBAlxCVq0xrqdE81iPrzlv88NqDRTHUIBLSKo5xvVM8yP/c+A5ztuQS0m7Ppx241NqDRRHUIBLSMrMyqO8opwJke8wpcXrHCaSByrGs/rACP6l8BaHUIBLSGq35yv+f8u/0zNiB+9WnscDFeP5gdMwew4FujSRZqMAl9BS1Rq4pNU8Smxbbi2/kxVH+gGu1kCNcRUnUYBL6KjRGvjdWWO4Zvswdh05dpFSY1zFaXQfsQS//bth0USYPxoiW8L45XS+8RnuHZ1CfGwMBoiPjWHm6F6aCiiOohW4BK9GpgZqjKs4nVcBboyJBZ4BegIWuMla+6kP6hKn09RAkUZ5uwJ/FFhhrf2tMaYlcJIPahInqzwMa5+A1TM1NVCkER4HuDGmLTAQGA9grS0Hyn1TljhSramBw+HyTE0NFGmANyvwXwHFwPPGmN5ANvBHa+1+n1QmzqGpgSIe8ebvpS2AvsCT1tokYD8w7fgXGWMmGmPWG2PWFxcXe/FxEpa2vg/zUuDTx6Hv9fCHz6D7KIW3SBN4E+D5QL61dl3Vz6/jCvRarLVPW2uTrbXJcXFxXnychJVarYGtXFMDRzyqka8ibvB4C8Va+4Mx5ntjTBdrbR6QCnzlu9IkLFkLX7ziGvmqqYEiXvG2C+UOYEFVB8p24EbvS5KwpdZAEZ/yKsCttRuBZN+UImFLrYEifqE7McW/CnNg6WS1Bor4gQJc/EOtgSJ+pwAX39v6Piy7E0q/0wOFRfxIAS6+U+uBwmfrgcIifqYAF+8d3xp48VS4aIpaA0X8TAEu3lFroEjAKMDFM2oNFAk4Bbi4T62BIkFBAS5Np9ZAkaCiAJemUWugSNBRgEvD6mkNXJxTQGbWBgpLy+gYG0P60C56PqVIM1OAS90aaA1cnFNAxqJcyioqASgoLSNjUS6AQlykGSnA5USNtAZmZuVVh/dRZRWVZGblKcBFmpECXI5pYmtgYWlZnW+v77iI+IcCXFzcaA3sGBtDQR1h3TE2xt9VikgNuuvC6cr3Q9Z98PfBsO9HV2vgtQsa7OtOH9qFmKjIWsdioiJJH9rF39WKSA1agTuZh62BR/e5M7Py1IUiEkAKcCfywdTAtKR4BbZIgCnAnURTA0XCigLcKTQ1UCTsKMBDjOsOSDf2njU1UCRsKcBDiNt3QGpqoEhY0zIshDR0B2QtHrQGikjo0Qo8hDTpDkhNDRRxDAV4CGnwDkg9UFjEcbSFEkLqvgMygse6bYHH+8G/F7laA2/9l8JbxAG0Ag8hx98Bmdx2D0+c8hKn53yi1kARB1KAh5i0pHjSzj3jWGvgHrUGijiVAjzUqDVQRKoowEOFHigsIsdRgIcCtQaKSB0U4MFMrYEi0gAFeDDS1EARaQIFeLDR1EARaSKvA9wYEwmsBwqstVd4X5JDaWqgiLjJFyvwPwKbgbY+OJczqTVQRDzg1fLOGJMADAee8U05DqOpgSLiBW9X4I8AfwLa1PcCY8xEYCJAp06dvPy4MLL1Pddet1oDRcRDHq/AjTFXALustdkNvc5a+7S1NtlamxwXF+fpx4WP/bvhjQkw/yqIbOVqDRzxqMJbRNzmzQq8PzDSGHM5EA20NcbMt9aO801pYUatgSLiYx4HuLU2A8gAMMZcAtyj8K6HWgNFxA/UB+5PlYfh08dhzSy1BoqIz/kkwK21a4A1vjhX2FBroIj4mVbgvqapgSLSTBTgvqTWQBFpRgpwX9i/G1ZkQO4/NTVQRJqNAtwbag0UkQBSgHvqp+1VrYFr1BooIgGhAHeXWgNFJEgowJtocU4BS99Zzt0H/0aPiJ0U/cdgOlz3uFoDRSRgtGxsgrc+/4af30zn74fSiTN7uLX8TgYXTGTx9kBXJiJOphV4Y7a+x3nLb6NjxC4WHh7MrMPX8Qut4cgRMrPySEvSClxEAkMBXp8arYFlRzpwdcUDfGZrX6QsLC0LUHEiIgrwE1kLX7xc1Rq4Dy6eys1rz2NH+eETXtoxNiYABYqIuGgPvKaftsNLabD4NtcNObd+BIPu5c7LehITFVnrpTFRkaQP7RKYOkVE0ArcpZHWwKP73JlZeRSWltExNob0oV20/y0iAaUAL8yBpXfAD7kNTg1MS4pXYItIUHFugGtqoIiEOGcGeK2pgTfCpTM0NVBEQo6zArzm1MD258CN78BZvwl0VSIiHnFGgNfRGsiAu6FFq0BXJiLisfAP8JpTA8+8AEY8qqmBIhIWwjfANTVQRMJceAZ4E1sDRURCWdAH+OKcgqbfQFOrNfB0uPol6DZCrYEiEpaCOsAX5xSQsSiXsopKAApKy8hYlAtwYoirNVBEHCaoN4Qzs/Kqw/uosopKMrPyjh3YvxvemADzr4IW0a7WwBGPKLxFJOwF9Qq8vnGthaVlag0UEccL6gDvGBtDQR0h3q9tqWtqoFoDRcTBgnoLJX1ol1pjXCOpZFLLZSysnAL52a7WwBtXKLxFxJGCegVec4zrqXu+ZE70s5xjv4X/VGugiEhQBzhUjXH9ZSGs+YurNfDyl6D7yECXJSIScEEf4AC0S4S+N6g1UESkhtAI8HPHuP4nIiLVgvoipoiI1E8BLiISojwOcGPMmcaY1caYzcaYL40xf/RlYSIi0jBv9sAPA3dbazcYY9oA2caYd621X/moNhERaYDHK3BrbZG1dkPVP+8FNgNqzBYRaSY+2QM3xiQCScA6X5xPREQa53WAG2NOBt4A7rTW/lLHr080xqw3xqwvLi729uNERKSKVwFujInCFd4LrLWL6nqNtfZpa22ytTY5Li7Om48TEZEajLXWszcaY4AXgZ+stXc28T3FwE6PPhDaA7s9fG840vdxjL6L2vR91BYO38dZ1toTVsDeBPhFwEdALnCk6vC91trlHpfY8Oett9Ym++PcoUjfxzH6LmrT91FbOH8fHrcRWmv/BehhkyIiAaI7MUVEQlQoBfjTgS4gyOj7OEbfRW36PmoL2+/D4z1wEREJrFBagYuISA0hEeDGmMuMMXnGmK3GmGmBridQNECsbsaYSGNMjjFmWaBrCTRjTKwx5nVjzJaq/04uDHRNgWKMuavq98m/jTEvG2OiA12TrwV9gBtjIoEngGFAd+A6Y0z3wFYVMEcHiHUDUoA/OPi7qOmPuGbxCDwKrLDWdgV649DvxRgTD0wGkq21PYFI4NrAVuV7QR/gwPnAVmvtdmttOfAKMCrANQWEBoidyBiTAAwHngl0LYFmjGkLDASeBbDWlltrSwNaVGC1AGKMMS2Ak4DCANfjc6EQ4PHA9zV+zsfhoQUaIFbDI8CfOHYzmZP9CigGnq/aUnrGGNM60EUFgrW2AJgNfAcUAXustSsDW5XvhUKA13WzkKNbZxobIOYUxpgrgF3W2uxA1xIkWgB9gSettUnAfsCR14yMMe1w/U29M9ARaG2MGRfYqnwvFAI8Hzizxs8JhOFfhZqqKQPEHKQ/MNIYswPX1tpgY8z8wJYUUPlAvrX26N/KXscV6E50KfCttbbYWlsBLAJ+E+CafC4UAvxz4GxjTGdjTEtcFyKWBrimgKgaIPYssNlaOzfQ9QSatTbDWptgrU3E9d/FKmtt2K2ymspa+wPwvTGmS9WhVMCpT8j6DkgxxpxU9fsmlTC8oOvNI9WahbX2sDFmEpCF60ryc9baLwNcVqD0B/4fkGuM2Vh1zG8DxCQk3QEsqFrsbAduDHA9AWGtXWeMeR3YgKt7K4cwvCNTd2KKiISoUNhCERGROijARURClAJcRCREKcBFREKUAlxEJEQpwEVEQpQCXEQkRCnARURC1P8BbEBDC903/1IAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "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/czJ6Y/Pwr+47kAMpUybYqVbGtDiqVj1OVtia/v9uZjBskfh/VXFW+VKYJ7zeD8HEWds1CoIE4tPnV4amk7aHWVT7GDB4uF4gDwWqQdL9C5JGRXu8K8QgAuAfBaRLCvgsPpb5Sp9d7a5zeDwBYJaH00R4z3lXwMPh+UKThGln3tvNzwVIFDIIzNPAbvJ2kBYhC3ttI9/7t+h7a97ntdydtkffHPS2eeiEKm2ew4iSv7JykfTvB+pgXEfNuAXghi+0yyr2RYkFBh5x4mwAA27p8xG2O52QM4t10K/Z9A3f3LwC42UfnMQBPT/7+NIAPzXRVIYQQB+Z2vwO/190vAsDkz/S1UQghxB3ljgfymNk5AOcAoFvw2ndCCCFm53bfwC+Z2X0AMPmT51cF4O5PuvtZdz/bznkKVSGEELNzuwb8GQCPT/7+OIDPHs50hBBCTMu+X6GY2acAPArgHjO7AOBjAD4O4NNm9rMAXgPw4amuZgBy8pnBhNcoj3qRdg7EamCUXmtY8iWPAy8Upsi/VXFPiBtEja+Dz8io2MDVMk1A/3p5gvZdJblQ+oEHAfOEOJ1zj4fS+YauV+lXYFdIGwCcyNOcM1Gej83A44jl2Ii8UN6u0n3rWnQwUmJvE+5hQc9FuUb75kg9L6KCHlGRhvWa9E+HBQAsk7ww0dyujNJiE72c36cIlptkNSj+0K9SD5Cdij8LdWAEKnKvmGcKAOTE62wc5KypWXtwLiJYXYmqy8fIyvQGWhXc1IB9Dbi7fyT4px+c6UpCCCEOFUViCiFEQ5EBF0KIhiIDLoQQDWW+BR0CMhIKHwmTWUlC6dvTh9Jv3ODi0c4yF4+YQDMKBBMWNs8EJQB4dcjDh9fLdH6vjHhfVom9a1yAahPRLaoGf5wIkAAXLK+QiusR58d83KgIxfkyLW4RiVWs/XpQVOBUKxVvI7Eyqkr/jeH9tJ3xJinSwQRaALgy5qLwLGkB1qt03dG+rZfpmR2SFAZAHGIftTOOFumZZeHuANAOUj3URFiM5nC8lT7XRcaFwqJIxyA6I4BbOE6Q25qVgX0ixRvqZZ6aIkJv4EII0VBkwIUQoqHIgAshREORARdCiIYiAy6EEA1lvl4o7rBhKt+yaNVACAcT05lnCgAUG+kgvsKV7Y0h90JgYb7Pb7+L9mVhyS0WWwvg6pAr71d20pDwtRYvNvDNKk3SH4VBHyMFC6KCAJdz7gnx9f7ppI15zQDAmd71pC3ypDhV8JB+lprg5T7PXNwihQJOdzZoX+aREaVHaBk/L8zjKPL0uE68LLYqvhfv7vC8cCyFwCyFFCLYeTvZC7yQyPkGeJGGTuhBkr4zXtzhe3+kFYTjj4n3VJAnb4eE7kclE8Zjcv8iB7fA8YYVoqlb3D4xu2ej6dM/AHoDF0KIxiIDLoQQDUUGXAghGooMuBBCNJQ5h9IbvEiFAi/Yt/nREEQkCIpMV8upuNLtcJHg9DKvuM3EmO9YeoP23e6mwlQ/mNxyfoy3k5zZD3V5VXqWZ3opCN1nguWRjItEEWUvvXdllwt3a3kqmn5P9zU+tyAMepXc60FQ+Z2lLIjSGNxfpALrwwU/F9vOY6lZLvWB88fp9SLN5x4Juu9tX6LtLL3BS8NUVI7YCPKPn1leT9p6efCMzCAKR7Dq8estHj4eCaFrvTQcP5rzSp6egdymL+HrLa5iBho0iq3UcBlJFQJQUzYzegMXQoiGIgMuhBANRQZcCCEaigy4EEI0lH0NuJk9ZWaXzeyFPW3/3czeMLPnJ//96J2dphBCiJuZxgvlkwB+D8Af3tT+O+7+m7NdzoGaVWKeYQTykWPjIFS1nV6LJW0HgG7ghcASxTPvD4BXYl/OeNGEWZLjRyr/SVKYgIUqA8CpPPWyiTxWTgTtGVl3VCX+oXbqOfOugsvuLfA9ulanHjmRZwnb+/sL7jXBCgJ0jD8KJfi5WMvSub0+5qHtbJ+Zlw4AlMH96xOvlZNBCgIe5s/3jaU8YIVCdsfg55CdWeZtAgBLRTqP9Rb3kHmwe422XxunqQnOtNN1AMCFUertdXo52LdBusfr69zrqepEMfbp2fLApa4m3nfe5c9CxL5v4O7+BQB8J4UQQtw1DvId+C+Y2f+bfMXCnZqFEELcMW7XgP8BgPcAeATARQC/FXU0s3Nm9pyZPTeq+P+aCSGEmJ3bMuDufsndK3evAfwPAO+7Rd8n3f2su59t57MV7BRCCBFzWwbczO7b8+NPAHgh6iuEEOLOsK8Xipl9CsCjAO4xswsAPgbgUTN7BLvpzs8D+LmpruaAjVPFuthOcx4UfT61Ypsot8HHULaVquabOS+k8FKeFkcAgP44VaEvDXnBA5aPYb3k/9cxDibNEuy/e/Vt2vd8kc6ZeVgAQEYSLzzQ4cp95CHDvBBYsQIAuDBK839c6/EcMoOaK+/rJH/Ha8N0XACoyLvIe7q8OAIv6MD7VoEHwRvl8aTtOvGOAICNcXoGjrW2ad+3Cu7Vw/bonwc8pwcrYjAMKqRsk/PdyXkOknaQm2SznL6wxCop0vD6NpfQ3hrwAhIMlocI4M/f+Y303gHA5lbatyA2BADaG8FzNkqfM5KOZbdvSTzyBrMVdNjXgLv7R0jzJ2a6ihBCiENHkZhCCNFQZMCFEKKhyIALIURDmXNBBwBGvvwnHyPtLR6uPtpORYWogj2upQMPgyVfz1PxEAC2d0h47VEuTO6UqdCUZ3wdm2RcABgN0zE2h0H2eEJ0vZyImC+2eEGAXpBWgAmv2yUPNV5tp8rNqztcKD5STB8f8Oo2FzELsu6rQ35PayJMXuyu0b47FRdYbxBhkgmCALA+Svt2gwIEx9p8L9g8rg15CHpZp89IVfN3te1ROufoDEVjGDlb44r3bZNUFjf6XAQ9ssQLjiy30zQGG2SPAWA4Tp/39U3et76e7kUrECCj9B/sKLf6wTM5Iu2k4M2t0Bu4EEI0FBlwIYRoKDLgQgjRUGTAhRCiociACyFEQ5mvF4o7UKbhuNkOaQvU/9ZWkEidQTxeymUeAjve5lsxbqVK8foOV837/dRbJCogUY749WqSQmAjSAKWZeleeLA9zPnnSs29NLod7iGxs5Oq9NWYvwN0l1JPgcs9fr3VDpf6WVqAi+s8jQHbi06Lr+P4UuoqcLHPx10f8L0fjVNvAeaNAQBO1jEi3hEAcO+RtPAGAFy6kYaVDwb8GWFpExDMrSqJ1wN/RODBvWZnrr3E935jRLzI2BwAjIP2q56eo243OLP99MzWwbOe76Tra23xzehc5/tZDEgRmcALJRumtsFu8BQLEXoDF0KIhiIDLoQQDUUGXAghGooMuBBCNJQ5h9I7ME4Fy3w7Fbxam0E16Hb6mUPSH++OO0gFiN4VLkqMhnwrRlXaf3M7qBxNtIrKg8kFolKxns6jLAPxKE/HMDLfqC+I8AcAoyyY84jMIxC8djbSMXZWeUqAq0QoBoCaCF6RSAsisPWDcTdukLzdkQAZ7CcT3qzg13O2b6wNwKtbfI98SK5HRDcAIJH0sDI4F+QMRGcoHwV7QSq0l4FQyK6XDfk6ysH075fbwZm1cTrn1hYft3WD9eXXy0nebwDIy7S9GHBHhnyLiPfEPt4KvYELIURDkQEXQoiGIgMuhBANRQZcCCEayr4G3MweNLO/MbMXzeyrZvbRSftxM3vWzF6a/MkrkwohhLgjTOOFMgbwS+7+ZTNbBfAlM3sWwM8A+Ly7f9zMngDwBIBfueVItcOHqccJI2PJzgG0NlNFN28H6niRtufB5UPlfZBK+uXRILE9EZCzQP1HFPJO2tkcgMCzJCAjajzzCAAAq6f3WAAJE98dI22rbgTvC8EWGfOmmEGkZ94RAFD1yOSCOeTbgccC2c8IUgcjJBqWeVPk/eAcsuXxxwmepePWgZNVFj265AxEZ4id7+gZqTrT73G0PuYZVvSD8Pj1dHKtILK9sxHYpy2SFoSEzAOA7ZANZTkvbsG+R8vdL7r7lyd/3wTwIoAHADwG4OlJt6cBfGimKwshhDgQM30HbmYPAfhuAH8P4F53vwjsGnkApw59dkIIIUKmNuBmtgLgzwD8orvfmOH3zpnZc2b23Mh5jTshhBCzM5UBN7MWdo33H7v7n0+aL5nZfZN/vw/AZfa77v6ku59197Nt42lYhRBCzM40XigG4BMAXnT3397zT88AeHzy98cBfPbwpyeEECJiGi+UDwD4aQBfMbPnJ22/CuDjAD5tZj8L4DUAH953JAOQk88M4pnSuhYMUaVv8VUvKI5AvFBCj4eKf5blJF1BPuSDFDvM64F7Qox7fIyM5KWPvBiMqP9ZkKOhpqkiAk8BLppThTwfBp4ezDMoUNgrnv6D5rjJIi8UMo26FeWFIRsaeQUF3g3MmylwyJlJaao2p09PxM4KwL1FPBg2vNeMYI+Y10o0NzaPIvD0iO/fdG0AkLHnN3hGij5p2wm8Tbb5Qcz76cKzYXBos/Rg+IAXN4nY97S4+98hNHv4wZmuJoQQ4tBQJKYQQjQUGXAhhGgoMuBCCNFQ5lvQoarhW6likXVTYdLJF/wAUFxPlYa8z2N/80Gqgo3Wgkre4yChPxHeqEA3I1EYOxNerQ5C3qePpIdVJJE+aQOA0XIg6JJk9ZEIxsL863x6UQoAMnK9cO9JczQuC68OVZ6gnYl0rJACwO9TJI7WwRPJ2iOhcJbrsfZoHVlwryvySOWRiEmOFrvPu5OLxiBpIaL7RJ7rFnU24A4AxQ5fdLHOxcZsSBZeBiLm+mba5lFOAI7ewIUQoqHIgAshREORARdCiIYiAy6EEA1FBlwIIRrKfL1QAliRB2PhzgAPx664cpt7qip3R1xVHq9w75S6lc4jH00fgh4lpY+8N5hqHinsVTedW9HnAzvxAGGeKbtjRC4LpC3qymo/RF4awV4w7xvmgQDwfY68G+heBJ4+YYg92Tt2P4A43QCDpn+IrtcLrsfCvyPnHdJ1Fq8gIAh5j5ZMukbjjpeC9Q3S/vkgOPfkvGSBvSjW02ypdYebyHx9i1+P2a0N3hfD1JPFic26FXoDF0KIhiIDLoQQDUUGXAghGooMuBBCNJS5ipiO4Ev6PknEW0bl48lnThVUfe6kiaazLk8+nddRbHPa7kdWeN8pfx8AUARK0ZispRXcJraXUVVrJtwEgmA0Bk1vEITHg4huFgg0LF0BABjLjdyhic05wd5bPxWrPDgX0dlie2Ts3gH0Pnl0/6P71w5KxbPLkXuSbfFyhkx08yW+Fzbm+0lF4X6Q15qtI9q3WZ6RWfoG59C2UjtkI26HojNL7/UO33sr0ue63iTh9bdAb+BCCNFQZMCFEKKhyIALIURDkQEXQoiGIgMuhBANxWYN3TzQxcyuAPjm5Md7AFyd28XnzyKvb5HXBmh9TWcR1/cudz95c+NcDfi/u7DZc+5+9q5cfA4s8voWeW2A1td0Fn19e9FXKEII0VBkwIUQoqHcTQP+5F289jxY5PUt8toAra/pLPr6/pW79h24EEKIg6GvUIQQoqHM3YCb2QfN7Otm9rKZPTHv6x82ZvaUmV02sxf2tB03s2fN7KXJn8fu5hwPgpk9aGZ/Y2YvmtlXzeyjk/bGr9HMumb2D2b2fydr+7VJe+PXthczy83sH83sf09+Xpj1mdl5M/uKmT1vZs9N2hZmffsxVwNuZjmA3wfwIwDeC+AjZvbeec7hDvBJAB+8qe0JAJ9394cBfH7yc1MZA/gld/92AO8H8POTe7YIaxwC+AF3/y4AjwD4oJm9H4uxtr18FMCLe35etPV9v7s/ssd1cNHWFzLvN/D3AXjZ3V9x9xGAPwXw2JzncKi4+xcAXLup+TEAT0/+/jSAD81zToeJu1909y9P/r6JXUPwABZgjb7LOwULW5P/HAuwtncwszMAfgzA/9zTvDDrC1j09f0r8zbgDwB4fc/PFyZti8a97n4R2DWAAE7d5fkcCmb2EIDvBvD3WJA1Tr5eeB7AZQDPuvvCrG3C7wL4r/j35acXaX0O4HNm9iUzOzdpW6T13ZJ5V6VnWdDlBtMAzGwFwJ8B+EV3v2FRQvuG4e4VgEfMbA3AZ8zsO+/ylA4NM/txAJfd/Utm9uhdns6d4gPu/qaZnQLwrJl97W5PaJ7M+w38AoAH9/x8BsCbc57DPLhkZvcBwOTPy3d5PgfCzFrYNd5/7O5/PmleqDW6+zqA/4NdPWNR1vYBAP/FzM5j9+vKHzCzP8LirA/u/ubkz8sAPoPdr2kXZn37MW8D/kUAD5vZt5pZG8BPAnhmznOYB88AeHzy98cBfPYuzuVA2O6r9icAvOjuv73nnxq/RjM7OXnzhpn1APwQgK9hAdYGAO7+39z9jLs/hN1n7a/d/aewIOszs2UzW33n7wB+GMALWJD1TcPcA3nM7Eex+71cDuApd/+NuU7gkDGzTwF4FLsZ0C4B+BiA/wXg0wC+BcBrAD7s7jcLnY3AzP4TgL8F8BX82/eov4rd78EbvUYz+w/YFbly7L7MfNrdf93MTqDha7uZyVcov+zuP74o6zOzd2P3rRvY/Tr4T9z9NxZlfdOgSEwhhGgoisQUQoiGIgMuhBANRQZcCCEaigy4EEI0FBlwIYRoKDLgQgjRUGTAhRCiociACyFEQ/kXq61d/zGlqiEAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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+mh4nOuzki6PSLeLekDkj5bXLMutPGkpA9HxPslbZV0o+0PqBtt63ebpIN9n7vWvg9FxNa+qYNda1+ptnvg10l6NiKei4hTkr4naXvLdWhURDwsaelzXbdL2lO83yPppjbr1KSIOBIRjxfv59ULgsvUgTZGz+IzSyeLX6EOtG2R7csl/bWkf+0r7kz7SnS9fX/QdoBfJumlvs+HirKuuTQijki9AJS0ecz1aYTtKyVdK+kRdaSNxe2FJyTNSdoXEZ1pW+Hrkr4oqf8B3V1qX0j6ke3HbM8WZV1q3zm1vaFD6knxTIPJgO31kr4v6fMRccIu23ggLxGxIGmr7XdIesD2e8dcpcbY/rikuYh4zPa2MVdnVK6PiMO2N0vaZ/vpcVeoTW33wA9JuqLv8+WSDrdchzYctb1FkorXuTHXpxbbk+qF93ci4v6iuFNtjIjXJf1MvfGMrrTtekmfsP28ercrP2z72+pO+xQRh4vXOUkPqHebtjPtG6TtAH9U0jW2r7I9JekWSXtbrkMb9kraUbzfIenBMdalFve62ndLOhgRX+v7UfZttH1J0fOW7T+R9BFJT6sDbZOkiPi7iLg8Iq5U79/aTyLib9SR9tleZ3vD4ntJH5V0QB1p3zBaX8hj+2Pq3ZebkLQ7Ir7aagUaZvteSdvUewLaUUl3Svp3SfdJ+lNJL0q6OSKy3MDS9l9I+i9JT+rsfdQvqXcfPOs22n6feoNcE+p1Zu6LiL+3fZEyb9tSxS2UL0TEx7vSPttXq9frlnq3g78bEV/tSvuGwUpMAMgUKzEBIFMEOABkigAHgEwR4ACQKQIcADJFgANApghwAMgUAQ4Amfp/1Ut+qRi9fuUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "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+XC47mfc13faYX7yodTH3SZXxAmE7HvUlZJnlfkwBILXVL9zVW8ljV9Dixy4b0Qu1t8SXj+UO8qJUYRSPd2GNL/9rn8XjYx3+cS15k8sWO+f6ErbK5/yiTwBw9ed/Sce+/+Ry157ayyWa3aN96Vv5Vn6pd1T7EtKubj7vRIQMMkz032/8auqCc299zbU/c//Z1Ket1j8fJ53r9/kFgNbHeOGuHFEapy/xJb4AkPyZLzVuO7+D+hS94een7GLeX/ZY0R26EELEBCV0IYSICUroQggRE5TQhRAiJiihCyFETFBCF0KImDD4skVCdowvIWsv4xXf0OrLt0Zv5hXsMqTw3XUXvUx9Hv/+BXSsebYvE0u283kniVoum+Q9BbtG+2PFFbzyXsHqcjrW47c1RPWkPdRnVc1s116yjkv5Ri6vc+2pb5JyhgA6b2507fX7KqgPkxkCQPdY/xgV7+YyyMw436d9GpeJFu70ZX7Jq/39AYDH7ryQv9/HfAlpVy2X5TJ5op0fUbW02T9+oYRLS9HJ7wWLtxC54/X11OdX31js2tsv4RU2U7v97Rx4gUsTOxby41dS6UsN9xOZMQDYVN9esIkfI1Y1NBfRz/dY0R26EELEBCV0IYSICUroQggRE5TQhRAiJiihCyFETDiqysXM7gZwOYCGEMLcvK0SwEMApgLYCeDjIQS+jH6YnKGgxVcjFJBef4l5vHdh99t+gaAD1cffm++HL55Hx1If5EVzrNaXi5TWcKVN6ylkfj3887WnypfGVN3LFSbT/+5NOvbCutNc+9onT6c+Z1z4jmvf0sR7Rda9Nc61Jy7g8Um+5vskZ/OCR+ktXLHy+5940bX/7LvnUx+Yf56esnQ3daldO9W1N4zkip7EB/gUss2+iiPRE6HiImqoMcW8h2vld311xf4zeHpomx1R2Gy/fyz2H+KqK5zh71OihhcIC0TY1FvK1WKjxrfQscwq/zilI9JJZmG7a0+u4yqXikp/Do11vKfwsXIsd+j3ALj0XbbbATwbQpgJ4Nn8ayGEEEPIURN6COElAAffZb4SwL353+8FcNWJnZYQQojjpb/P0MeHEGoBIP/T/z+yEEKIQeN9XxQ1sxvNbI2Zrcm1+c+bhBBCDJz+JvR6M6sCgPxP2tIjhLAihFAdQqhOlPGFAiGEEAOjvwn9cQA35H+/AcBjJ2Y6Qggh+suxyBYfALAUwBgz2wvgywC+BuBHZvZJALsBfOyYtmZArtCXFGUn+5K0nq6IPpLkvcq2893qLvd9eqfyIlclL/KiOT3LfVllawH3SWR8iVY2xeVWqTo/Dnsu4QWU9qzlEsRUq/9Zft5l66nPvo/7/V0v+glvFvn8Q74uL3c2l6N2kl6RJetJRTEAp16xlY498W1fnnhoIdejpcf652P25pHU529+fL9rv/O711GfGVe9Tcdq7/LloHUX8wJT2fH++XDolZOoT/dy32d29Q6+nVsq6NikFZtc+6pH5lOf3kV+IbLyX/DrqPvKJtfetoMfo861XELaNcuXdlpE0byqn/qyymV/y4v9PfltXyKdOyOiGNoxctSEHkK4lgzxMnFCCCEGHX1TVAghYoISuhBCxAQldCGEiAlK6EIIERMGtQWdpXJInuSrB9Kv+avZHRP4CvOURTWufVduAvXJVva69tP+mX/pafq9G+nYytcXuPZQ6m8HAAoO+oqVZCdvoxbIR28RiScA5LZwhUDvFF/Vs/ssHoen9v3KtZ9691/wOZBiUYmNfmE1AEj6nQXRPpsXmHpj9Qw6VniRX1wt0cVP/6mf9s+tm1e/Sn1u+84nXXv4YBP12fz0TDqW+YCveijeRQIEAHN8tUjnRH4+pir8uLZ8YzL1mfcDroZa862Frr39fK7OQZu/T5nFXIk04nm/n2T6Aq6g6kjxa8JaUq69pJbf97b8UZNr/8V3eLG/lmV+CzrL8qJrx4ru0IUQIiYooQshRExQQhdCiJighC6EEDFBCV0IIWKCEroQQsSEQZUtIpMAtvsldBPL3t0UqQ/b6EuTAGDvrye69qK2iJ6L5u9y7t+5XO+1FYvomM3xZXnl67m0rH2SL0djhcsAIDnGl5Zl9vGSxAURH9fjVvrzu3XbW9TnQ39MZHmf4rFLv+7LxHK+QgwA0DHSl6pVVPLttO7n50lXs7+vZVt54bdlL/u9Q//lzz9Bfb5y13+69ju/dj31gfFjniEyza5SftmWr/L7dk66bA/1qXvKlyee9eVV1Of5u86mY+1X+jLR9DreU7T4gB+HQ9Vc6ti6yJfeVq7kktgZ12+nY1t/6RdDq1xWS30y9/lFzxov5D1XQ4t/3lmvZItCCCHyKKELIURMUEIXQoiYoIQuhBAxQQldCCFighK6EELEBAuBy6ZONEWTJ4eJn/+cO5Yd7cuTikq5/Cf1ui+DGrm8jvo0vDGeT5CwZNkGOrb24TP8gYiw9p7b4toznRFaPqJoKlmfpi7tU3mlukkzGlx71wO89+Rln3/RtT98/1Lqw+hewCWIPURmiCK+P8U7uUw0QZRvn75uJfV59AsXufaz/on3T330iXNce8/kCAlbRIXNwoP+WPGcJurTuteX7I3+Db93W3KTv08vf8/vBwsAPR+JmEO9L1UtruHn95Slu1z7viem8O3M8itIJkfweFc+za+X6pvecO0vP8xly50L/GqnuR4e7/J1fh/SzsX8mnjnmr9fG0Kopn+QR3foQggRE5TQhRAiJiihCyFETFBCF0KImKCELoQQMeGoKhczuxvA5QAaQghz87Z/APBpAI35P/tSCIFLBvIUTZ8YJtxxkzuWa/VXwIsbeCGi7Kn+qnDJKl6wKv3hetdev30M9SndyZUIHRP8Qlvnn8OLXL306hzXnivlKo6iej8OmQgFRaqWF58at8iPQ+NargJKtfpSmyzfDP7yD3/u2r/+ysXUx7r9+4zS3fw4FF+wn44d3Fbp2pMdvBhS7wQ/rnaQKzVGr/ffr/x6vz8pAHT18vO7/k3/WIzcRl0w4up9rn3f67zPbk+Ffw6HJM8NJWO5IqP4KV9pM/G6HdRnyyvTXHv3uIheqI2k0N5Uv2cnAGQ7eLxTjf6x7RnFr0sk/BgVkusVAErmHXLtbZt5gbntt916wlQu9wC41LF/I4SwIP/vqMlcCCHE+8tRE3oI4SUAfm1bIYQQvzMM5Bn6Z83st2Z2t5nx/ysIIYQYFPqb0O8CcAqABQBqAfwr+0Mzu9HM1pjZmmwLf+4mhBBiYPQroYcQ6kMI2RBCDsD3ASyO+NsVIYTqEEJ1cgRfrBRCCDEw+pXQzazqiJcfBcCLnQghhBgUjkW2+ACApQDGAKgH8OX86wXoK0G1E8BnQgi88d7/vlcjgMNVeMYA4Fqz4YPi0IfioBgcRnF4bwymhBDGHs1pUKst/p8Nm605Fl1l3FEc+lAcFIPDKA79j4G+KSqEEDFBCV0IIWLCUCb0FUO47d8lFIc+FAfF4DCKQz9jMGTP0IUQQpxY9MhFCCFiwpAkdDO71My2mNk2M7t9KOYwFOTLJDSY2YYjbJVm9oyZvZ3/GesyCmY22cyeN7NNZvaWmd2Stw+3OBSb2etmtj4fh3/M24dVHADAzJJm9oaZPZF/PRxjsNPM3jSzdWa2Jm877jgMekI3sySAbwP4MIDTAVxrZqcP9jyGiHvw3sqVtwN4NoQwE8Cz+ddxphfArSGE0wCcDeCm/PEfbnHIAPhQCGE++r7TcamZnY3hFwcAuAXApiNeD8cYAMCyfPXaw3LF447DUNyhLwawLYSwPYTQDeBBAFcOwTwGHVK58koA9+Z/vxfAVYM5p8EmhFAbQvhN/vdW9F3IEzH84hBCCG35l6n8v4BhFgczmwTgMgA/OMI8rGIQwXHHYSgS+kQAe454vTdvG66MP/wt2/zPcUM8n0HDzKYCWAjg1xiGccg/algHoAHAMyGE4RiHbwK4DcCRXTaGWwyAvg/zp81srZndmLcddxx4W433D6+ti6Q2wwwzKwPwEwCfCyG0mPHuQXElhJAFsMDMKgA8YmZzh3hKg4qZHe6EttbMlg7xdIaaJSGEfWY2DsAzZra5P28yFHfoewFMPuL1JAB+36zhQf3hYmf5nw1DPJ/3HTNLoS+Z/1cI4ad587CLw2FCCE0AXkDf+spwisMSAFeY2U70PXr9kJndh+EVAwBACGFf/mcDgEfQ92j6uOMwFAl9NYCZZjbNzAoBXAPg8SGYx+8KjwO4If/7DQAeG8K5vO9Y3634fwDYFEL4+hFDwy0OY/N35jCzNIDlADZjGMUhhPDFEMKkEMJU9OWB50IIn8AwigEAmFmpmZUf/h3AxeirYHvccRiSLxaZ2UfQ9+wsCeDuEMJXB30SQwCpXPkogB8BOBnAbgAfCyHEtuWfmZ0H4GUAb+J/n5t+CX3P0YdTHOahb6Erib4bqx+FEL5iZqMxjOJwmPwjly+EEC4fbjEws+nouysH+h6D3x9C+Gp/4qBvigohREzQN0WFECImKKELIURMUEIXQoiYoIQuhBAxQQldCCFighK6EELEBCV0IYSICUroQggRE/4bdYmDYkl2gyoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "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+XC47mfc13faYX7yodTH3SZXxAmE7HvUlZJnlfkwBILXVL9zVW8ljV9Dixy4b0Qu1t8SXj+UO8qJUYRSPd2GNL/9rn8XjYx3+cS15k8sWO+f6ErbK5/yiTwBw9ed/Sce+/+Ry157ayyWa3aN96Vv5Vn6pd1T7EtKubj7vRIQMMkz032/8auqCc299zbU/c//Z1Ket1j8fJ53r9/kFgNbHeOGuHFEapy/xJb4AkPyZLzVuO7+D+hS94een7GLeX/ZY0R26EELEBCV0IYSICUroQggRE5TQhRAiJiihCyFETFBCF0KImDD4skVCdowvIWsv4xXf0OrLt0Zv5hXsMqTw3XUXvUx9Hv/+BXSsebYvE0u283kniVoum+Q9BbtG+2PFFbzyXsHqcjrW47c1RPWkPdRnVc1s116yjkv5Ri6vc+2pb5JyhgA6b2507fX7KqgPkxkCQPdY/xgV7+YyyMw436d9GpeJFu70ZX7Jq/39AYDH7ryQv9/HfAlpVy2X5TJ5op0fUbW02T9+oYRLS9HJ7wWLtxC54/X11OdX31js2tsv4RU2U7v97Rx4gUsTOxby41dS6UsN9xOZMQDYVN9esIkfI1Y1NBfRz/dY0R26EELEBCV0IYSICUroQggRE5TQhRAiJiihCyFETDiqysXM7gZwOYCGEMLcvK0SwEMApgLYCeDjIQS+jH6YnKGgxVcjFJBef4l5vHdh99t+gaAD1cffm++HL55Hx1If5EVzrNaXi5TWcKVN6ylkfj3887WnypfGVN3LFSbT/+5NOvbCutNc+9onT6c+Z1z4jmvf0sR7Rda9Nc61Jy7g8Um+5vskZ/OCR+ktXLHy+5940bX/7LvnUx+Yf56esnQ3daldO9W1N4zkip7EB/gUss2+iiPRE6HiImqoMcW8h2vld311xf4zeHpomx1R2Gy/fyz2H+KqK5zh71OihhcIC0TY1FvK1WKjxrfQscwq/zilI9JJZmG7a0+u4yqXikp/Do11vKfwsXIsd+j3ALj0XbbbATwbQpgJ4Nn8ayGEEEPIURN6COElAAffZb4SwL353+8FcNWJnZYQQojjpb/P0MeHEGoBIP/T/z+yEEKIQeN9XxQ1sxvNbI2Zrcm1+c+bhBBCDJz+JvR6M6sCgPxP2tIjhLAihFAdQqhOlPGFAiGEEAOjvwn9cQA35H+/AcBjJ2Y6Qggh+suxyBYfALAUwBgz2wvgywC+BuBHZvZJALsBfOyYtmZArtCXFGUn+5K0nq6IPpLkvcq2893qLvd9eqfyIlclL/KiOT3LfVllawH3SWR8iVY2xeVWqTo/Dnsu4QWU9qzlEsRUq/9Zft5l66nPvo/7/V0v+glvFvn8Q74uL3c2l6N2kl6RJetJRTEAp16xlY498W1fnnhoIdejpcf652P25pHU529+fL9rv/O711GfGVe9Tcdq7/LloHUX8wJT2fH++XDolZOoT/dy32d29Q6+nVsq6NikFZtc+6pH5lOf3kV+IbLyX/DrqPvKJtfetoMfo861XELaNcuXdlpE0byqn/qyymV/y4v9PfltXyKdOyOiGNoxctSEHkK4lgzxMnFCCCEGHX1TVAghYoISuhBCxAQldCGEiAlK6EIIERMGtQWdpXJInuSrB9Kv+avZHRP4CvOURTWufVduAvXJVva69tP+mX/pafq9G+nYytcXuPZQ6m8HAAoO+oqVZCdvoxbIR28RiScA5LZwhUDvFF/Vs/ssHoen9v3KtZ9691/wOZBiUYmNfmE1AEj6nQXRPpsXmHpj9Qw6VniRX1wt0cVP/6mf9s+tm1e/Sn1u+84nXXv4YBP12fz0TDqW+YCveijeRQIEAHN8tUjnRH4+pir8uLZ8YzL1mfcDroZa862Frr39fK7OQZu/T5nFXIk04nm/n2T6Aq6g6kjxa8JaUq69pJbf97b8UZNr/8V3eLG/lmV+CzrL8qJrx4ru0IUQIiYooQshRExQQhdCiJighC6EEDFBCV0IIWKCEroQQsSEQZUtIpMAtvsldBPL3t0UqQ/b6EuTAGDvrye69qK2iJ6L5u9y7t+5XO+1FYvomM3xZXnl67m0rH2SL0djhcsAIDnGl5Zl9vGSxAURH9fjVvrzu3XbW9TnQ39MZHmf4rFLv+7LxHK+QgwA0DHSl6pVVPLttO7n50lXs7+vZVt54bdlL/u9Q//lzz9Bfb5y13+69ju/dj31gfFjniEyza5SftmWr/L7dk66bA/1qXvKlyee9eVV1Of5u86mY+1X+jLR9DreU7T4gB+HQ9Vc6ti6yJfeVq7kktgZ12+nY1t/6RdDq1xWS30y9/lFzxov5D1XQ4t/3lmvZItCCCHyKKELIURMUEIXQoiYoIQuhBAxQQldCCFighK6EELEBAuBy6ZONEWTJ4eJn/+cO5Yd7cuTikq5/Cf1ui+DGrm8jvo0vDGeT5CwZNkGOrb24TP8gYiw9p7b4toznRFaPqJoKlmfpi7tU3mlukkzGlx71wO89+Rln3/RtT98/1Lqw+hewCWIPURmiCK+P8U7uUw0QZRvn75uJfV59AsXufaz/on3T330iXNce8/kCAlbRIXNwoP+WPGcJurTuteX7I3+Db93W3KTv08vf8/vBwsAPR+JmEO9L1UtruHn95Slu1z7viem8O3M8itIJkfweFc+za+X6pvecO0vP8xly50L/GqnuR4e7/J1fh/SzsX8mnjnmr9fG0Kopn+QR3foQggRE5TQhRAiJiihCyFETFBCF0KImKCELoQQMeGoKhczuxvA5QAaQghz87Z/APBpAI35P/tSCIFLBvIUTZ8YJtxxkzuWa/VXwIsbeCGi7Kn+qnDJKl6wKv3hetdev30M9SndyZUIHRP8Qlvnn8OLXL306hzXnivlKo6iej8OmQgFRaqWF58at8iPQ+NargJKtfpSmyzfDP7yD3/u2r/+ysXUx7r9+4zS3fw4FF+wn44d3Fbp2pMdvBhS7wQ/rnaQKzVGr/ffr/x6vz8pAHT18vO7/k3/WIzcRl0w4up9rn3f67zPbk+Ffw6HJM8NJWO5IqP4KV9pM/G6HdRnyyvTXHv3uIheqI2k0N5Uv2cnAGQ7eLxTjf6x7RnFr0sk/BgVkusVAErmHXLtbZt5gbntt916wlQu9wC41LF/I4SwIP/vqMlcCCHE+8tRE3oI4SUAfm1bIYQQvzMM5Bn6Z83st2Z2t5nx/ysIIYQYFPqb0O8CcAqABQBqAfwr+0Mzu9HM1pjZmmwLf+4mhBBiYPQroYcQ6kMI2RBCDsD3ASyO+NsVIYTqEEJ1cgRfrBRCCDEw+pXQzazqiJcfBcCLnQghhBgUjkW2+ACApQDGAKgH8OX86wXoK0G1E8BnQgi88d7/vlcjgMNVeMYA4Fqz4YPi0IfioBgcRnF4bwymhBDGHs1pUKst/p8Nm605Fl1l3FEc+lAcFIPDKA79j4G+KSqEEDFBCV0IIWLCUCb0FUO47d8lFIc+FAfF4DCKQz9jMGTP0IUQQpxY9MhFCCFiwpAkdDO71My2mNk2M7t9KOYwFOTLJDSY2YYjbJVm9oyZvZ3/GesyCmY22cyeN7NNZvaWmd2Stw+3OBSb2etmtj4fh3/M24dVHADAzJJm9oaZPZF/PRxjsNPM3jSzdWa2Jm877jgMekI3sySAbwP4MIDTAVxrZqcP9jyGiHvw3sqVtwN4NoQwE8Cz+ddxphfArSGE0wCcDeCm/PEfbnHIAPhQCGE++r7TcamZnY3hFwcAuAXApiNeD8cYAMCyfPXaw3LF447DUNyhLwawLYSwPYTQDeBBAFcOwTwGHVK58koA9+Z/vxfAVYM5p8EmhFAbQvhN/vdW9F3IEzH84hBCCG35l6n8v4BhFgczmwTgMgA/OMI8rGIQwXHHYSgS+kQAe454vTdvG66MP/wt2/zPcUM8n0HDzKYCWAjg1xiGccg/algHoAHAMyGE4RiHbwK4DcCRXTaGWwyAvg/zp81srZndmLcddxx4W433D6+ti6Q2wwwzKwPwEwCfCyG0mPHuQXElhJAFsMDMKgA8YmZzh3hKg4qZHe6EttbMlg7xdIaaJSGEfWY2DsAzZra5P28yFHfoewFMPuL1JAB+36zhQf3hYmf5nw1DPJ/3HTNLoS+Z/1cI4ad587CLw2FCCE0AXkDf+spwisMSAFeY2U70PXr9kJndh+EVAwBACGFf/mcDgEfQ92j6uOMwFAl9NYCZZjbNzAoBXAPg8SGYx+8KjwO4If/7DQAeG8K5vO9Y3634fwDYFEL4+hFDwy0OY/N35jCzNIDlADZjGMUhhPDFEMKkEMJU9OWB50IIn8AwigEAmFmpmZUf/h3AxeirYHvccRiSLxaZ2UfQ9+wsCeDuEMJXB30SQwCpXPkogB8BOBnAbgAfCyHEtuWfmZ0H4GUAb+J/n5t+CX3P0YdTHOahb6Erib4bqx+FEL5iZqMxjOJwmPwjly+EEC4fbjEws+nouysH+h6D3x9C+Gp/4qBvigohREzQN0WFECImKKELIURMUEIXQoiYoIQuhBAxQQldCCFighK6EELEBCV0IYSICUroQggRE/4bdYmDYkl2gyoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "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+f7q3eyuayci1MM3lkzmpjk5JERpp66IP9NOWzmtvn4v6/c18os39rO1upWUmHA+PT+PZTOzKErXhVdE/JEKX86q32vZeKCZh9dXsq68kfBQD/cum8pVUzJIi4twOp6InAMVvgxZTUsXX3l6K6UHW/AYuKgwhR/dVEJucrTT0URkCFT4ck76+r1sqW5l1bZafrfxEJGhHm5bcAGfv6xIY/wiLqfCl/O2v6GDn7+xn1Xba4kI9fC5S4v4uyuLdMUtEZfSgVdy3orS4/jlrbN44fMXMz07gZ+9vpe/fXwzRzt7nI4mIudIhS9DMiM3kf9euYBb5+Xy6q56vvDE+3T39jsdS0TOgQpfhszjMfzgEyXcu2wqGyqbufYXf+SN3fVOxxKRIdIFUOSc3b4wn7S4CL7/ym5WPFbK9OwE7lhcwJwLkshJ0mweEbfSTls5b8d7+vn3tft4cuMhWroGTs9ww8wsfnRTCZFhIQ6nEwlOmqUjo+pEXz/lRzr4703VPLHxEBPSY1k+L48ViwqcjiYSdFxxTVtjTCHwTSDBWnvzWK1XRl9EaAglOYmU5CRyYX4y//nHSr63ahdZCZEsnZ7pdDwR8RnSTltjzKPGmAZjzM5Tli8xxpQbY/YbY+4+23NYayuttSuGE1bc78ZZ2fzqtjkAfO6J91lX3uBwIhE5aaizdH4LLBm8wBgTAjwALAWKgVuNMcXGmOnGmFWnfKSPaGpxtdzkaF7+4iIiQj3c8Vgp332pjCc2HqTp2Amno4kEtSGP4Rtj8oFV1tppvvsLgO9aa6/x3b8HwFr7g494nufONqRjjFkJrATIy8ubc/DgwSHlE/fZVdvOt1/cyc7aNrp7vSREhXHT7Bw+szCf3OQoHa0rMgpGaww/G6gedL8GmH+WECnA94FZxph7zvSLwVr7MPAwDOy0HUY+cVhxVjzPfW4hXq9lz5EOfvpaOf+1oYpH3zlASkw4P/jEdD5WnKHiFxkjwyn80/0vPWNBW2ubgbuGsT7xUx6PoTgrnkduv5Cqpk7e3t/E4xsOsvLxzaxYVMA/X1fsdESRoDCcwq8BcgfdzwFqhxdngDHmeuD6oqKikXg6cZH81BjyU2O4YWYW9768i1+/fYD69m6+c/1UnXtfZJQNZww/FNgLXAkcBjYBn7LWlo1UOM3DD2x9/V5+8tpeHnqrAq+FcfGRfLwkk1l5iSydlkmIR0M9Iudq2AdeGWOeAi4DUoF64DvW2l8bY64F7gdCgEettd8fqdCgwg8W++o7eGNPA2vKjrCrtp0TfV6K0mO5fWE+n5qXp+IXOQc60lb8htdreWHrYX7+xj4ONndRkBrD1VMzuG3+BeQkaWaPyEfxu8IfNIZ/5759+5yOIw6w1vLQ+krW7mlgU9VRvBYSo8O4Y1EB15VkkZ8a43REEVfyu8I/SVv4AlDZeIx39jexansdGw8cBWDyuDiunjqOq4szmDQujrAQnelbBFT4EkCqj3bx6q56Vu+oY8uhFrwWYsJD+KsF+VxXksm07ASnI4o4SoUvAammpYvSqhae3VzNuxXNWAtTs+JZVJTK+LRYlkwfR3ykLrouwcXvCl9j+HKuWjp7eHl7Lc+W1lBe30FPn5dQjyE/NYZFRancOCubmbmJTscUGXV+V/gnaQtfzoe1ls0HW3htVz07a9vYfLAFrxc+e3E+STHhfHx6JrnJujKXBCYVvgS1ls4e7vyvUrZUt9LvtUSGeZg0Lp5b5ubwF7OyiQ7XlT4lcKjwRRjY8q8+epxH3q5k/d5Gqpq7CAsxXDs9k4vHpzI5M46JGXG6PKP4NRW+yCn6vZZ3K5p4aWsta8qO0N7dB0B4iIeLxqdwx6IC5uYnaetf/I7fFb522spY8not+xqOUdl4jM0HW/jdxoN093qJCPVw6cQ0Fo5PYX5hClMy452OKvKR/K7wT9IWvjihrauXrTWtvLm7nj+UHaG+/QShHsPKSwpZ6Bv6SY3VmT3FnVT4IufJ67XUth3nWy/sZF15IwCxEaF85aoJXDs9k6zEKIcTinyYCl9kBFQf7WJLdSv3v76XysZOosNDWFCYQklOInddVkhEqHb2ivNU+CIjbG99Bw+uq2BXbTvl9R0kRYcxLTuBBeNTKEyNZXZeIunxkU7HlCDkd4WvnbbiT9aVN7B6Rx3ba9rYc6QDgFCP4V8/MZ3ZeYkUpsbi0Tn9ZYz4XeGfpC188TdHO3vYWNnM/1u1i7q2bgBSYsK5eU4OKxYXkB6nrX4ZXSp8kTHW2+9le00rFQ2dvLKjjrf2DuzwLUqPZdmMLBZPSGVWXpLDKSUQqfBFHLartp215Q28vrueLYdagYFz+hemxTC/IIVZeYkUpsUSG6EDvWR4VPgiLtLa1cNzm2t4fXc9NS3HqWk5/sHXLkiJJjcpmk/OzeG6kixdz1fOmQpfxMUONHVSfqSdisZOdtS0saW6hfr2EyRFhzEzN5HJmfEsnTaOaVkJ2vkrH+lsha+/H0UcVpAaQ8Gga/R29/bz+u561u5ppKy2jbf2NvLgugqSY8LJTY7mry+6gGUzs3RZRzlnrtzC17RMkT9r6OhmTVk9ZYfb2FR1lIrGTsJDPOSlRFOSncDsC5KYkhlHVmIUmQk68jfYaUhHJEB09/bz6q56ymrbqGjoZGt1K03HTnzw9WumZjA1K4HxabHMykskMyESYzQMFEw0pCMSICLDQlg2I4tlM7KAgXP8HzraRWVTJ69sr+ON3fWsKav/4PHpcRHMyE3k4vEpLJ2eSYaO/g1q2sIXCTDdvf2UH+lgy6EWtla3sqmqhcOtAzOBkqLDuHxSOuPTY5lXkMycvCTtCA4w2sIXCSKRYSHMyE1khu+i7V6vpay2nY0HmimrbWfd3kae33IYgOLMeGbkJnL5pDQWT0gjKlwngAtkKnyRAOfxGKbnJDA9J+GDZUc7e1i9o45nN9ewekcdT713CI8ZOBJ4WtbAY68rySItTuf9DyQa0hEJcn39Xtbva2RrdRs7Dw98NHQM7AiOCQ/h4qJUrpySzpTMeLITo0jRxV9cTUM6InJGoSEerpicwRWTMz5YVtF4jNXb66jv6GbV9jpe3TWwIzgsxHD7gnwun5zOgsIUjf/7GW3hi8hZeb2WA82d7Ks/xu931rFqex39Xkt6XASXTkxjSmY8xVnxTMmMJyEqzOm4Qc/v5uHrwCsR9+ro7uXNPQ28sOUw22raONrZ88HX5hcks6golUsmplGSk6BjABzgd4V/krbwRdyv+mgX22paee/AUd7c0/DByeDyU6JZNCGV+QUpXDIhjYRobf2PBRW+iIyZtq5e/lBWx+odRyitOkpnTz/pcRGsvKSQovRYZuUlaehnFKnwRcQRff1e3qlo5jsv7qSquQuA8BAPf7OogL+7sojocM0bGWkqfBFxXPOxE5TXd/DMpmpe2FpLRKiHovRYlk4b98HF35Niwp2O6fc0LVNEHJcSG8HC2AgWjk9l2cwsXt/dwI6aNu57de8Hj8lJiuKuS8dzy9xcwkN1+ueRpi18EXGMtZaaluPsa+igsrGTJzceorKpk4hQD/MKklkybRxXF4/TEb/nQEM6IuIX+r2W9XsbeXNPA2/tbeTQ0YFx/6L0WC4qTGZGTiKFaTFMzUogMkzn/TkdFb6I+B2v17L9cBt/qmzmT5XNbDowMOMHIDzUw/yCZJZfmMfSaeN0xO8gKnwR8Xt9/V6qW46zr76DjQeOsmp7LfXtJ5iRk8DNcwYu+q6dvip8EQlAXq/lp6/tZfWOOiqbOoGBoZ9vLJnMvPzkoD3QS4UvIgHLWsvuug5e3HqYxzZU0d3rxWMYOL1DdgKTM+O5bFJa0Mz597vC17l0ROR8HDvRx6aqo6zb08A7Fc1UNh7DayE+MpSSnEQuzE9mXkEy03MSiI0IzF8Aflf4J2kLX0SGo7u3n9d31/PCllpqW4+z+0g7JytvQWEKP7llBlmJUc6GHGEqfBERBs7zs6GymV117fxqXQV9Xi+FabHcvWQyVxVnfPQT+AEVvojIKcqPdPDK9lr+d+thqo8e55e3zuK6kky/P6WzCl9E5AwONXfxmd+8R2VTJ7PyEvnrBRdwXUkWYSH+eWoHFb6IyFl09fTx0FuV/NeGKlq6eslPieaWC3O5ZW4uqX52DV8VvojIEPT1e1lTVs+/r93P7rp20uIi+NrHJnLjrGy/OZXD2QrfP/9mEREZBaEhHj5eksnvv7yYVV9aREZ8BHc/v4O//vV71LYedzresKnwRUROY1p2Ai9+YRF3L51M6cGjXH7fOrZVtzoda1hU+CIiZxDiMdx16Xje+ofLSYkJ58tPb+GwH2/pq/BFRD5CbnI09y+fRV1bN5fft457nt9BV0+f07HOmQpfRGQI5hUk8+bXL2PptHE89d4hbnloA/1e9056OR0VvojIEGUnRvHz5bP4h2smsfNwO1/5761U+y7S4g9U+CIi5+hzl45n8YRUXt5Wy5U/eYsH1u53OtKQBObp4kRERpHHY3h8xXxqW4/z7RfL+Lc15bR09vDJubkUpMa49gLs7kwlIuIHshKj+OlfzuCGmVk88vYBrrl/PTc+8A7t3b1ORzstFb6IyDDER4bx8+Wz+P2XF/OlK4rYVdfOJT9ey4tbDzsd7f9Q4YuIjIApmfF87epJvPTFi7kgJYZ/fG47++o7nI71IWNW+MaYG40x/2mMedEYc/VYrVdEZCyV5CRy380l9HktS3/+Ryoajzkd6QNDKnxjzKPGmAZjzM5Tli8xxpQbY/YbY+4+23NYa1+w1t4JfAb4y/NOLCLichMy4njurgV4reXvn9nGjpo2pyMBQ9/C/y2wZPACY0wI8ACwFCgGbjXGFBtjphtjVp3ykT7oW7/l+z4RkYA1Ky+Jf7p2CpUNx7jxP97ha89so76929FMQz49sjEmH1hlrZ3mu78A+K619hrf/XsArLU/OMP3G+CHwGvW2tfPsp6VwEqAvLy8OQcPHhzyP0ZExG2aj53gB7/fw3ObawjxGH72lzNZNiNr1NY3WqdHzgaqB92v8S07ky8BVwE3G2PuOtODrLUPW2vnWmvnpqWlDSOeiIjzUmIjuO+TM/jDVxYzOy+Rrz+7zbGzbg6n8E934ccz/rlgrf2FtXaOtfYua+2vhrFeERG/M3lcPA//1VwSosL44e/3OJJhOIVfA+QOup8D1A4vzgBjzPXGmIfb2tyxo0NEZCQkxYRz16Xj2VDZzCvb6xjrKw4Op/A3AROMMQXGmHBgOfDSSISy1r5srV2ZkJAwEk8nIuIan56fR3ZiFF948n1+9IfyMV33UKdlPgVsACYZY2qMMSustX3AF4E1wG7gGWtt2ehFFRHxf5FhIbz4xYuZmZvIr96qYHtN65ite0iFb6291Vqbaa0Ns9bmWGt/7Vu+2lo70Vo73lr7/dGNKiISGFJjI/jWx6cA8Pkn3h+z9bry1AoawxeRQDc3P5l7l02lpuX4mJ1T35WFrzF8EQkGC8enAPDQ+gp6+ryjvj5XFr6ISDAoSo9lenYCv/vTIb7z0ujvAlXhi4g4xBjD859fyNXFGazaVjvqW/muLHyN4YtIsAgL8XDjrGw6TvSN+owdVxa+xvBFJJjMK0gmPMTDL98c3WvjurLwRUSCSWpsBMvn5bJ+X+OoXh5RhS8i4gJLp2ViLTywdvS28l1Z+BrDF5Fgs2B8CvMKknmrvHHU1uHKwtcYvogEo8smpbHnSAcHmjpH5fldWfgiIsHo6uJxAPzk1dE5qZoKX0TEJYrSY/nE7Gz+uK+Jfu/Inzo5dMSfUUREztu10zLBwrETfSREhY3oc6vwRURc5KriDK4qzhiV53blkI5m6YiIjDxXFr5m6YiIjDxXFr6IiIw8Fb6ISJBQ4YuIBAkVvohIkFDhi4gECVcWvqZlioiMPGPtyB++O1KMMY3AwfP89lSgaQTjjAQ3ZgJ35lKmoXNjLmUaupHOdYG1Nu10X3B14Q+HMabUWjvX6RyDuTETuDOXMg2dG3Mp09CNZS5XDumIiMjIU+GLiASJQC78h50OcBpuzATuzKVMQ+fGXMo0dGOWK2DH8EVE5MMCeQtfREQGUeGLiASJgCt8Y8wSY0y5MWa/MebuMV73o8aYBmPMzkHLko0xrxlj9vk+Jw362j2+nOXGmGtGKVOuMWatMWa3MabMGPNlp3MZYyKNMe8ZY7b5Mt3rdKZB6wkxxmwxxqxyUaYqY8wOY8xWY0ypG3IZYxKNMc8ZY/b43lsLXJBpku81OvnRboz5igtyfdX3Pt9pjHnK9/53JpO1NmA+gBCgAigEwoFtQPEYrv8SYDawc9CyHwN3+27fDfzId7vYly8CKPDlDhmFTJnAbN/tOGCvb92O5QIMEOu7HQZsBC5y+rXyrevvgSeBVW74+fnWVQWknrLM6ffVY8AdvtvhQKLTmU7JFwIcAS5w+L2eDRwAonz3nwE+41SmUXvBnfgAFgBrBt2/B7hnjDPk8+HCLwcyfbczgfLTZQPWAAvGIN+LwMfckguIBt4H5judCcgB3gCu4M+F7/jrxOkL37FcQLyvxIxbMp0m49XAO07nYqDwq4FkBi4pu8qXzZFMgTakc/LFPanGt8xJGdbaOgDf53Tf8jHPaozJB2YxsEXtaC7f0MlWoAF4zVrreCbgfuAfAe+gZU5nArDAq8aYzcaYlS7IVQg0Ar/xDX89YoyJcTjTqZYDT/luO5bLWnsYuA84BNQBbdbaV53KFGiFb06zzK3zTsc0qzEmFvgf4CvW2vazPfQ0y0Y8l7W231o7k4Gt6nnGmGlOZjLGXAc0WGs3D/VbTrNstH5+F1trZwNLgS8YYy45y2PHIlcoA0OXD1prZwGdDAxLOJnpzyszJhxYBjz7UQ89zbKRfl8lATcwMDyTBcQYY25zKlOgFX4NkDvofg5Q61CWk+qNMZkAvs8NvuVjltUYE8ZA2T9hrX3eLbkArLWtwDpgicOZLgaWGWOqgKeBK4wxv3M4EwDW2lrf5wbgf4F5DueqAWp8f5UBPMfALwDHXyufpcD71tp6330nc10FHLDWNlpre4HngYVOZQq0wt8ETDDGFPh+yy8HXnI400vA7b7btzMwhn5y+XJjTIQxpgCYALw30is3xhjg18Bua+1P3ZDLGJNmjEn03Y5i4D/FHiczWWvvsdbmWGvzGXjfvGmtvc3JTADGmBhjTNzJ2wyM/+50Mpe19ghQbYyZ5Ft0JbDLyUynuJU/D+ecXL9TuQ4BFxljon3/F68EdjuWaTR3nDjxAVzLwEyUCuCbY7zupxgYp+tl4Df1CiCFgR2B+3yfkwc9/pu+nOXA0lHKtIiBPwm3A1t9H9c6mQsoAbb4Mu0Evu1b7uhrNWhdl/HnnbZO//wKGZi1sQ0oO/medkGumUCp72f4ApDkdCbfeqKBZiBh0DKnX6t7Gdig2Qk8zsAMHEcy6dQKIiJBItCGdERE5AxU+CIiQUKFLyISJFT4IiJBQoUvIhIkVPgiIkFChS8iEiT+Pypxs02vXZ7eAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "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/ddbaXVLMuB8qvVSasreiU/1WkOriC/VE6OBQBYI8fDWqC6WiPHw0qPf9ejRVr5NZ3z2qUs8RsALJfp2HCpxX10ZS3t8846/661jbRhg9Klu2Y3d+hPAXjoXW2PA3jR3e8B8OLwvRBCiAmyY0B395cA3HhX8yMAnh6+fhrAp/d3WkIIIfbKqM/Qz7j7AgAM/57evykJIYQYhVu+KGpm583sgpld6LX4c04hhBDjMWpAv2JmtwPA8O8i+6C7P+Hu59z9XG2KL2oJIYQYj1ED+nMAHhu+fgzA1/dnOkIIIUZlN7LFrwJ4EMBJM7sI4IsA/gDAM2b2KwDeBvCZ3ezM+kBzOa2R6s2mry2DOpf/LCItJ9oIEuNcnk6Pma2zrGFALUgK1SMyqCgx1gqRGm4EEsRumyTG6nJpIkj9VgAAS6gVJGRqFem+9Saf92qT1H1tBvK/Zvp/cscbR+mYuYLXhJ3KSXKuoBIpky1GcsslkqhpqcNrxUYSxPWN9L5KciwAADrpeVuXHwtGJKxRWcxeLZB8NtKSxtYUt93NqbQdrk9z280WaXlio8YTh0WJ5Ng5G/loZSXdZ8s8BhUr6TnU18YvNrpjQHf3R0nXT429dyGEEPuGfikqhBAVQQFdCCEqggK6EEJUBAV0IYSoCAdagi7rOZo30ivgvSZZac/4NadDEj+tbvCvtTaVXjXP6jzZUFSWzcmqeb/k8/Y2USJ0+JiM9BVB2aogpxgnuMQPiFkHDW7v9WZaObBB1C8AcKNJkiQ1eZKkZlASrZ7v3RCsBF0nSHLFyo512lzxMAiOVSNlGWstrtRgeamyMpCsjJCcy4PI0a+nD6J+MzhOptI22pjix0lOztkocRg7XwGgT5LZ9QMfZavpMY0lfiI1ltLxpLk0ygn7rvmMvQUhhBDvCxTQhRCiIiigCyFERVBAF0KIiqCALoQQFUEBXQghKsKByhatN0DjSrpwnudpqZr1+RRzIt8qV/mY/lRaZsQkeQAQlCGkkq8iUCBlnfS8g1KIyKgcjY8JcopRou86IOq7QT2QgjWYhI0nFes30jtaa3D56GoRfNl8hKRHJGEV+kGSKyItzdt8TBH01dKlUJHzPGT0OMnLwAYjyBaj84UdD/1GcJw007brkfMV4MdjGSQOC8rIIiPJ7IK8fSjW0mMay3xHU9fTBm/cCALALtEduhBCVAQFdCGEqAgK6EIIUREU0IUQoiIooAshREVQQBdCiIpwsLLFsof88vVk35SnZT55m9fzq6+lp1/O8OtUn5Q17BeBRivoYjKorMdlS3k33ZcF0jJSFjPcTyTRokNC2SKRowW2o/aOJGykb8DLtGJQ4/I2Uh40hEk+jSd1pBLSPJC91dqBz4k8sdbmEk16bEXHySjZFvPAfyPIFlm2VXb8RPuJ5haREd/mHW672ka6r7HCdcvNa2l5Yu3KTT65XaI7dCGEqAgK6EIIUREU0IUQoiIooAshREVQQBdCiIqwo8rFzJ4E8DCARXe/d9j2uwB+FcDV4ce+4O7P77QtL3voLVxO9uVlWiLQWJ+n2yuW0wqY/jSXQ7BkUYMiuLaNonLp85XxrJuWFWQlVy9YmV41t36QlCpKzsW+blDDldloUA8SKLH6ksQPQKCSCJKARUoWz4gaIhI2EduF6iWiUopUEnknUKy00j7POlxqk3WJuqIXHFsDJukJ7F0LauYWJAFeUHuWJWuLjhOWBGwwqsqFnLNZ4L9iI+2LfIUn2sqXVpPtg8Vrwex2x27u0J8C8FCi/Y/d/b7hvx2DuRBCiFvLjgHd3V8CcOMA5iKEEGIMxnmG/jkz+08ze9LMju3bjIQQQozEqAH9zwB8GMB9ABYA/CH7oJmdN7MLZnahRPCTOSGEEGMxUkB39yvu3nf3AYA/B3B/8Nkn3P2cu58rEPyOVwghxFiMFNDN7PZtb38RwKv7Mx0hhBCjYk6SYv3fB8y+CuBBACcBXAHwxeH7+7Ap2nsLwK+5+8KOOzO7CuD7w7cnAYyv0/nhR3bYRHaQDbaQHd5rgw+6+6mdBu0Y0G8VZnbB3c9NZOfvI2SHTWQH2WAL2WF0G+iXokIIUREU0IUQoiJMMqA/McF9v5+QHTaRHWSDLWSHEW0wsWfoQggh9hc9chFCiIowkYBuZg+Z2ffM7A0ze3wSc5gEwzQJi2b26ra242b2gpn9z/BvpdMomNldZvYtM3vNzL5rZp8fth82OzTN7F/N7D+Gdvi9YfuhsgMAmFluZv9uZn8/fH8YbfCWmX3HzL5tZheGbXu2w4EHdDPLAfwpgJ8F8FEAj5rZRw96HhPiKbw3c+XjAF5093sAvDh8X2V6AH7T3T8C4OMAfn3o/8Nmhw6AT7r7x7D5m46HzOzjOHx2AIDPA3ht2/vDaAMA+MQwe+2WXHHPdpjEHfr9AN5w9zfdvQvgawAemcA8DhySufIRAE8PXz8N4NMHOaeDxt0X3P3fhq9XsXkin8Xhs4O7+9rwbTH85zhkdjCzOwH8PIC/2NZ8qGwQsGc7TCKgnwXwzrb3F4dth5UzW7+yHf49PeH5HBhmdjeAHwfwLziEdhg+avg2gEUAL7j7YbTDnwD4bfxgOZbDZgNg82L+TTN7xczOD9v2bIcdKxbdAlLlRCS1OWSY2SyAvwHwG+6+YkF1nKri7n0A95nZUQDPmtm9E57SgWJmW5XQXjGzByc8nUnzgLtfMrPTAF4ws9dH2cgk7tAvArhr2/s7AVyawDzeL1zZSnY2/Ls44fnccsyswGYw/yt3/9th86Gzwxbuvgzgn7C5vnKY7PAAgF8ws7ew+ej1k2b2lzhcNgAAuPul4d9FAM9i89H0nu0wiYD+MoB7zOxDZlYH8FkAz01gHu8XngPw2PD1YwC+PsG53HJs81b8ywBec/c/2tZ12OxwanhnDjObAvDTAF7HIbKDu/+Ou9/p7ndjMw78o7v/Mg6RDQDAzGbMbG7rNYBPYTOD7Z7tMJEfFpnZz2Hz2VkO4El3/9KBT2ICkMyVfwfgGQAfAPA2gM+4e2VL/pnZTwL4ZwDfwf8/N/0CNp+jHyY7/Bg2F7pybN5YPePuv29mJ3CI7LDF8JHLb7n7w4fNBmb2I9i8Kwc2H4N/xd2/NIod9EtRIYSoCPqlqBBCVAQFdCGEqAgK6EIIUREU0IUQoiIooAshREVQQBdCiIqggC6EEBVBAV0IISrC/wL089BcaKPlCQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "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 +} diff --git a/Task 4/Assignment_4_eng.pdf b/Task 4/Assignment_4_eng.pdf new file mode 100644 index 0000000..bf72f79 Binary files /dev/null and b/Task 4/Assignment_4_eng.pdf differ diff --git a/Task 4/Assignment_4_ru.pdf b/Task 4/Assignment_4_ru.pdf new file mode 100644 index 0000000..1310f1f Binary files /dev/null and b/Task 4/Assignment_4_ru.pdf differ diff --git a/Task 4/Lecture 4.ipynb b/Task 4/Lecture 4.ipynb new file mode 100644 index 0000000..eeca9db --- /dev/null +++ b/Task 4/Lecture 4.ipynb @@ -0,0 +1,1641 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from IPython.display import display, HTML\n", + "display(HTML(\"\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Lecture 4: Conditioning. Floating point arithmetic and stability. Systems of linear equations. " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "## Introductory example\n", + "Consider solving a linear systems of equations $$Ax = f$$ with\n", + "\n", + "$$A_{ij} = \\frac{1}{i + j + 1}, \\quad i,j = 0, \\ldots, n-1.$$\n", + "\n", + "For example, at $n=5$:\n", + "\n", + "$$\n", + "A=\\begin{bmatrix}1&{\\frac {1}{2}}&{\\frac {1}{3}}&{\\frac {1}{4}}&{\\frac {1}{5}}\\\\{\\frac {1}{2}}&{\\frac {1}{3}}&{\\frac {1}{4}}&{\\frac {1}{5}}&{\\frac {1}{6}}\\\\{\\frac {1}{3}}&{\\frac {1}{4}}&{\\frac {1}{5}}&{\\frac {1}{6}}&{\\frac {1}{7}}\\\\{\\frac {1}{4}}&{\\frac {1}{5}}&{\\frac {1}{6}}&{\\frac {1}{7}}&{\\frac {1}{8}}\\\\{\\frac {1}{5}}&{\\frac {1}{6}}&{\\frac {1}{7}}&{\\frac {1}{8}}&{\\frac {1}{9}}\\end{bmatrix}.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def hilbert_fl(n):\n", + " return np.array([[1.0/(i + j + 1) for i in range(n)] for j in range(n)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Lets try to solve a small system numerically" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "n = 3\n", + "A = hilbert_fl(n)\n", + "f = np.random.randn(n)\n", + "x = np.linalg.solve(A, f) # we have not discussed yet what this function is doing but apparently it solves a system of equations" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.7277615519878704e-15\n" + ] + } + ], + "source": [ + "err = np.linalg.norm(A @ x - f) / np.linalg.norm(f)\n", + "print(err)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### So far so good, what about a larger system?" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "n = 15\n", + "A = hilbert_fl(n)\n", + "f = np.random.randint(-3,3,n)\n", + "x = np.linalg.solve(A, f)" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.5838742081211422\n" + ] + } + ], + "source": [ + "err = np.linalg.norm(A @ x - f) / np.linalg.norm(f)\n", + "print(err)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "The error grows with increase of $n$, and we have to find out why. It is not a problem with the algorithm!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Lets solve the system exactly symbolically " + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import Rational\n", + "from sympy.matrices import Matrix\n", + "def hilbert_sym(n):\n", + " return Matrix([[Rational(1, i+j+1) for i in range(n)] for j in range(n)])" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 0.5 ],\n", + " [0.5 , 0.33333333]])" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hilbert_fl(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}1 & \\frac{1}{2}\\\\\\frac{1}{2} & \\frac{1}{3}\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 1, 1/2],\n", + "[1/2, 1/3]])" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hilbert_sym(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "n = 12" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": {}, + "outputs": [], + "source": [ + "A_sym = hilbert_sym(n)\n", + "A_sym_inv = A_sym.inv(method=\"LU\")" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{array}{cccccccccccc}144 & -10296 & 240240 & -2702700 & 17297280 & -68612544 & 176432256 & -299304720 & 332560800 & -232792560 & 93117024 & -16224936\\\\-10296 & 981552 & -25765740 & 309188880 & -2061259200 & 8409937536 & -22076086032 & 38044955520 & -42800574960 & 30263032800 & -12206089896 & 2141691552\\\\240240 & -25765740 & 721440720 & -9018009000 & 61837776000 & -257554337040 & 686811565440 & -1198416098880 & 1361836476000 & -970938969000 & 394350596640 & -69604975440\\\\-2702700 & 309188880 & -9018009000 & 115945830000 & -811620810000 & 3434057827200 & -9271956133440 & 16342037712000 & -18725251545000 & 13443770340000 & -5492740453200 & 974469656160\\\\17297280 & -2061259200 & 61837776000 & -811620810000 & 5771525760000 & -24725216355840 & 67432408243200 & -119841609888000 & 138278780640000 & -99868008240000 & 41012462050560 & -7308522421200\\\\-68612544 & 8409937536 & -257554337040 & 3434057827200 & -24725216355840 & 106992754412544 & -294230074634496 & 526565596677120 & -611192210428800 & 443680271274240 & -183018111900624 & 32742180446976\\\\176432256 & -22076086032 & 686811565440 & -9271956133440 & 67432408243200 & -294230074634496 & 814790975910912 & -1466861305029120 & 1711338189200640 & -1247850762958800 & 516757021837056 & -92769511266432\\\\-299304720 & 38044955520 & -1198416098880 & 16342037712000 & -119841609888000 & 526565596677120 & -1466861305029120 & 2654320456719360 & -3110531785218000 & 2276990587872000 & -946216088737920 & 170392979877120\\\\332560800 & -42800574960 & 1361836476000 & -18725251545000 & 138278780640000 & -611192210428800 & 1711338189200640 & -3110531785218000 & 3659449159080000 & -2688113888460000 & 1120519052452800 & -202341663604080\\\\-232792560 & 30263032800 & -970938969000 & 13443770340000 & -99868008240000 & 443680271274240 & -1247850762958800 & 2276990587872000 & -2688113888460000 & 1980715496760000 & -827939077645680 & 149882713780800\\\\93117024 & -12206089896 & 394350596640 & -5492740453200 & 41012462050560 & -183018111900624 & 516757021837056 & -946216088737920 & 1120519052452800 & -827939077645680 & 346945899203904 & -62950739787936\\\\-16224936 & 2141691552 & -69604975440 & 974469656160 & -7308522421200 & 32742180446976 & -92769511266432 & 170392979877120 & -202341663604080 & 149882713780800 & -62950739787936 & 11445589052352\\end{array}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 144, -10296, 240240, -2702700, 17297280, -68612544, 176432256, -299304720, 332560800, -232792560, 93117024, -16224936],\n", + "[ -10296, 981552, -25765740, 309188880, -2061259200, 8409937536, -22076086032, 38044955520, -42800574960, 30263032800, -12206089896, 2141691552],\n", + "[ 240240, -25765740, 721440720, -9018009000, 61837776000, -257554337040, 686811565440, -1198416098880, 1361836476000, -970938969000, 394350596640, -69604975440],\n", + "[ -2702700, 309188880, -9018009000, 115945830000, -811620810000, 3434057827200, -9271956133440, 16342037712000, -18725251545000, 13443770340000, -5492740453200, 974469656160],\n", + "[ 17297280, -2061259200, 61837776000, -811620810000, 5771525760000, -24725216355840, 67432408243200, -119841609888000, 138278780640000, -99868008240000, 41012462050560, -7308522421200],\n", + "[ -68612544, 8409937536, -257554337040, 3434057827200, -24725216355840, 106992754412544, -294230074634496, 526565596677120, -611192210428800, 443680271274240, -183018111900624, 32742180446976],\n", + "[ 176432256, -22076086032, 686811565440, -9271956133440, 67432408243200, -294230074634496, 814790975910912, -1466861305029120, 1711338189200640, -1247850762958800, 516757021837056, -92769511266432],\n", + "[-299304720, 38044955520, -1198416098880, 16342037712000, -119841609888000, 526565596677120, -1466861305029120, 2654320456719360, -3110531785218000, 2276990587872000, -946216088737920, 170392979877120],\n", + "[ 332560800, -42800574960, 1361836476000, -18725251545000, 138278780640000, -611192210428800, 1711338189200640, -3110531785218000, 3659449159080000, -2688113888460000, 1120519052452800, -202341663604080],\n", + "[-232792560, 30263032800, -970938969000, 13443770340000, -99868008240000, 443680271274240, -1247850762958800, 2276990587872000, -2688113888460000, 1980715496760000, -827939077645680, 149882713780800],\n", + "[ 93117024, -12206089896, 394350596640, -5492740453200, 41012462050560, -183018111900624, 516757021837056, -946216088737920, 1120519052452800, -827939077645680, 346945899203904, -62950739787936],\n", + "[ -16224936, 2141691552, -69604975440, 974469656160, -7308522421200, 32742180446976, -92769511266432, 170392979877120, -202341663604080, 149882713780800, -62950739787936, 11445589052352]])" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A_sym_inv" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}225 & -25200 & 928200\\\\-25200 & 3763200 & -155937600\\\\928200 & -155937600 & 6892441920\\end{matrix}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 225, -25200, 928200],\n", + "[-25200, 3763200, -155937600],\n", + "[928200, -155937600, 6892441920]])" + ] + }, + "execution_count": 51, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A_sym_inv[:3,:3]" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": {}, + "outputs": [], + "source": [ + "x_sym = A_sym_inv @ f" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([178587533700, -36343113031680, 1840126802426640], dtype=object)" + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x_sym[:3]" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([2, -2, -1, 2, 2, -2, 1, -2, -1, -2, 0, -2, -2, -1, -2],\n", + " dtype=object)" + ] + }, + "execution_count": 54, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A_sym @ x_sym" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 2, -2, -1, 2, 2, -2, 1, -2, -1, -2, 0, -2, -2, -1, -2])" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Thats cool: we found an exact solution. How does it relate to the one we found numerically above?" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 1.68481673e+09, -2.49810169e+11, 9.23256364e+12, -1.47865831e+14,\n", + " 1.26548177e+15, -6.35346650e+15, 1.90893727e+16, -3.10903988e+16,\n", + " 8.47002490e+15, 7.74591773e+16, -1.83635048e+17, 2.14030465e+17,\n", + " -1.44209897e+17, 5.37157771e+16, -8.60260839e+15])" + ] + }, + "execution_count": 56, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([178587533700, -36343113031680, 1840126802426640,\n", + " -40763720391095520, 494478649301036400, -3692652241725326880,\n", + " 18178212566190156480, -61465382245798399440, 146065938596015157720,\n", + " -246060765797279312880, 292261743908214880800,\n", + " -239289867888435655200, 128509153948247763600,\n", + " -40739658302071890000, 5777758685388081600], dtype=object)" + ] + }, + "execution_count": 57, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x_sym" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1494.0414789748738" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.linalg.norm(x - x_sym.astype(np.float64))/np.linalg.norm(x)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So we have two errors: error on the solution, $|x-x_0|/|x_0|$ and error on the residual, $|A x-f|/|f|$. Note that we can easily estimate the second type of error but estimating the first one is tricky: you may think you have to know the true solution $x_0$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Lets try to further increase the system size $n$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Whats going on? Why dows np.linalg.solve fail so miserably?" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{array}{ccccccccccccccc}1 & \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15}\\\\\\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16}\\\\\\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17}\\\\\\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18}\\\\\\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19}\\\\\\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20}\\\\\\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21}\\\\\\frac{1}{8} & \\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22}\\\\\\frac{1}{9} & \\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23}\\\\\\frac{1}{10} & \\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23} & \\frac{1}{24}\\\\\\frac{1}{11} & \\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23} & \\frac{1}{24} & \\frac{1}{25}\\\\\\frac{1}{12} & \\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23} & \\frac{1}{24} & \\frac{1}{25} & \\frac{1}{26}\\\\\\frac{1}{13} & \\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23} & \\frac{1}{24} & \\frac{1}{25} & \\frac{1}{26} & \\frac{1}{27}\\\\\\frac{1}{14} & \\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23} & \\frac{1}{24} & \\frac{1}{25} & \\frac{1}{26} & \\frac{1}{27} & \\frac{1}{28}\\\\\\frac{1}{15} & \\frac{1}{16} & \\frac{1}{17} & \\frac{1}{18} & \\frac{1}{19} & \\frac{1}{20} & \\frac{1}{21} & \\frac{1}{22} & \\frac{1}{23} & \\frac{1}{24} & \\frac{1}{25} & \\frac{1}{26} & \\frac{1}{27} & \\frac{1}{28} & \\frac{1}{29}\\end{array}\\right]$" + ], + "text/plain": [ + "Matrix([\n", + "[ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15],\n", + "[ 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16],\n", + "[ 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17],\n", + "[ 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18],\n", + "[ 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19],\n", + "[ 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20],\n", + "[ 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21],\n", + "[ 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22],\n", + "[ 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23],\n", + "[1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23, 1/24],\n", + "[1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23, 1/24, 1/25],\n", + "[1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23, 1/24, 1/25, 1/26],\n", + "[1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23, 1/24, 1/25, 1/26, 1/27],\n", + "[1/14, 1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23, 1/24, 1/25, 1/26, 1/27, 1/28],\n", + "[1/15, 1/16, 1/17, 1/18, 1/19, 1/20, 1/21, 1/22, 1/23, 1/24, 1/25, 1/26, 1/27, 1/28, 1/29]])" + ] + }, + "execution_count": 60, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A_sym" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 0.5 , 0.33333333, 0.25 , 0.2 ,\n", + " 0.16666667, 0.14285714, 0.125 , 0.11111111, 0.1 ,\n", + " 0.09090909, 0.08333333, 0.07692308, 0.07142857, 0.06666667],\n", + " [0.5 , 0.33333333, 0.25 , 0.2 , 0.16666667,\n", + " 0.14285714, 0.125 , 0.11111111, 0.1 , 0.09090909,\n", + " 0.08333333, 0.07692308, 0.07142857, 0.06666667, 0.0625 ],\n", + " [0.33333333, 0.25 , 0.2 , 0.16666667, 0.14285714,\n", + " 0.125 , 0.11111111, 0.1 , 0.09090909, 0.08333333,\n", + " 0.07692308, 0.07142857, 0.06666667, 0.0625 , 0.05882353],\n", + " [0.25 , 0.2 , 0.16666667, 0.14285714, 0.125 ,\n", + " 0.11111111, 0.1 , 0.09090909, 0.08333333, 0.07692308,\n", + " 0.07142857, 0.06666667, 0.0625 , 0.05882353, 0.05555556],\n", + " [0.2 , 0.16666667, 0.14285714, 0.125 , 0.11111111,\n", + " 0.1 , 0.09090909, 0.08333333, 0.07692308, 0.07142857,\n", + " 0.06666667, 0.0625 , 0.05882353, 0.05555556, 0.05263158],\n", + " [0.16666667, 0.14285714, 0.125 , 0.11111111, 0.1 ,\n", + " 0.09090909, 0.08333333, 0.07692308, 0.07142857, 0.06666667,\n", + " 0.0625 , 0.05882353, 0.05555556, 0.05263158, 0.05 ],\n", + " [0.14285714, 0.125 , 0.11111111, 0.1 , 0.09090909,\n", + " 0.08333333, 0.07692308, 0.07142857, 0.06666667, 0.0625 ,\n", + " 0.05882353, 0.05555556, 0.05263158, 0.05 , 0.04761905],\n", + " [0.125 , 0.11111111, 0.1 , 0.09090909, 0.08333333,\n", + " 0.07692308, 0.07142857, 0.06666667, 0.0625 , 0.05882353,\n", + " 0.05555556, 0.05263158, 0.05 , 0.04761905, 0.04545455],\n", + " [0.11111111, 0.1 , 0.09090909, 0.08333333, 0.07692308,\n", + " 0.07142857, 0.06666667, 0.0625 , 0.05882353, 0.05555556,\n", + " 0.05263158, 0.05 , 0.04761905, 0.04545455, 0.04347826],\n", + " [0.1 , 0.09090909, 0.08333333, 0.07692308, 0.07142857,\n", + " 0.06666667, 0.0625 , 0.05882353, 0.05555556, 0.05263158,\n", + " 0.05 , 0.04761905, 0.04545455, 0.04347826, 0.04166667],\n", + " [0.09090909, 0.08333333, 0.07692308, 0.07142857, 0.06666667,\n", + " 0.0625 , 0.05882353, 0.05555556, 0.05263158, 0.05 ,\n", + " 0.04761905, 0.04545455, 0.04347826, 0.04166667, 0.04 ],\n", + " [0.08333333, 0.07692308, 0.07142857, 0.06666667, 0.0625 ,\n", + " 0.05882353, 0.05555556, 0.05263158, 0.05 , 0.04761905,\n", + " 0.04545455, 0.04347826, 0.04166667, 0.04 , 0.03846154],\n", + " [0.07692308, 0.07142857, 0.06666667, 0.0625 , 0.05882353,\n", + " 0.05555556, 0.05263158, 0.05 , 0.04761905, 0.04545455,\n", + " 0.04347826, 0.04166667, 0.04 , 0.03846154, 0.03703704],\n", + " [0.07142857, 0.06666667, 0.0625 , 0.05882353, 0.05555556,\n", + " 0.05263158, 0.05 , 0.04761905, 0.04545455, 0.04347826,\n", + " 0.04166667, 0.04 , 0.03846154, 0.03703704, 0.03571429],\n", + " [0.06666667, 0.0625 , 0.05882353, 0.05555556, 0.05263158,\n", + " 0.05 , 0.04761905, 0.04545455, 0.04347826, 0.04166667,\n", + " 0.04 , 0.03846154, 0.03703704, 0.03571429, 0.03448276]])" + ] + }, + "execution_count": 62, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Condition number of a problem" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Conditioning of a problem" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- Conditioning of a problem is defined without reference to any particular algorithm." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- A well-conditioned problem $f(x)$: all small perturbations of $x$ lead to small changes in $f(x)$." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- An ill-conditioned problem $f(x)$: some small perturbation of $x$ leads to large changes in $f(x)$." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "The absolute condition number of a problem $f(x)$ is:\n", + "\n", + "$$\\kappa = \\lim _{\\varepsilon \\rightarrow 0}\\,\\sup _{\\|\\delta x\\|\\,\\leq \\,\\varepsilon }{\\frac {\\|\\delta f(x)\\|}{\\|\\delta x\\|}}$$\n", + "\n", + "and the relative condition number is:\n", + "\n", + "$$\\kappa = \\lim _{\\varepsilon \\rightarrow 0}\\,\\sup _{\\|\\delta x\\|\\,\\leq \\,\\varepsilon }{\\frac {\\|\\delta f(x)\\|/\\|f(x)\\|}{\\|\\delta x\\|/\\|x\\|}}.$$\n", + "\n", + "Relative $\\kappa$ is normally more important (its at least dimensionless)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "If the function $f(x)$ is differentiable:\n", + "$$\n", + "\\kappa = \\frac{\\Vert J(x) \\Vert}{\\Vert f(x)\\Vert / \\Vert x \\Vert},\n", + "$$\n", + "where $J$ is Jacobian matrix: $J _{ij}={\\frac {\\partial f_{i}}{\\partial x_{j}}}$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Well- and ill-conditioned problems" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- A problem is well-conditioned if $\\kappa$ is small ($\\lesssim 10^2$) and ill-conditioned if $\\kappa$ is large ($\\gtrsim 10^{6}$)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Consider a problem of computing $f(x)=\\sqrt{x}$:\n", + "$$\\kappa = \\frac{\\Vert J(x) \\Vert}{\\Vert f(x)\\Vert / \\Vert x \\Vert} = \\frac{1/(2\\sqrt{x})}{\\sqrt{x}/x}=\\frac12.$$\n", + "This problem is well-conditioned by any standard." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Consider computing $f(x)=\\tan x$ for $x=10^{100}$. \n", + "The problem is ill-conditioned. Why?\n", + "\n", + "[Surely You're Joking, Mr. Feynman](https://edisciplinas.usp.br/pluginfile.php/4420924/mod_resource/content/1/LuckyNumbers_Feynman.pdf)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Condition numbers in linear algebra computations" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "- In matrix computations, the basic operations are computing $f(x) = A x$ and computing $g(x) = A^{-1}x$ (that is, solving a linear system)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- They have the condition numbers\n", + "$\\kappa_f(x) = \\Vert A \\Vert \\frac{\\Vert x \\Vert}{\\Vert A x \\Vert}$ and $\\kappa_g(x)=\\Vert A^{-1} \\Vert \\frac{\\Vert A x \\Vert}{\\Vert x \\Vert}$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- The worst-case bounds (independent on $x$):\n", + "$$\n", + "\\kappa_{f,g}(x) \\leq \\kappa(A),\n", + "$$\n", + "where the right-hand-side is known as a condition number of a matrix $A$: \n", + "$$\\kappa(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert.$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- Note, that the condition number is different for different norms. Which are worst-case vectors $x$ for $f(x)$ and $g(x)$ and $L_2$ norm?" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Condition number of a linear system\n", + "Now consider the perturbation of a system of linear equations:\n", + "from\n", + "$$ A x = f$$\n", + "to\n", + "$$ (A + \\Delta A) (x + \\Delta x) = f + \\Delta f.$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "Perturbation to the solution:\n", + "\n", + "$$\n", + "\\frac{\\Vert \\Delta x \\Vert}{\\Vert x \\Vert} \\leq \\frac{\\mathrm{cond}(A)}{1 - \\mathrm{cond}(A)\\frac{\\|\\Delta A\\|}{\\|A\\|}} \\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big)\n", + "$$\n", + "\n", + "The crucial role is played by the condition number $\\kappa(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert$." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "The larger the condition number, the less number of digits we can recover." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Note, that if $\\Delta A = 0$, then\n", + "\n", + "$$ \\frac{\\Vert \\Delta x \\Vert}{\\Vert x \\Vert} \\leq \\mathrm{cond}(A) \\frac{\\|\\Delta f\\|}{\\|f\\|} $$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "The condition number in spectral norm is equal to the ratio of the largest singular value and the smallest singular value.\n", + "\n", + "$$ \\mathrm{cond}_2 (A) = \\|A\\|_2 \\|A^{-1}\\|_2 = \\frac{\\sigma_{\\max}}{\\sigma_{\\min}} $$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Hilbert matrix (again)" + ] + }, + { + "cell_type": "code", + "execution_count": 87, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2.495951750009794e+17" + ] + }, + "execution_count": 87, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n = 15\n", + "A = hilbert_fl(n)\n", + "np.linalg.cond(A)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "The Hilbert matrix is very poorly conditioned, indeed. This implies that solution of the equation $Ax = b$ for x is very sensitive to the small changes in $A$ and $b$. Here it is time to recall how exactly are real numbers represented in a CPU." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Floating point\n", + "The numbers in computer memory are typically represented as floating point numbers \n", + "\n", + "A floating point number is represented as \n", + "\n", + "$$\\textrm{number} = \\textrm{significand} \\times \\textrm{base}^{\\textrm{exponent}},$$\n", + "\n", + "where *significand* is integer, *base* is positive integer and *exponent* is integer (can be negative), i.e.\n", + "\n", + "$$ 1.2 = 12 \\cdot 10^{-1}.$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## IEEE 754\n", + "In modern computers, the floating point representation is controlled by [IEEE 754 standard](https://en.wikipedia.org/wiki/IEEE_floating_point) which was published in 1985. Before that point different computers behaved differently with floating point numbers. \n", + "\n", + "IEEE 754 has:\n", + "- Floating point representation (as described above), $(-1)^s \\times c \\times b^q$.\n", + "- Two infinities, $+\\infty$ and $-\\infty$\n", + "- Rules for assigning kinds of NaN\n", + "- Rules for rounding\n", + "- Rules for $\\frac{0}{0}, \\frac{1}{-0}, \\ldots$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## The two most common format, single & double\n", + "\n", + "The two most common formats, called binary32 and binary64 (called also single and double formats). Recently, the format binary16 plays important role in learning deep neural networks.\n", + "\n", + "| Name | Common Name | Base | Digits | Emin | Emax |\n", + "|------|----------|----------|-------|------|------|\n", + "|binary16| half precision | 2 | 11 | -14 | + 15 |\n", + "|binary32| single precision | 2 | 24 | -126 | + 127 | \n", + "|binary64| double precision | 2 | 53 | -1022 | +1023 | \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "For all real $x$ (in an admissible range, such as $[1.1\\cdot 10^{-38}, 3.4\\cdot 10^{38}]$ for single precision) there exists floating point number $\\tilde x$ such that $|x-\\tilde x| \\leq |x|\\epsilon_{machine}$\n", + "- For single precision: $$\\epsilon_{machine} = 2^{-24}\\approx 5.96\\cdot 10^{-8}$$\n", + "- For double precision: $$\\epsilon_{machine} = 2^{-53}\\approx 1.11\\cdot 10^{-16}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For single precision, we have $8388608$ numbers in $(1,2)$ which is quite dense. But we still have the same $8388608$ numbers in $(2^{126}, 2^{127})$ -- much larger interval! At each arithmetic operation, the result will typically not fit into the grid and will be replaced by some proxy from the grid. This can be disastrous for the final result if accumulates as a result of careless programming." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Accuracy and memory\n", + "\n", + "The relative accuracy of single precision is $10^{-7}-10^{-8}$, while for double precision is $10^{-14}-10^{-16}$.\n", + "\n", + "- A float16 takes 2 bytes, float32 takes 4 bytes, float64, or double precision, takes 8 bytes\n", + "\n", + "- These are the only two floating point-types supported in hardware (float32 and float64).\n", + "\n", + "- Normally, double precision is for Scientific Computing and float is for GPU/Data Science." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Accuracy demos" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.5 0.6 0.7 0.8]\n", + "0.7999999999999999\n" + ] + } + ], + "source": [ + "z = np.arange(0.5, 0.9, 0.1)\n", + "print(z)\n", + "print(z[-1])" + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4.3670121385730454e-07\n" + ] + } + ], + "source": [ + "a = np.float32(6.0)\n", + "b = np.sqrt(a)\n", + "print(b ** 2 - a)" + ] + }, + { + "cell_type": "code", + "execution_count": 102, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n" + ] + } + ], + "source": [ + "a = np.array(0.4147369269524216, dtype=np.float32)\n", + "b = np.exp(a)\n", + "print(np.log(b) - a)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Loss of significance\n", + "\n", + "- Many operations lead to the [loss of significance](https://en.wikipedia.org/wiki/Loss_of_significance)\n", + "- For example, it is a bad idea to subtract two big numbers that are close, the difference will have fewer correct digits\n", + "- This is related to algorithms and their properties (stability)" + ] + }, + { + "cell_type": "code", + "execution_count": 103, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "0.0\n", + "2.0\n" + ] + } + ], + "source": [ + "import math\n", + "test_list = [1, 1e20, 1, -1e20]\n", + "print(np.sum(test_list))\n", + "print(1 + 1e20 + 1 - 1e20)\n", + "print(math.fsum(test_list))" + ] + }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.0\n" + ] + } + ], + "source": [ + "print(1e20 - 1e20 + 1 + 1)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Some disasters attibuted to bad numerical computing\n", + "- [Vancouver Stock Exchange Index computation error](https://www5.in.tum.de/~huckle/Vancouv.pdf): in January 1982 the index was initialized at 1000 and subsequently updated and truncated to three decimal places on each trade. Such a thing was done about 3000 times each day. The accumulated truncations led to an erroneous loss of around 25 points per month. On 25 Nov 1983, the error was corrected, raising the value of the index from its Friday closing figure of 524 to 1098.\n", + "- [The Patriot Missile failure](https://www-users.cse.umn.edu/~arnold/disasters/patriot.html), happened In Dharan, Saudi Arabia, on February 25, 1991, resulted in 28 deaths, is ultimately attributable to poor handling of rounding errors.\n", + "- [The explosion of the Ariane 5 rocket](https://www-users.cse.umn.edu/~arnold/disasters/ariane.html) just after lift-off on its maiden voyage off French Guiana, on June 4, 1996, was ultimately the consequence of a simple overflow.\n", + "- [The sinking of the Sleipner A offshore platform](https://www-users.cse.umn.edu/~arnold/disasters/sleipner.html) in Gandsfjorden near Stavanger, Norway, on August 23, 1991, resulted in a loss of nearly one billion dollars. It was found to be the result of inaccurate finite element analysis. " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Stability of algorithms\n", + "\n", + "- Let $x$ be an object (for example, a vector) \n", + "- Let $f(x)$ be the function (functional) you want to evaluate \n", + "- You have a numerical algorithm working in finite precision arithmetics will deliver an approximation $alg(x)$ to the correct result $f(x)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "The algorithm is called forward stable, if $$\\Vert alg(x) - f(x) \\Vert \\leq \\varepsilon $$ \n", + "\n", + "The algorithm is called backward stable, if for any $x$ there is a close vector $x + \\delta x$ such that\n", + "\n", + "$$alg(x) = f(x + \\delta x)$$\n", + "\n", + "and $\\Vert \\delta x \\Vert$ is small.\n", + "\n", + "The algorithm is called stable if it gives approximately correct answers for approximately correct problems: for each $x$ there exists $\\tilde x$ such that $\\Vert x-\\tilde x\\Vert = \\mathcal{O}(\\epsilon_{machine})$ such that\n", + "$$\n", + "\\Vert alg(x) - f(\\tilde x) \\Vert = \\mathcal{O}(\\epsilon_{machine}) \n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The best algorithms usually compute exact answers for slightly perturbed data. If the output is not too sentitive to input data (characterized by the condition number), then this is sufficient in practice." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Examples of instability - 1\n", + "\n", + "How to compute the following function in numerically stable manner?\n", + "\n", + "$$\\log(1 - \\tanh^2(u))$$" + ] + }, + { + "cell_type": "code", + "execution_count": 107, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Original function: -inf\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\Zenon\\AppData\\Local\\Temp\\ipykernel_21092\\1730504921.py:2: RuntimeWarning: divide by zero encountered in log\n", + " print(\"Original function:\", np.log(1 - np.tanh(u)**2))\n" + ] + } + ], + "source": [ + "u = 30\n", + "print(\"Original function:\", np.log(1 - np.tanh(u)**2))" + ] + }, + { + "cell_type": "code", + "execution_count": 109, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Use more numerically stable form: -58.61370563888011\n" + ] + } + ], + "source": [ + "print(\"Use more numerically stable form:\", np.log(4) - 2 * np.log(np.exp(-u) + np.exp(u)))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Examples of instability - 2\n", + "How to compute the following function in numerically stable manner?\n", + "$$SoftMax(x)_j = \\dfrac{e^{x_j}}{\\sum\\limits_{i=1}^n e^{x_i}}$$" + ] + }, + { + "cell_type": "code", + "execution_count": 118, + "metadata": {}, + "outputs": [], + "source": [ + "n = 5\n", + "x = np.random.randn(n)\n", + "x[0] = 10" + ] + }, + { + "cell_type": "code", + "execution_count": 119, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[9.99869717e-01 4.25877087e-05 1.71219662e-05 4.14169156e-05\n", + " 2.91559416e-05]\n" + ] + } + ], + "source": [ + "print(np.exp(x) / np.sum(np.exp(x)))" + ] + }, + { + "cell_type": "code", + "execution_count": 120, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([10. , -0.06381458, -0.97501805, -0.09169088, -0.44272154])" + ] + }, + "execution_count": 120, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x" + ] + }, + { + "cell_type": "code", + "execution_count": 121, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[9.99869717e-01 4.25877087e-05 1.71219662e-05 4.14169156e-05\n", + " 2.91559416e-05]\n" + ] + } + ], + "source": [ + "print(np.exp(x - np.max(x)) / np.sum(np.exp(x - np.max(x))))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Linear systems " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "## Scales of linear systems\n", + "\n", + "In different applications, the typical size of the linear systems can be different. \n", + "\n", + "- Small: $n \\leq 10^4$ (full matrix can be stored in memory, dense matrix)\n", + "- Medium: $n = 10^4 - 10^6$ (typically, sparse or structured matrix)\n", + "- Large: $n = 10^8 - 10^9$ (typically sparse matrix + parallel computations)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will talk only about approaches to small systems today." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Main questions about linear systems\n", + "\n", + "1. What is the accuracy we get from the solution (due to rounding errors)?\n", + "2. How we compute the solution? (using the Cramer's rule with dets is good only for $2 \\times 2$ matrices)\n", + "3. What is the complexity of the solution of linear systems?" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Solving $Ax=b$ via QR factorization\n", + "- Factor $A$ into $QR$: $A = QR$\n", + "- $Ax=b$ is equivalent to $QRx=b$ that is $Rx = y$ where $y = Q^*b$\n", + "- Solve triangular system $Rx=y$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- The algorithm is stable. It means that the solution $\\tilde x$ computed by this algorithm satisfies:\n", + "$$\n", + "\\frac{\\Vert\\tilde x - x\\Vert}{\\Vert x\\Vert} = \\mathcal{O}(\\kappa(A)\\epsilon_{machine})\n", + "$$\n", + "- Solving linear system with this method takes $2n^3$ flops [dominated by complexity of QR]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Solving $Ax=b$ via Gaussian elimination\n", + "- Gaussian elimination proceeds via the computation of one of the most important matrix decompositions: LU-decomposition.\n", + "\n", + "- LU-decomposition of the square matrix $A$ is the representation\n", + "\n", + "$$A = LU,$$\n", + "\n", + "where \n", + "- $L$ is lower triangular (elements strictly above the diagonal are zero)\n", + "- $U$ is upper triangular matrix (elements strictly below the diagonal are zero)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$$\n", + "\\begin{bmatrix}a_{11}&a_{12}&a_{13}\\\\a_{21}&a_{22}&a_{23}\\\\a_{31}&a_{32}&a_{33}\\end{bmatrix}=\\begin{bmatrix}\\ell _{11}&0&0\\\\\\ell _{21}&\\ell _{22}&0\\\\\\ell _{31}&\\ell _{32}&\\ell _{33}\\end{bmatrix}\\begin{bmatrix}u_{11}&u_{12}&u_{13}\\\\0&u_{22}&u_{23}\\\\0&0&u_{33}\\end{bmatrix}\n", + "$$\n", + "\n", + "This factorization is non-unique, so it is typical to require that the matrix $L$ has ones on the diagonal." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "LU decomposition is useful to solve a linear system since\n", + "\n", + "$$ A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, $$\n", + "\n", + "and this reduces to the solution of two linear systems forward step\n", + "\n", + "$$ L y = f, $$\n", + "\n", + "and backward step\n", + "\n", + "$$ U x = y. $$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Complexity of the Gaussian elimination/LU decomposition\n", + "\n", + "- The complexity is dominated by LU decomposition, which takes $\\frac23 n^3$: several times faster than QR \n", + "- In many cases, when solution with multiple right hand sides is required, computing LU-decomposition once is a good idea\n", + "- Once the decomposition is found then solving linear systems with $L$ and $U$ costs only $\\mathcal{O}(n^2)$ operations." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Computing LU decomposition\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## When LU fails\n", + "\n", + "- What happens, if pivots in the Gaussian elimination are really small or zero?. \n", + "\n", + "- There is classical $2 \\times 2$ example of a matrix with a bad LU decomposition. \n", + "\n", + "- The matrix we look at is \n", + "\n", + "$$\n", + " A = \\begin{pmatrix}\n", + " \\varepsilon & 1 \\\\\n", + " 1 & 1 \n", + " \\end{pmatrix}\n", + "$$\n", + "\n", + "- If $\\varepsilon$ is sufficiently small, we might fail." + ] + }, + { + "cell_type": "code", + "execution_count": 170, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "L:\n", + " [[1.e+00 0.e+00]\n", + " [1.e+18 1.e+00]]\n", + "U:\n", + " [[ 1.e-18 1.e+00]\n", + " [ 0.e+00 -1.e+18]]\n", + "A:\n", + " [[1.e-18 1.e+00]\n", + " [1.e+00 1.e+00]]\n", + "L * U - A:\n", + " [[ 0.00000000e+00 0.00000000e+00]\n", + " [-1.11022302e-16 -1.00000000e+00]]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "eps = 1e-18\n", + "A = np.array([[eps, 1],[1.0, 1]])\n", + "A0 = A.copy()\n", + "n = A.shape[0]\n", + "L = np.zeros((n, n))\n", + "U = np.zeros((n, n))\n", + "for k in range(n): \n", + " L[k, k] = 1\n", + " for i in range(k+1, n):\n", + " L[i, k] = A[i, k] / A[k, k]\n", + " for j in range(k+1, n):\n", + " A[i, j] = A[i, j] - L[i, k] * A[k, j]\n", + " for j in range(k, n):\n", + " U[k, j] = A[k, j]\n", + "print('L:\\n', L)\n", + "print('U:\\n', U)\n", + "print('A:\\n', A0)\n", + "print('L * U - A:\\n', np.dot(L, U) - A0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## The concept of pivoting\n", + "\n", + "- We can do pivoting, i.e. permute rows and columns to maximize $A_{kk}$ that we divide over. \n", + "\n", + "- The simplest but effective strategy is the row pivoting: at each step, select the index that is maximal in modulus, and put it onto the diagonal. \n", + "\n", + "- It gives us the decomposition \n", + "\n", + "$$A = P L U,$$\n", + "\n", + "where $P$ is a permutation matrix.\n", + "\n", + "Pivoting ensures $| L_{ij}|<1,$ but the elements of $U$ can grow, up to $2^n$. Theoretically, gaussian elimination is unstable, though (in practice, this is almost never encountered)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Linear equations: Summary\n", + "- Linear systems can be solved by Gaussian elimination, complexity is $\\mathcal{O}(n^3)$.\n", + "- Linear systems can be solved by LU-decomposition, complexity is $\\mathcal{O}(n^3)$ for the decomposition, $\\mathcal{O}(n^2)$ for each solve\n", + "- Linear least squares can be solved by normal equations (bad)\n", + "- Linear least squares can be solved by QR-decomposition (good)\n", + "- Without structure, we can solve up to $10^4$ linear systems on a laptop (memory restrictions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will discuss iterative methods for linear systems in what follows." + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Python 3.10.8 64-bit", + "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": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e" + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Task 4/tex/Assignment_4_eng.tex b/Task 4/tex/Assignment_4_eng.tex new file mode 100644 index 0000000..2859126 --- /dev/null +++ b/Task 4/tex/Assignment_4_eng.tex @@ -0,0 +1,161 @@ +\documentclass{article} +\usepackage[utf8]{inputenc} +\usepackage{biblatex} +\addbibresource{library.bib} +\usepackage{listings} +\usepackage{amssymb} +\usepackage{comment} +\usepackage{graphicx,amsmath} +\newcommand{\norm}[1]{\left\lVert#1\right\rVert} +\usepackage{hyperref} +\hypersetup{ + colorlinks=true, + linkcolor=blue, + filecolor=magenta, + urlcolor=cyan, + pdftitle={Overleaf Example}, + pdfpagemode=FullScreen, + } + +\title{Numerical Methods: Lecture 4. Conditioning. Floating point arithmetic and stability. Systems of linear equations.} +\author{Konstantin Tikhonov} + +\begin{document} + +\maketitle + +\section{Suggested Reading} + +\begin{itemize} +\item Lectures 12-19, 20-23 of \cite{trefethen1997numerical} +\item Lectures 6-7 of \cite{tyrtyshnikov2012brief} +\end{itemize} + +\section{Exercises} + +Deadline: 18 Nov + +\begin{enumerate} + +\item (3) Propose a numerically stable way to compute the function $f(x,a)=\sqrt{x+a}-\sqrt{x}$ for positive $x,\;a$. + +\item (2) Consider numerical evaluation $\mathcal{C}=\tan(10^{100})$ with the help of arbitrary-precision arithmetic module \lstinline{mpmath}, which can be called as follows: +\lstset{language=Python} +\lstset{frame=lines} +% \lstset{label={lst:code_direct}} +\lstset{basicstyle=\ttfamily} +\begin{lstlisting} +from mpmath import * +mp.dps = 64 # precision (in decimal places) +mp.pretty = True ++pi +\end{lstlisting} +What is the relative condition number of evaluating $\mathcal{C}$ w.r.t the input number $10^{100}$? How many digits do you need to keep at intermediate steps to evaluate $\mathcal{C}$ with 7-digit accuracy? +\begin{comment} +\item (3) Check, that the following function +\lstset{language=Python} +\lstset{frame=lines} +\lstset{label={lst:code_direct}} +\lstset{basicstyle=\ttfamily} +\begin{lstlisting} +import math +def round_to_n(x, n): + if x == 0: + return x + else: + return round(x, -int(math.floor(math.log10(abs(x)))) + (n - 1)) +\end{lstlisting} +rounds $x$ to $n$ significant digits. +A sample program to compute $\sum_{k=1}^{3000}k^{-2}\approx 1.6446$ via consequent summation with rounding of intermediate results to 4 digits looks as follows: +\lstset{language=Python} +\lstset{frame=lines} +\lstset{label={lst:code_direct}} +\lstset{basicstyle=\ttfamily} +\begin{lstlisting} +res = 0 +for k in range(1,3001): + res = round_to_n(res+1/k**2, 4) +\end{lstlisting} +Despite the absence of subtractions (and related precision loss), this code allows to get only two significant digits. Explain, why this happens and propose a more accurate way to compute this sum (maintaining the restriction of keeping only 4 digits of intermediate result). +\end{comment} +\item (4) Implement the function \lstinline{solve_quad(b, c)}, receiving coefficients $b$ and $c$ of a quadratic polynomial $x^2 + b x + c$, and returning a pair of equation roots. Your function should always return two roots, even for a degenerate case (for example, a call \lstinline{solve_quad(-2, 1)} should result into \lstinline{(1, 1)}). Additionally, your function is expected to return complex roots. + +After checking ensuring that your algorithm sort of works, try it on the following 5 tests. Make sure that all of them pass. +\lstset{language=Python} +\lstset{frame=lines} +\lstset{label={lst:code_direct}} +\lstset{basicstyle=\ttfamily} +\begin{lstlisting} +tests = [{'b': 4.0, 'c': 3.0}, + {'b': 2.0, 'c': 1.0}, + {'b': 0.5, 'c': 4.0}, + {'b': 1e10, 'c': 3.0}, + {'b': -1e10, 'c': 4.0}] +\end{lstlisting} + +\item (5) Consider the polynomial $$ +w(x)=\Pi_{r=1}^{20}(x-r)=\sum_{i=0}^{20} a_i x^i +$$ and investigate the condition number of roots of this polynomial w.r.t the coefficients $a_i$. Perform the following experiment, using \texttt{numpy} root-finding algorithm. Randomly perturb $w(x)$ by replacing the coefficients $a_i\to n_i a_i$, where $n_i$ is drawn from a normal distribution of mean $1$ and variance $\exp(-10)$. Show the results of $100$ such experiments in a single plot, along with the +roots of the unperturbed polynomial $w(x)$. Using one of the experiments, estimate the relative and absolute condition number of the problem of finding the roots of $w(x)$ w.r.t. polynomial coefficients. + +\item (10) +Consider the least squares problem $Ax\approx b$ at +$$ +A = \begin{bmatrix} +1 & 1\\ +1 & 1.00001\\ +1 & 1.00001 +\end{bmatrix},\quad b = \begin{bmatrix} +2 \\ +0.00001 \\ +4.00001 +\end{bmatrix}. +$$ + +\begin{itemize} + \item +Formally, solution is given by +\begin{equation} +\label{ex} +x = ( A^T A )^{-1} A^T b. +\end{equation} + Using this equation, compute the solution analytically. +\item Implement Eq. (\ref{ex}) in \lstinline{numpy} in single and double precision; compare the results to the analytical one. +\item Instead of Eq. (\ref{ex}), implement SVD-based solution to least squares. Which approach is numerically more stable? +\item Use \lstinline{np.linalg.lstsq} to solve the same equation. Which method does this function use? +\item +What are the four condition numbers of this problem, mentioned in Theorem 18.1 of Ref. \cite{trefethen1997numerical}? Give examples of perturbations $\delta b$ and $\delta A$ that approximately attain those condition numbers? +\end{itemize} +\item (7) +Let $$A = \begin{bmatrix} + \epsilon & 1 & 0\\ +1 & 1 & 1\\ +0 & 1 & 1 +\end{bmatrix}$$ +\begin{itemize} + \item Find analytically LU decomposition with and without pivoting for the matrix $A$. + \item Explain, why can the LU decomposition fail to approximate factors $L$ and $U$ for $|\epsilon|\ll 1$ in finite-precision arithmetic? +\end{itemize} + +\item (6) Consider computing the function $f(n, \alpha)$ defined by $f(0,\alpha)=\ln(1+1/\alpha)$ and recurrent relation +\begin{equation} + f(n,\alpha)=\frac{1}{n}-\alpha f(n-1,\alpha). +\end{equation} +Compute $f(20, 0.1)$ and $f(20, 10)$ in standard (double) precision. Now, do the same exercise in arbitrary +precision arithmetic: +\lstset{language=Python} +\lstset{frame=lines} +\lstset{label={lst:code_direct}} +\lstset{basicstyle=\ttfamily} +\begin{lstlisting} +from mpmath import mp, mpf +mp.dps = 64 # precision (in decimal places) +f = mp.zeros(1, n) +f[0] = mp.log(1+1/mpf(alpha)) +for i in range(1, n): + f[i] = 1/mpf(i) - mpf(alpha)*f[i-1] +\end{lstlisting} +\end{enumerate} +Plot the relative difference between exact and approximate results, in units of machine epsilon \texttt{np.finfo(float).eps} for $\alpha=0.1$ and $\alpha=10$ as function of $n$. How would you evaluate $f(30, 10)$ without relying on the arbitrary precision arithmetic? +\printbibliography +\end{document} diff --git a/Task 4/tex/Assignment_4_ru.tex b/Task 4/tex/Assignment_4_ru.tex new file mode 100644 index 0000000..4d1cdd3 --- /dev/null +++ b/Task 4/tex/Assignment_4_ru.tex @@ -0,0 +1,162 @@ +\documentclass[prb, notitlepage, aps, 11pt]{revtex4-2}% +\usepackage[utf8]{inputenc} +\usepackage[T2A]{fontenc} +\usepackage[english, russian]{babel} +\usepackage{amsmath} +\usepackage{enumitem} +\usepackage{amsmath} +\usepackage{delimset} +\usepackage[pdftitle = a]{hyperref} +\usepackage{datetime} +\usepackage{minted} +\usemintedstyle{friendly} +\usepackage[a]{esvect} +\hypersetup{ + colorlinks=true, + linkcolor=blue, + filecolor=magenta, + urlcolor=cyan, + pdfpagemode=FullScreen, +} +\usepackage{microtype} + +\newcommand{\framesep}{0.6em} +\BeforeBeginEnvironment{minted}{\vspace{-1.6em}} +\AfterEndEnvironment{minted}{\vspace{-0.5em}} + +\begin{document} + +\begin{center} + версия от \today \quad \currenttime +\end{center} + +\title{\texorpdfstring{ + Численные методы, осень 2022\\ + Задание 4 [Число обусловленности. Числа с плавающей точкой и вычислительная устойчивость] \\ + Всего баллов: 37 \ Срок сдачи: 18 ноября + }{} +} + +\maketitle + +\section*{Рекомендованная литература} + +\begin{itemize} +\item Лекции 12--19, 20--23 из \cite{trefethen1997numerical} +\item Лекции 6--7 из \cite{tyrtyshnikov2012brief} +\end{itemize} + +\section*{Упражнения} + +\begin{enumerate} + +\item (3) Предложите вычислительно устойчивый способ вычислить функцию +$f(x,a)=\sqrt{x+a}-\sqrt{x}$ +при положительных $x$ и $a$. + +\item (2) Вычислите $\mathcal{C}=\tan(10^{100})$ с помощью модуля \mintinline{python}{mpmath}, предназначенного для арифметики произвольной точности. Пример использования: +% +\begin{minted}[frame = lines, framesep = \framesep]{python} + from mpmath import * + mp.dps = 64 # точность (число десятичных цифр) + mp.pretty = True + +pi # pi — переменная из mpmath +\end{minted} +% +Чему равно относительное число обусловленности при вычислении $\mathcal{C} = \mathcal{C}(10^{100})$? Сколько цифр нужно хранить в памяти при промежуточных вычислениях, чтобы получить $\mathcal{C}$ с~точностью в 7 значащих цифр? + +\item (4) Реализуйте функцию \mintinline{python}{solve_quad(b, c)}, возвращающую корни приведённого квадратного уравнения $x^2 + b x + c = 0$. Корни могут повторяться или быть комплексными. + +Когда вам покажется, что функция работает, запустите её на следующих пяти тестах. Добейтесь того, чтобы она правильно работала на каждом из них. +% +\begin{minted}[frame = lines, framesep = \framesep]{python} + tests = [{'b': 4.0, 'c': 3.0}, + {'b': 2.0, 'c': 1.0}, + {'b': 0.5, 'c': 4.0}, + {'b': 1e10, 'c': 3.0}, + {'b': -1e10, 'c': 4.0}] +\end{minted} + +\item (5) Рассмотрите многочлен +$$ +w(x) = \prod_{r=1}^{20} (x-r) = \sum_{i=0}^{20} a_i x^i +$$ +и исследуйте число обусловленности его корней, выступающих в роли функций от коэффициентов $a_i$. Проведите эксперимент: случайным образом измените коэффициенты и найдите новые корни с помощью алгоритма из \texttt{numpy}. +Коэффициенты изменяйте по правилу $a_i \to n_i a_i$, где $n_i$ подчиняются нормальному распределению с математическим ожиданием, равным 1, и дисперсией, равной $\exp(-10)$. Проведите 100 таких экспериментов и изобразите результаты на одном графике вместе с корнями исходного многочлена. + +Оцените по одному из экспериментов абсолютное и относительное число обусловленности корней многочлена как функций его коэффициентов. + +\item (10) +Рассмотрим задачу наименьших квадратов --- $Ax\approx b$: +$$ +A = \begin{bmatrix} +1 & 1\\ +1 & 1.00001\\ +1 & 1.00001 +\end{bmatrix},\quad b = \begin{bmatrix} +2 \\ +0.00001 \\ +4.00001 +\end{bmatrix} +$$ + +\begin{itemize} + \item Формально решение можно найти как + % + \begin{equation} + \label{ex} + x = ( A^T A )^{-1} A^T b. + \end{equation} + % + Вычислите его по этой формуле аналитически. + + \item Вычислите (\ref{ex}) с помощью + \mintinline{python}{numpy}, используя числа одинарной и двойной точности; сравните результат c аналитическим. + + \item Помимо формулы (\ref{ex}), реализуйте решение, основанное на сингулярном разложении. Какой способ вычислительно более стабильный? + + \item Решите эту же задачу с помощью \mintinline{python}{np.linalg.lstsq}. Какой алгоритм использует эта функция? + + \item Какие четыре числа обусловленности, относящиеся к этой задаче, упоминаются в теореме 18.1 из~\cite{trefethen1997numerical}? (Возможно, их требуется вычислить --- прим. пер.). + Приведите примеры таких $\delta b$ и $\delta A$, при которых приблизительно достигаются оценки на $\norm{\delta x}$, даваемые числами обусловленности. +\end{itemize} + +\item (7) +Пусть +$$ +A = \begin{bmatrix} + \epsilon & 1 & 0\\ +1 & 1 & 1\\ +0 & 1 & 1 +\end{bmatrix} +$$ +\begin{itemize} + \item Аналитически найдите LU-разложение матрицы $A$ с применением выбора главного элемента и без него. + \item Объясните, почему при $|\epsilon|\ll 1$ мы можем неправильно оценить множители $L$ и $U$ в арифметике конечной точности. +\end{itemize} + +\item (6) Пусть функция $f(n, \alpha)$ определена следующим образом: +% +\begin{align*} + f(n,\alpha) &= \frac{1}{n} - \alpha f(n-1,\alpha) \\ + f(0,\alpha) &= \ln(1+1/\alpha) +\end{align*} +% +Вычислите $f(20, 0.1)$ и $f(20, 10)$ с помощью арифметики обычной (двойной) точности. Теперь сделайте то же самое в арифметике произвольной точности: +% +\begin{minted}[frame = lines, framesep = \framesep]{python} + from mpmath import mp, mpf + mp.dps = 64 # precision (in decimal places) + f = mp.zeros(1, n) + f[0] = mp.log(1 + 1/mpf(alpha)) + for i in range(1, n): + f[i] = 1/mpf(i) - mpf(alpha) * f[i-1] +\end{minted} +% +Постройте в единицах машинного эпсилон график относительной разности между точными и приближёнными результатами как функции от $n$. Сделайте это при $\alpha=0.1$ и при $\alpha=10$. Машинный эпсилон можно получить как \mintinline{python}{np.finfo(float).eps}. \\ +Как бы вы стали вычислять $f(30, 10)$ без арифметики произвольной точности? + +\end{enumerate} + +\bibliography{library.bib} +\end{document} diff --git a/Task 4/tex/library.bib b/Task 4/tex/library.bib new file mode 100644 index 0000000..38f6d94 --- /dev/null +++ b/Task 4/tex/library.bib @@ -0,0 +1,118 @@ +@book{trefethen1997numerical, + title={Numerical linear algebra}, + author={Trefethen, Lloyd N and Bau III, David}, + volume={50}, + year={1997}, + publisher={Siam} +} +@book{robert2013monte, + title={Monte Carlo statistical methods}, + author={Robert, Christian and Casella, George}, + year={2013}, + publisher={Springer Science \& Business Media} +} +@book{murphy2012machine, + title={Machine learning: a probabilistic perspective}, + author={Murphy, Kevin P}, + year={2012}, + publisher={MIT press} +} +@book{boyd2004convex, + title={Convex optimization}, + author={Boyd, Stephen and Boyd, Stephen P and Vandenberghe, Lieven}, + year={2004}, + publisher={Cambridge university press} +} +@book{trefethen2019approximation, + title={Approximation Theory and Approximation Practice, Extended Edition}, + author={Trefethen, Lloyd N}, + year={2019}, + publisher={SIAM} +} +@book{devroye:1986, + author = {Devroye, Luc}, + date = {1986)}, + description = {Non-Uniform Random Variate Generation}, + keywords = {algorithms generation random simulation}, + location = {New York}, + publisher = {Springer-Verlag}, + title = {Non-Uniform Random Variate Generation(originally published with}, + year = 1986 +} +@book{williams2006gaussian, + title={Gaussian processes for machine learning}, + author={Williams, Christopher K and Rasmussen, Carl Edward}, + volume={2}, + number={3}, + year={2006}, + publisher={MIT press Cambridge, MA} +} +@book{hairer1993solving, + title={Solving ordinary differential equations. 1, Nonstiff problems}, + author={Hairer, Ernst and N{\o}rsett, Syvert P and Wanner, Gerhard}, + year={1993}, + publisher={Springer-Vlg} +} +@book{hairer1993solving2, + title={Solving ordinary differential equations. 2, Stiff and differential-algebraic problems}, + author={Hairer, Ernst and N{\o}rsett, Syvert P and Wanner, Gerhard}, + year={1993}, + publisher={Springer-Vlg} +} +@book{tyrtyshnikov2012brief, + title={A brief introduction to numerical analysis}, + author={Tyrtyshnikov, Eugene E}, + year={2012}, + publisher={Springer Science \& Business Media} +} +@book{amosov2003, + title={Numerical Methods for Engineers}, + author={Amosov, AA and Dubinskiy YuA and Kopchenova, NV}, + year={2003}, + publisher={Izdatelstvo MEI} +} +@article{arbenz2012lecture, + title={Lecture notes on solving large scale eigenvalue problems}, + author={Arbenz, Peter and Kressner, Daniel and Z{\"u}rich, DME}, + journal={D-MATH, EHT Zurich}, + volume={2}, + year={2012} +} +@article{trefethen1996finite, + title={Finite difference and spectral methods for ordinary and partial differential equations}, + author={Trefethen, Lloyd Nicholas}, + year={1996}, + publisher={Cornell University-Department of Computer Science and Center for Applied~…} +} +@book{boyd2018introduction, + title={Introduction to applied linear algebra: vectors, matrices, and least squares}, + author={Boyd, Stephen and Vandenberghe, Lieven}, + year={2018}, + publisher={Cambridge university press} +} +@article{halko2011finding, + title={Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions}, + author={Halko, Nathan and Martinsson, Per-Gunnar and Tropp, Joel A}, + journal={SIAM review}, + volume={53}, + number={2}, + pages={217--288}, + year={2011}, + publisher={SIAM} +} +@book{demmel1997applied, + title={Applied numerical linear algebra}, + author={Demmel, James W}, + year={1997}, + publisher={SIAM} +} +@article{dahlquist198533, + title={33 years of numerical instability, Part I}, + author={Dahlquist, Germund}, + journal={BIT Numerical Mathematics}, + volume={25}, + number={1}, + pages={188--204}, + year={1985}, + publisher={Springer} +}