phys512/polynomial_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d6186669",
   "metadata": {},
   "source": [
    "# Polynomial fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "c9b00a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c00302",
   "metadata": {},
   "source": [
    "Let's consider the example of fitting a polynomial to some measurements of $y_i(x_i)$ (ie. $y$ measured at a set of points $x_i$). An obvious way to solve the normal equations\n",
    "\n",
    "$$\\mathbf{A^T} \\mathbf{d} = \\mathbf{A^T}\\mathbf{A}\\mathbf{a}$$ \n",
    "\n",
    "is to calculate the inverse of $\\mathbf{A^T}\\mathbf{A}$ and write\n",
    "\n",
    "$$\\mathbf{a} = (\\mathbf{A^T}\\mathbf{A})^{-1} \\mathbf{A^T} \\mathbf{d}.$$\n",
    "\n",
    "For simplicity, let's take a constant error $\\sigma_i$ for each point, so the $\\sigma_i$'s drop out of the equations and we'll just ignore them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "5473e41a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Polynomial coeffs =  [-0.91784457 -0.51906032  0.74800661 -0.42672864  0.14414031  0.03095307\n",
      " -0.85288205 -0.6404737   0.16328377 -0.54767293 -0.32993672]\n",
      "Fitted coeffs =  [-1.75823981e+00 -5.19060326e-01  4.17157825e+01 -4.26728618e-01\n",
      " -1.81910971e+01  3.09530634e-02  3.05132449e+00 -6.40473700e-01\n",
      " -3.98823737e-01 -5.47672929e-01 -2.79319819e-01]\n",
      "frac error =  [ 9.15618254e-01  5.59540506e-09  5.47692697e+01 -5.16825993e-08\n",
      " -1.27204095e+02 -1.33385372e-07 -4.57766294e+00  1.53910324e-09\n",
      " -3.44251920e+00 -5.25259157e-10 -1.53413968e-01]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABD9UlEQVR4nO3de3gU5d3/8c8mJJsDZAkkJEHCWaEKEkXFxIpQKbHFaqq1iFbBJ4pYsZyKgrSIbRULUrXWSnlaoU9/HhBFbBHRNIhWiSgIImhSkXAQSDgmCyHkOL8/ht1kIQlJyGZ2Nu/Xdc01uzP37n6HEPbDPffc4zAMwxAAAIBNhVhdAAAAwLkgzAAAAFsjzAAAAFsjzAAAAFsjzAAAAFsjzAAAAFsjzAAAAFsjzAAAAFtrZ3UBraG6ulr79u1Thw4d5HA4rC4HAAA0gmEYOnbsmLp27aqQkPr7X9pEmNm3b5+Sk5OtLgMAADTDnj171K1bt3r3t4kw06FDB0nmH0ZMTIzF1QAAgMZwu91KTk72fo/Xp02EGc+ppZiYGMIMAAA2c7YhIgwABgAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAtkaYAQAAttYmbjQZEEqPSv99VzpZLJ0sksqPS4YhyTD3J6VIA39S0/7rf0vRcVL7BCk6XgrlRwUAQF34hmxpR/Kl/A+kbz+VuqdKl9xubi87Lr0xvv7XpfysJsyUl0gv3lxrp8MMNjFdpdieUp/vSYPH1eyuLJfahbfwgQAAYA+EmZZw6Gvp85el3FXSwa9qtpcerQkzMV2lPtdKzg5ShEsKby85HOYiSeddVvO68hIp8WLp+AGp5IBkVEslB81l/+eSM6YmzFSUSo93lVzdpM7nS3HnS537nlqfb37uWW6dDgCAnRFmzkVludmDkv9BzTZHqNT9SnPpeXXN9pBQ6Y7ljXvf9l2kCf8xH1dXSScOS8cKpOJvpaJdUtwFNW2Ldpthp2i3uXyT7ftel9wh3fgn83FVhbRthRTX1ww8zg5NPmQAAAINYeZctAs3w4sjROr7fWnAzdL535eiOrXcZ4SEmuGmfRcp6eIz98ddIP3ya+nwdrOH6PDX0qHt5vroTvO0lMeRfGn53TXP2yeaoaZzb6lTH6nXUOm8S1uudgAAWgFh5lxdN9c8ZdQx2ZrPdzhqwk6PNN99VRVSVXmt52VS9zQz6JQclI4XmMuuD839wx6uCTNHd0r/fMAMOZ371Kxje0rtnK1xZAAANAph5lx1+Y7VFdQvNMxcPBIHSv/ztvm4tMjszTn8jXTkG3Pdrda4nYP/NU+f1T6FJpm9UK5u0vBfSYNGm9tOFpu9Qh17mAOVGaMDAGhFhJm2KrKjGV5qB5jaEgdIGc/7hp3D30gVJebYnNr2fFpz9VVYlNSxuxlsYnuY6/O/L8X38+vhAADaLsIM6hbTVUq5zXebYUjHC81Q07lvzfbKUqlDV+nYfqnihHQw11w8ouNqwkz+B9I7s04FnjqWCJf/jw0AEFQIM2g8h0PqkGgutX3nR+ZSWWZecXV056mrq3ZJR3dJXS6saXt4u1SwxVzqctP/Shf/1Hx8ME/anu0bdiI7+uPIAAA2RphBy2nnNAcJd+5Tf5t+P5RizqsJO55Lyov2SCcOmeNxPHatk96Z6ft6p+tUsEmWrv6l1G2wub3smDngOTKWMTsA0MYQZtC66urZ8SgvkUJrzWTsSpYu+nFN4Ck5KJUVS4VfmEvqxJq2X7wmrZwshXeo+/RVj6uk6M5+PTQAgDUIMwgc4dG+z88fYS4e5SWnJg481atT+0qykkOn2hyTDmwzl9rGvSVFf9d8/OU/pc0v1gQdV3LNoOWoTvTsAIDNEGZgH+HR5kDiuq6Muma6lHp/Tdgp3l3rFNZu38kD92+W/ru67s8Ii5bG/qvm9NXBPKl4j9Spt+Tqzg0/ASAA8S/zOdhfXKr8QyXqFRetJFek1eUgPEqKv8BcGnLRj82xOUWnBZ7jheal5+3ja9pueVX6z5Pm45B2Zu9N5z5muOnUx7w5aEvO+AwEG8OQytzmvepOHDHXnscdEqQLb6xpt/we835zFSfMW7nUljhQSn+s5nnWbKmq0rwtS/suUvsEc+mQIHVI8p1jC0GPMNNMSz/drZnLv1C1IYU4pLk3DdToy7tbXRYaI3GguZyuotTs2Yk5r2ZbhEuK/450NF+qPGnOuXPkm5r9F4ysCTPrF0nbs2rNmtzLfOxKpkcHwccwzEBybL957zjPOrZHzRWJ5SXSE92l6sq63+P89Jow43BIX/3L/D2r8/OqfZ9vXGJO2FmXLhdJP19X83z7v83f67gLzFvEIOjwL2wz7C8u9QYZSao2pIeXb9XQC+Jt2UNDD9MpYZHm3cZru+oX5lJdLR3bd2oSwR2nQk2+GVQ89qyXvn73zPcNCTP/gR+3yvxfo2S+1uEwX88/rghUJ4vN6RUcIeZEmpJUfkL63+Fmb2bFiTNfc356TZgJj5baRZpj2dpFSJGdzPAfGWsuyVf4vjb9MfN+d2GRZk9obdFxvs+vmiSddEsni6TjB82e1eMHzHWnXjXtqqulZf9jXjwQ3l5KSjHvQ9dnuNT1Uv6jEST4KTZD/qESb5DxqDIM7Tx0wnZhgB6mRgo5dRsHVzep9zV1t0m9X+r53Zqgc/ibmh6dozulqFpXU62dK21Zal69FduzVm9Ob3Pd4yq6ydF6qqukjYvN4OKZH+roTjMoSGZAuf1V83F4lOTeXxNkouLM0zqeKxW7XuL73g9sMHs4wxrxb+Pld5+9jcfV0+o5lmrzdLHHySIziO3bLJUfN+9Ft+tDae3j5lQPg8dKI3/b+M9FQCLMNEOvuGiFOOQTaEIdDvWMi7KuqGYIth4my5136Zl3Ha+ultx7zaX2/wCrK80gU1UuHfqvuXg4QqRZhTXPP15oDkKufcPPDl3NgAU0pLpKcu+rCSi11537SDc+Z7ZzhEj/ftQc23K6qM5nXmn4s9fM7a5uZ7/xbH1TMfhLSIg5jsYjqpN01yrzz+LQf6XdH0s73pN2vG8GnbBa/25XVUp7N0jJQ7iq0WYIM82Q5IrU3JsG6uHlW1VlGAp1OPT4TQNsFwCCqYcpYIWEmBP8nX5X9Z+8YP7jWvxtzb2vjuww11VlUrta8+1sW26ewqqtXaTZlR53vvSTJTXB5qTb/Iecf4jbBsOQThw+FVB2mqdoLsqo2f/k+eb+upQW1Tx2OGpuX9Kxh9lbGNvDnLKgdjDwOP30kB2EhJrTOXT5jnTZXebv377NUkxSTZu8VdKrd0gJA6TLM6VBt0lhEZaVjMYjzDTT6Mu7a+gF8dp56IR6xkXZ8ss/WHqYbCsk1PzCiO0h9fle/e0GjzPP7XtCT9Eu835YB76Uyo779tC89FNp/+enrrbqXdOb4/liqn2JOgKfYZgD08Nr/U5mPSId+KrmKrzap1Tiv+MbZjokmeNeXMm1Asqpv3OdTpup+we/9+eRBJ6Q0JopGDzce83/KBRulVZOkdb+Xrp6qnTpWEJNgHMYhmGcvZn1nnvuOc2fP18FBQUaNGiQnn32WV1xReP+d+B2u+VyuVRcXKyYmBg/V2ovSz/dfUYPk13HzLSZgcxVFeaX2JEd5hfdhTfU7FvwHXOgcl1c3aUpX9Q8f3+eeZrLO2lg98adNkDL27fJ7F1x76s1XcCp2324kn2vzPlz2mmTQjrM0NKxu9nr8KOna3adOGKOV2GQeeOVHpU2vyzlPCe5vzW3tU+UvvcrKeV2Tu+2ssZ+f9sizCxdulR33nmnFi5cqCFDhujpp5/WsmXLlJeXpy5dupz19YSZhu0vLrV1D5MUfAOZmx3MKstPBR3PqatTp6+Kdptfdne8UdO2zuBz6mai3S6XRv+jZnPeavPUl2egZ0THRp3KajMBsy5lx80ra0oOmiHFva9m/FRImPSTv9W0ff4qszegLuEdpJl7av68N79sDir39LQRQP2jslza/P+kDxaYoSZhgDT+fa5+amVBFWaGDBmiyy+/XH/
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# first generate some data from a random polynomial\n",
    "\n",
    "k = 10   # order of polynomial\n",
    "n = 10   # number of data points\n",
    "\n",
    "x = np.linspace(-2,2,n)\n",
    "a0 = -1.0 + 2*np.random.rand(k+1)\n",
    "print(\"Polynomial coeffs = \", a0)\n",
    "poly = np.polynomial.Polynomial(a0)\n",
    "y = poly(x) + np.random.normal(scale=0.0, size=n)\n",
    "plt.plot(x,y,'.')\n",
    "\n",
    "# Compute the design matrix\n",
    "# A_{ij} = f_k(x_i)\n",
    "A = np.zeros((n, k+1))\n",
    "for kk in range(k+1):\n",
    "    A[:, kk] = x**kk\n",
    "\n",
    "# Now do the linear algebra part \n",
    "rhs = np.transpose(A)@y\n",
    "lhs = np.transpose(A)@A\n",
    "a = np.linalg.inv(lhs)@rhs\n",
    "\n",
    "print(\"Fitted coeffs = \", a)\n",
    "print(\"frac error = \", (a-a0)/a0)\n",
    "\n",
    "xx = np.linspace(-2,2,1000)\n",
    "poly = np.polynomial.Polynomial(a)\n",
    "yy = poly(xx)\n",
    "plt.plot(xx,yy,'--')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13fd394a",
   "metadata": {},
   "source": [
    "If you run this with no noise (`scale=0` on line 10), you should see that for low polynomial order the coefficients are fit to machine precision, but the errors become large pretty quickly as you increase the polynomial order. Something is going horribly wrong with the matrix inversion as the polynomial order increases!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aaa27d6",
   "metadata": {},
   "source": [
    "## Using SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac127b3",
   "metadata": {},
   "source": [
    "Singular value decomposition (SVD) is extremely useful in these situations where you are dealing with a close-to-singular matrix. The ratio of the largest to smallest singular values is known as the **condition number** of the matrix, and measures how close to singular it is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "e8d9067e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range of singular values = 37411\n",
      "np.linalg.cond gives 37411\n"
     ]
    }
   ],
   "source": [
    "U, Sdiag, VT = np.linalg.svd(A,0)\n",
    "\n",
    "print(\"Range of singular values = %g\" % (max(abs(Sdiag))/min(abs(Sdiag))))\n",
    "\n",
    "print(\"np.linalg.cond gives %g\" % (np.linalg.cond(A)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "762f72b0",
   "metadata": {},
   "source": [
    "Let's use SVD to rewrite the normal equations:\n",
    "\n",
    "$$\\mathbf{A^T} \\mathbf{d} = \\mathbf{A^T}\\mathbf{A}\\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{VSU^T} \\mathbf{d} = \\mathbf{VSU^T}\\mathbf{USV^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{VSU^T} \\mathbf{d} = \\mathbf{VS^2V^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{V^TVS}\\mathbf{U^T} \\mathbf{d} = \\mathbf{V^T}\\mathbf{VS^2V^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{SU^T} \\mathbf{d} = \\mathbf{S^2V^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{U^T} \\mathbf{d} = \\mathbf{SV^T} \\mathbf{a}$$ \n",
    "\n",
    "$$\\Rightarrow \\mathbf{VS^{-1}}\\mathbf{U^T} \\mathbf{d} =  \\mathbf{a}$$ \n",
    "\n",
    "The matrix $\\mathbf{VS^{-1}}\\mathbf{U^T}$ is called the **pseudo-inverse**. We can use it to map from the data vector $\\mathbf{d}$ to the model parameters $\\mathbf{a}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "1595e75e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted coeffs =  [-0.9225663  -0.51906032  0.86120186 -0.42672864 -0.23053608  0.03095307\n",
      " -0.46748529 -0.6404737   0.01658669 -0.54767293 -0.31193299]\n",
      "frac error =  [ 5.14436612e-03 -1.24056747e-14  1.51329211e-01  2.34153658e-15\n",
      " -2.59938661e+00 -1.92566047e-13 -4.51875798e-01 -1.31741475e-14\n",
      " -8.98417984e-01  4.45976152e-15 -5.45672264e-02]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABB3ElEQVR4nO3dfVxUdd7/8feAMIAKoiBgouJNqGnelmFlubli2RbbbptambuU2Wppepm3P9O6yi7NWtfazGtLd/fqxiyz1kwlzdqSLE0zLSgT7wXvZ7zl9vz+ODIwCsggw5mB1/PxOI9zZs53zny+TMi77znzPTbDMAwBAAD4qQCrCwAAALgchBkAAODXCDMAAMCvEWYAAIBfI8wAAAC/RpgBAAB+jTADAAD8GmEGAAD4tXpWF1ATioqKdODAATVs2FA2m83qcgAAQCUYhqGTJ0+qWbNmCggof/ylToSZAwcOKD4+3uoyAABAFezdu1fNmzcvd3+dCDMNGzaUZP4wwsPDLa4GAABUhtPpVHx8vOvveHnqRJgpPrUUHh5OmAEAwM9c6hIRLgAGAAB+jTADAAD8GmEGAAD4NcIMAADwa4QZAADg1wgzAADArxFmAACAXyPMAAAAv0aYAQAAfo0wAwAA/BphBgAA+DXCDAAA8Gt14kaTPuHMMWnzv6RzDqkwT9L5m2bVC5FCG0lxXaRWN5S0NwzpEjfWAgAAhJnLctBxVllHTishqr7iIkKl/LPSznXSri+kPV9J7W+TbhxnNj7nkNKmlX+wnn8qCTO5p6Tn20kR8VJkKyk6UYpuLzVtL0UlSvYG3u4aAAB+gzBTRYu/2aNJS79XkWEoKeBHPddqk1oe+Y+Ud6qkUePWJdv1o6SrB0khEVK9YHPkRZIKzklnT0hX9Cxp69gr5Z+RjmSay8+r3N/8uj9LA2aa20VF0ulDUsNYr/QTAABfR5ipgoOOs5q09HsFGXl6L3i6OgXskg6c3xneXGrXT2rRW2rZu+RF9obSXa9W7g2atJUe/VY6sUc69ot0OFM69KO5Pn1IahBT0vboDunla6QGseapquKlWVcp/ApOVQEAaj3CTBVkHTmtIkPKVbByjEi1Ng7q/cIbdPXtf1bnXrdcfoAIDJKatDGXNn3d95055v746A7JFiCdypZ+znYfxQltLCU/I3UdYj4uKpRkkwLcr/u+6HQZAAB+hDBTBQlR9RVgk4oMaXrBUDmN+jpla6gvOlzv/ZGQsMbuj9vfJk3aJ+Vslw5skQ5+Jx3cYo7knD1mntYqtuMT6d1UKe5q1wjOx0eb6rG0U8o3AhVgk2be1Vn3XNPCu30AAKAaEWaqIC4iVDPv6qzJS7dprxGjQJtNz97VybpRjeD6Uvy15lIs/5x0aLt5yqrYwe+kvJPS7i/NRdKtkr4PDtKPRkvNyB+qyUtt6nNltOIaBksBgTXbDwAAqoAwU0X3XNNCfa6M1q4jZ9QqKsz3Ts8EhUhX9HB/7obHpfYDz4/efCdH1iYF5nyvBrZz6mbbodMKUaFhaNeRM4rb9r/S+r+a36Iq/iZVdAcpqp1UP5prcQAAPoMwcxniIkJ9L8RUJDBIirnKXLoO0RnHWd3w3CdqoRx1smVppxGnQJtNraLCpO8zpdOHzWXXf9yPE1RfevhzKer8qE/299KZo1JkghTRnBEdAECN8psw8/LLL2v27NnKzs5Wly5dNG/ePF177bWXfiHKFRcRqmfv6qLJS7cpqyjO/XTZbbOka1KlwxnmcihDOvyjdGKvlH9aCo8rOdA3r0mbFprbtkDza+Lhzc4vV0g3PSGFRpr7z54wQ1Vw/WrvDxcyA0Dd5BdhZvHixRo7dqzmz5+vXr166S9/+YuSk5OVmZmppk2bWl2eXyv3dFlwfemK7uZSWkGu5NjnHkbCmkhRV0rHd5mzGzv3m0uxvlNKtldPkTb/n1Qv1DxdVb+JuQ6LMufi6TNeCgk32x7ZIeU6JHuEeSFzSLhUz15mP0rm/VGtuJCZYAYAlWczjOLZ23xXr169dM011+ill16SJBUVFSk+Pl6PPvqoJk6ceMnXO51ORUREyOFwKDw83Nvl1l1FheZpKed+yXnAXE4dkm75fyVt3r5Xylhe/jEmHygJSu8/In33pvv+eiGSPdwMNg9+IoVG6qDjrP46a6quDfhRZw27zsquXNl1f5/2atgwQgoKlTr/3pzrR5IO/ySdPCAFBl+wBJnrhnFS4PmcX1RkXh9Ug9cI1aZgVltCGf0ArFHZv98+PzKTl5enTZs2adKkSa7nAgIC1K9fP6Wnp5f5mtzcXOXm5roeO51Or9cJmdfKNIw1lwsvPi52z/9JeafPX49zRDpzpGT77HH3EZ+QcHMSwnMO81tYkjljcsE5c/LAwGBJ5rw/3W0/6beBX7q/1/pS24m3loSZb/5X+npB+f14bHPJ7M1rZkhf/sWcyyegnnkaLeD8YguU/rTSvN2EJG14Vdowv1Sbeudfd377jnlS0w5m2+3LzHt12QLclnMFRQrPPKK2+p1+UryKDGn5+2/pN7/8pDB70Pl2Nkm2ktdd+5B5HZQk7dskbXu37Ha2AKnTXSVtD2VIP/77fFgLKLU+v7T5VUm9J/ZIP60qFepsJceXpBZJ5kXiknQyW/o5TV/vOqZ3N+1Tkcx2d/eMV6+EJuaEjsXHPXNM+mVtyc/eVfP540Z3KDnuOaeU9dnF71/ctnHrks8i74y0e72rqVs72cxru6LamQ8L8qS9G8p+f9n0UVaBHl3lVJEh1bMVaX5fqV+H2At+Fuc3QyNL/tsxDPN6MtdxL6jD3lBqVCqkHv5JknHxz9ZmMwN5eLOStif2SEaRexvX/d7sUoNSI9anDklFhfrguwN6ZkWGigybAmzS5IEdlNK9pft0D2dPmMe98P1lM8N+EAEIvsnnw8yRI0dUWFiomJgYt+djYmKUkZFR5mtmzpypGTNm1ER58JTNZt5byt5AapxQcdtb/8dcJHPUJ9dp/kE75zC365n/sCZE1df8ouuVkR+vMOUq1JanMOXpD10bK0x55j2zioOMZP5D3/QqqTDXPC1WmG+uC/LM9fmQJMncJ5n/wBfmVVzvmaPSsZ3l788/U7J9PMuc9+cCIZJuDZQWFSS7nkuw7VdYxrvlHzfxtpKAcvhH6au/ld82tlOpMPOD9Ol/l9/2zr+VhI6cH6QV/1V+29ueLwkdR36SPhylayVdG1SqzXfnl18/VXLcY1nSe6nlH/fmSVLT86Ovjn3S4vvKb9v7Man/0+b26UPSG78rv+01D0oD55jb5xzSP24vt+nZwj4qMkZIkoKNPPVb/yf3oFxaxxTpD/8oefzqjeXX0K6/dO+SUm37SAVny27b6kZpWKkRzVdvMueRKkuz7tLwT0se/++vJMde3SnpztJnaT+R8r+9UkGPfVPy3Gu/Nj+/soQ3l8ZuL3m86HYpZ5sUEFRqZPP8dlgT6YEPS9r+5wXzNHRw/ZIl6Pw6JFxqf3tJyCsqumhizwsxwoQL+XyYqYpJkyZp7NixrsdOp1Px8fEWVoTLFhBo/l9v8YXEpcRFhGrgb4do8tJtKjQM14XMYeWdmukz3lwq41dTza+0G4VmoCoqOL9dZG5Htipp232oOZpxYTvj/OPGbUratks2b0thFJn/B28USUaRHGdy9fyqDGUZJffa2mwkytlnusLtga525v+VG+Zrm5Q6bsxV0g1jL2inku3SNUS2kro/4F6DSmpx61vDGKnDHecfGCX3FpPM7UYtSx6HNtaxK/pqy57jssk4Py5jtu/SPEKNSo9G2BtKCTe5H9N1bMO82WqxoFApvtf5/UbJuvg1Ec1L2gYGmxNDutqq5JiGYZ5KLBYQaN7A1a1f5vpsXoEOHy+ZeNKQtLuoqWLC7QqpZ7v4uBdOatkwzv2YpbdLT2gpmf9t59tLtS117AtHRILrmyOUF/4MZJiBojRbgAxbgIqKzM8iwFby2eUXShe0Ll/gBX8uzjnM0dSyNHD/n0/9tEra+1XZbeuFSlOzSx6/Pdi8WW9YY/NauuIvEpxfv3Oupya+/0OtOA2L6uPz18zk5eUpLCxM777
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "aSVD = VT.T @ np.diag(1/Sdiag) @ U.T @ y\n",
    "\n",
    "print(\"Fitted coeffs = \", aSVD)\n",
    "print(\"frac error = \", (aSVD-a0)/a0)\n",
    "\n",
    "plt.plot(x,y,'.')\n",
    "\n",
    "xx = np.linspace(-2,2,1000)\n",
    "poly = np.polynomial.Polynomial(aSVD)\n",
    "yy = poly(xx)\n",
    "plt.plot(xx,yy,'--')\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}