Skip to content

Instantly share code, notes, and snippets.

@tobydriscoll

tobydriscoll/TB-Lecture-19-julia.ipynb

Last active October 11, 2016 20:51
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save tobydriscoll/dfb794e2c6891944790e628f68058ba4 to your computer and use it in GitHub Desktop.
Save tobydriscoll/dfb794e2c6891944790e628f68058ba4 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
TB-Lecture-19-julia.ipynb
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
TB-Lecture-19-matlab.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lecture 19: Stability of least squares algorithms\n",
"\n",
"We're going to explore several alternatives for least squares by means of a fitting problem:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AoLFC4lHA34RQAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAxMS1PY3QtMjAxNiAxNjo0NjozN1/ER0IAACAA\nSURBVHic7Z19XBVV/scPAoKBD1wlRVrBIBB82M10d9MtFRWy3Fq1ftH+rHzaSmpBzXZtrUQ3W1Pb\nzCzNRQNK12031wwpTAHdNBTT35qmKA+XRFHMy4PJg8C9vz+OToe59w5zZ+acOTP3+375x3W4d+bM\nOTPnc77f8z3f4+NwOBAAAAAA6E0XvQsAAAAAAAiBIAEAAACcAIIEAAAAcAEIEgAAAMAFIEgAAAAA\nF/jpXQDvoq2t7dlnn7Xb7cKRLl26vPfeezoWCQAAgBNAkJhSVlZWWFj429/+NjAwEB/p0gWMVAAA\nAIRAkBhz+vTp7t27v/LKKz4+PvJ/deLEicLCQpd/CgwMjIuL++lPfxocHKxNEZ34/PPPS0tLEUJD\nhgwZO3asxDePHz++f/9+hJDFYvntb3+rY0k85ciRIyUlJdeuXRs8ePDo0aNZXhrQC1GjQxNzgQNg\nyIoVK6ZPn+7pr959913pRgwJCXn//fcplNfhcDgeeeQRfJWnnnpK+ptvvfUW/mZ8fLy+JZHP5cuX\nyd7n0UcfZXZpZWzYsGHVqlWrVq365ptv9C2JcXHZ6C6bGGqbMWAhMeX06dMhISHz5s07duxYjx49\nxowZ8+yzz3br1k3laWtra2fOnBkYGJicnKxJOb2HhQsXurM++WT58uXnzp1DCPXu3XvIkCF6F8eQ\nyG90qG3GgCAx5fTp09euXUtOTk5NTf3mm282b9585MiRrVu3yp9JeuWVVyIiIvDnCxcu5ObmfvXV\nV/i/aWlpNATp5ZdffuqppxBC4eHhmp9cdz799FP8ISUlJT09vWvXrvqWB2CAy0Y393NuFECQqHD1\n6tWcnBzhv6GhoRMmTGhvb587d+7w4cPxUGvatGlDhw7905/+tHv37vvuu0/mmR966KHhw4cL/33p\npZfuvvvuoqIihFBNTU15efntt99Ofr+9vb2srCwsLKx79+7uzllfX19dXd3W1hYWFta7d2/RX4cO\nHTp06FB3v62urnY4HP3795dZfmnsdntNTc2lS5cCAgL69evXq1cv9eesrKx0OBwRERHO83aNjY02\nmw1/fvTRR0NDQ5VdQk4ly7w16bZQgMwT1tTUNDU1CWMd9XjalOfPn7927VpMTAx5UP3T64y7Rpd+\nzgFG6O0z5Jq6urpx48Z98MEHLv/6xRdfpKWlPfDAA48++ujLL79stVqFP1mt1p8RzJo1y+UZ2tra\nhg4dunLlSulikHNIX3/9teivK1euFP5aUlIiHC8qKrrnnntuueUWhFCXLl2GDh366quv2u128rcf\nf/zx3XffTT4P8fHx69evJ7/zzDPPhIeHh4eH//GPfxQO2u32JUuWCP3XgAED/v73vzvPIU2bNg3/\ndvXq1eQ57777bnx8y5Yt+Mjly5fnzZt36623CiXx8fEZMWLEP//5T+FXHk3knD59etKkST179sQ/\n6d69+8SJE8mZgPDwcFJHQ0NDw8PD3377bZdnc3dpOZUs59YcMtpi9OjR4eHhvr6++K8hISHh4eEP\nPfSQuxqQ07gOh2P58uVRUVH4C3369Fm5cuX+/ftx64wbN074mrZNKZxtx44dR48eHTlyJEIoLCzM\no4qVeYMkEo0ues49rW1AE0CQpFi9enVMTMyGDRuc//Tqq6/GxMQMHTr00UcfnTRpUmxs7LBhw/7z\nn/9InO3y5cuHDh1qa2sjD955551LliyRLoa0IAluur59+woHV6xYIbxLJA8//HBjYyP+zgcffCAY\nDSLrYfHixcKpnPvi9vb23/zmN6Iz+/j4CKabIEj33HMPPvLyyy+TZR44cCA+vnHjRofDUV1dHRcX\nR56KPHN2dra7krgjKytLCKwn8ff3Fzos578ihF577TWXJ3R5aTmVLPPW5LSFIBskd999t8sCyzmh\n3W5/7LHHnM8ptBoZnKJtUwpnW7ZsmWD9CIKk4dMrwvmc6Gaji5rYo9oGtAIESUxbW1tFRUVhYeG8\nefNiYmJcClJRUVFMTMyECROqqqrwkc8//zwuLu7ee+9tampyd+bi4uKYmBhStP7v//4vJiZm27Zt\n0kVyJ0gXLlz429/+FhQUhP/09NNP4+MHDhzAr2hAQMDixYv//e9/z58/X3BNLFu2DH9txIgRwgvc\n0tJy+fLlVatW4SNdu3ZtaGjAX3Puiz/88EOyj9i2bdt77733s5/9TDjoqSC99NJL+L99+vQpKCho\naWm5cOHCmDFj8MGJEye6K4lLKisrhTpJSEjIzMzMzs5OSkrCRwICAk6fPu1wOHJycv71r38JZV69\nenVOTk5paanLczpfWmYly7w1OW1RUFCQk5PTp08ffDwtLS0nJ+fAgQMuCyznhOTtT5ky5e9///v6\n9evj4+Od21HzphTOJrRURETEyJEj5VeszKdXhESji5rYo9oGtAIEScyZM2diOuIsSC+88EJMTMzn\nn39OHpw/f35MTMyePXvcnbmlpWXy5Mn33nvv7t27L1y4sG/fvgkTJkycOLGlpUW6SKQg+RAggj/8\n4Q/t7e34+3feeSc+KHhRHA7HJ598gg+GhoY2Nja2tbX5+d2YQVyzZg3+jt1u/9Of/pSSkpKSknLq\n1Cl80LkvFobAZJB0XV3dT37yE1FHJrMXe+qpp0aPHj169Oh169YJ31m3bh3+jmD5yRSk6dOn46+N\nGjXq+vXr+GBbW9vEiRPx8QcffBAf/OGHH4QK/OqrryTO6XxpOZUs89bkt4XD4RAqefPmze5KK/OE\ngkU7bdo0wRt26dIlofdXIEgym1I4G0Lo17/+9eXLl4Uva/v0OuOu0V0+XXJqG9AQrwhqKCkpiY2N\ndfmn8vLy8PDwgIAA4UhYWNjatWvx5+PHj2dkZDj/6tChQ76+vsK4DzN+/Phdu3YVFRWNHz/e5bW6\ndu363nvvLVu27Pe//73D4fD397/33nuXLl3qUWSXw5XboV+/fgkJCThUr76+/tixYwihLl26/OQn\nP8HxDgih3r17h4SE1NbWXr58ee/evZMnT77jjjtOnTqFEJo3b96OHTuSkpJGjRr1yiuvkLXhTEND\nA/4VQigtLU043rNnz8cff/y1116Tfy8CouRJNTU1R44c2bBhA/4vmWlJDkLYYUpKir+/P/7s6+v7\n3HPPffHFF+QXFCO/kuXcmq+vr7K2cIecE167du3o0aP48/z584Xxza233vrb3/5WmA70FE+bsk+f\nPv/4xz+ElQ+0n16Ac7xCkDZv3hwZGTl37lzR8f/+978vvfRSVlYW+RAHBwcL7h1hFEbS2tpaU1Nz\n2223iWYpsNO5qqpKoiT9+/ffsGFDU1NTTU1NeHi4y/NL88gjj+ABrMPhKC4uxh68ixcv3n///Z99\n9lliYuKJEyfwN+12+7333uvyJNXV1Qihv/3tbw888EB9fT1CqLCwEK/MuOWWWx588MH09HR3Eo5X\ns2NI9w5CyN1P5HD8+PG//e1vRUVFJSUlV69eVXyexsbGiooK/JmczEBEaS9fvnz58mXFMXUIIfmV\njOTdmrK2kKDTE5aXlwtfFq2wGTZsmKeXI/GoKSdNmkSuw6P99AKc4xWC9Prrry9cuHD9+vWkJmE1\n2rx5s8Vi8ehsDQ0NdrtdCN8SwEfwGyJNt27dFMfXLlq0iAz7PnXq1NChQ9vb2+12+8qVKxMTE1ta\nWvCfunTpEh0d7fIkeKA6evRoq9WalZX12Wef7d+/v6mpCSHU2Ni4bdu2vLy8r776yuVbTS6ZErkN\n5edDam1tJf+bkZHxu9/9Dn++5ZZbfv7znw8ePLhbt26dpqhwpq2tTRiGiybGyf+KCuAp8itZ5q0p\nawsJOj3h9evXhS+LBkYemewqm1K06If20wtwjrdk9ly9enVZWdn69evxf//73/++8sormzdvVjBM\nxm+gs3GDj6js6TwlLi5u8uTJ+DP2bwwePBj/1+FwfP311yWuePrpp/F3evXqlZaW9vnnn9fW1u7d\nu3fu3Lm4166trc3KynJ5RbKbOHnyJPmnsrIy0ZcFjxlpVzkcjpqaGuG/169f//3vf48/L1q0yGaz\nHTp0aPPmzb/85S/xQY/y/vXo0UPw+5MXJf/bq1cvlQunZFayR7emoC2kkT7hoEGDhKuLKsq5HRG1\nphROi6H99AKc4y2ChAhN+uabb1555ZWMjAxlThv80DsLDz7iMlyVKoKD8dq1a/X19X379r3tttsQ\nQg6HY+fOncLX7HZ7Wlra1KlTp06devHixf/+979RUVFRUVF4JiwgICAhIeHdd9996KGH8PevXLni\n8nLBwcGC7+vtt98Wjre0tGzdulX0ZWExSnFxsXDwiy++IIfnJ06caG5uxp9ffPFFwX0qeG9cTptJ\ncNddd+EPmzZtIo8L04HCFxQjs5Jl3pqytiCVQIScEwYFBQlmOjlR2tLSQkZRCrBpStpPr2IkahvQ\nEh0CKXTlqaeeuvfee2tqauR8ec+ePc5RdlevXo2JiZk0aZLoy6WlpTExMTNmzNCsrDeRXoe0Y8cO\n4a94ca7QofTo0ePpp5/es2fPv/71rwcffBAfxKsdW1tbe/TogY+sXLkSR6OdPn1acKFkZWXh8ztH\nH5HDz9mzZ+/fvz83N1eIYUNEdNaSJUuEg+PGjcvOzp4xYwY5Z7Bx48bKykrhv/iibW1tn332mbA8\nxWKxuCuJS06cOCGMu6dPn753796CgoLZs2fjI126dDl8+DD+ppooOzmVLPPW5LeFw+EQRgMREREv\nvfRSZmamc2llnnDZsmVC8RYuXHjw4MGcnJyEhATndtS8Kd3F7MmsWI9qTIRHUXZyahvQEO8SpKNH\njz744IPPP//8u+++K+f7LgXJ4XDcddddQ4YMES1xxV9+8cUXNSvuTaQF6ZtvvhH+Kixy+p//+R/k\niltvvfXbb7/F3xGCcRFCt9xyS//+/QWPyp133llfX4+/5vyitra2Tpgwwfnkgr9F6MjOnTvn7NsM\nDw8XvoljhYU1JQih0NBQvDZFEBU/Pz93JXHHqlWr3KUHJJchqxEkmZUs89ZktoXD4Zg1axZ5LXdL\nNeWc8Pr16y4Thgq9MClI2jalhCDJrFj5NSbCI0GSWduAVniRy+7YsWPYUyeaT1JAXFzc9evXRTMo\nOIhWFHjGgLi4OGGgunv3bvzhH//4x/vvv0/OGHfp0mXKlCn79+8XYs+effbZv//973jZR2Nj44UL\nFxwOR1BQ0KxZs3bv3i2MQJ3x8/P77LPPUlNThe/4+/unpaW98cYbom/edtttn332Wb9+/YQjP/3p\nTw8cOCCK6fjXv/41btw4/Pny5cuNjY3Tp08/dOgQPtLW1nb8+HGP6mThwoX79+//2c9+Ri7mHzJk\nSH5+fnp6ukenkkBOJcu8Nflt8cYbb0yZMqXTuAM5J/T39y8sLHziiSeEXwUHB2/evFmYpCFh2ZRU\nn16PkFnbgFb4ODx00BuUo0ePLl26lJw3euGFFyIiIp577jmJX+3duzclJWXBggWiV/SDDz549dVX\nf/7zn3/wwQf4SHV19f3339/S0rJnzx6tMo1qwvfff3/y5MmgoKDbb7/dXTxhVVVVVVXVDz/80Ldv\n36ioKJxATA52u/3bb79tbGz86U9/KrH+w+FwfPfddxUVFXFxcX379nX3tcrKyvLy8l69et1xxx1a\n7Td47dq1b7/91m63Dx48mN4ehp1Wsvxbk9kWbW1tNputvb3dYrFIr7yRc8La2toTJ06EhobGxMR0\n6dJl7dq1eIVZfHy8aNTFuCmpPr3ykV/bgEoMI0g7duzYt2/fmTNnevToMWTIkDlz5ki8DyKOHDmy\nfPnyjRs3iqIY/vjHP4aHh6emprr7oTtBam1t/c1vflNaWjp48OAHHnjg0qVLu3bt+v7773/3u98t\nXLhQwd0BAD9ICBIAUMUALrv29vaUlJQ//vGP+/fv792795UrV7KzsxMTE48cOSLzDDt37nRWI4TQ\n66+/fu3aNbwBl0f4+/tnZWUlJSWdPn165cqVWVlZTU1Nzz///Pz58z09FQAAAIAxgIW0devWpUuX\n3nPPPevWrcPJEf75z3++9NJLffv23bt3r2gdA2NaW1utVmtQUFBYWJhHy2UAgFvAQgL0wgCClJiY\neP78+YKCAnKHlTlz5vznP//Zvn27EOQDAIAm2Gy2S5cuIYQCAgJE+z0CAFV4d9k5HI4LFy7ExsaS\naoQQwjmGpRPHAQCgAIvFEhcXFxcXB2oEMIb3XHbt7e1Lly51jl/A+UsiIyN1KBMAAABAAQO47Jw5\nePDgrFmzoqKidu7cyT5VDwAAAEAD3i0kZ7744osXXnghICDgz3/+s7QaQbpfAAAAZ0pKSvQugmuM\nJEgNDQ1/+ctftm/f3qdPn7feeovchcEd3NY7M2JjY6ESoBIQVAJCCCoBIcT3SJ33oAaB/Pz8+++/\nf/v27Q888EBOTg6ZMgsAAAAwAcawkN5///0VK1aEh4dnZWUJe6sAAAAAZsIAgpSfn79ixYoRI0a8\n99579NKRAQAAAPrCe5Rde3v7pEmTbDZbYWGhp2oE/mIAAAARPHeMvFtIpaWllZWVAwYMePPNN53/\nOnPmTLy/JAAAAGB0eBckvHvKd99953Jb5cmTJ4MgAQAAmAPeXXZq4NkyBQAA0AWeO0bDhH0DAAAA\n5gYECQAAAOACECQAAACAC0CQAAAAAC4AQQIAAAC4AAQJAAAA4AIQJAAAAIALvEKQrLZmq61Z71IA\nAAAAUvCeqUETBi4/iBCKtARWLB6ld1kAAAAA15jfQpq57RT+YLU1ZxZX61sYAAAAwB3mFyTSWVcJ\njjsAAABeMb8gkVhrQZAAAAA4xcsECSwkAAAAXvEuQQIAAAC4BQQJAAAA4ALvEiRrbZPeRQAAAABc\n412CBAAAAHCLdwkSBDUAAABwi/kFSeSmA00CAADgE/MLEgAAAGAIvE6QYG0sAAAAn3idIAEAAAB8\n4nWCZLVB5DcAAACPeJ0gAQAAAHzidYIECb8BAAD4xOsECQAAAOATrxMkiLIDAADgE/MLkmglLCyM\nBQAA4BPzCxIAAABgCLxOkCDhNwAAAJ94nSABAAAAfOJ1ggRzSAAAAHzidYIEAAAA8ImfXheur6+f\nMmXKrFmzpk+f3umXc3NzCwoKRAd79+69aNEiBZe22pojLYEKfggAAADQQzdBysjIOH/+/LVr1+R8\nedeuXfn5+UFBQeTB/v37K7u0tRYECQAAgDuYClJ7e/u5c+cqKyt37NiRm5sr/4clJSXDhw/fsmWL\np1eEGSMAAACjwFSQysvLJ0+e7OmvGhsbq6qqxo4dS6FEAAAAAC8wFaSwsLC1a9fiz8ePH8/IyJDz\nqzNnzjgcjkGDBmlVDKutCUX10upsAAAAgCYwFaTg4OCkpKQbF/aTe+mSkhKEUL9+/VatWnXixIlu\n3brFxsY+/vjjffr0UVYMSPgNAADAIboFNcgHC9K8efNaWloiIyPPnTtXUFCwdevWt956a9SoUXqX\nDgAAANAGA6xDKikp8fHxmTNnzpEjRz799NOvv/563rx5DQ0NixYt+uGHH6R/mzA+wfkgJPwGAMCr\niCXQuyxSGMBCWrNmTWtrqxDk7evrO3fu3PLy8p07d+bl5U2bNk3it/l78wcuP8ikmAAAAJyC/UwY\nnjXJABZSaGio85KjxMREhNDZs2cVnBBiwQEAADjEAIJkt9sdDofoYHBwMELoypUr0r8F7xwAAIBR\n4F2QrFZrXFzcggULRMexBRodHa3knLADBQAAAH9wJ0g2my0vLy8vL6+1tRUhFBERYbFY8vPzT58+\nLXynrq5u06ZNfn5+EydO1K+kAAAAgJZwF9Rw9uzZ1NRUhNDhw4d79uzp4+OTnp6empqanJz82GOP\nxcbGVldXb9my5fLly2lpabfffruCS8AcEgAAAIdwJ0jOJCUlrV+/fvXq1Zs3b0YI+fj4DBgwYN26\ndWAeAQAAmAkf53gB0xAbG/te7qFx7x51/lPF4lGQ8BsAAC8kNjaWjALnCu7mkNgA0XcAAAC84aWC\nBAAAAPCGlwqS1QaR3wAAAHxhckEC4QEAADAKJhckd8AOFAAAALzhpYIEAAAA8IaXChJE2QEAAPCG\ntwoSuOwAAAA4w0sFCQAAAOANLxUkSPgNAADAG14qSAAAAABvgCABAAAAXGByQSLXG5HZVCGoAQAA\ngDdMLkgSgCYBAABwhRcJ0tioEL2LAAAAALjFiwRJBKyNBQAA4ArvEiTYlA8AAIBbvEyQQroJnyER\nOAAAAFd4lyCRQMJvAAAArvBeQQIAAAC4wuSCJIpcgDkkAAAAbjG5IImIDCHWxkKUHQAAAE94lyCR\nMFgYm1lcnVlcTfsqAAAA5sBP7wKwIzIkMIKhy27cu8cKy2oRQpW25iVJA5ldFwAAwKB4l4UUaSHC\nvmnuQJFZXI3VCCGUvruisKyO3rUAAADMgXcJEjNEMeX7Smv1KgkAAIBRAEGiQuaRDlNHYCEBAAB0\nCggSFUQRE4VltaBJAAAA0niXIHUI+6YWZecysg68dgAAANKYXJB02fRonytjCCwkAAAAaUwuSLog\nxNeRUA3qAwAAMAFeJEjOi5Bo2E9WW7PL01ptzbBIFgAAQAIvEiQ2iMwjcpvarOKLzIsDAABgGAwm\nSPX19QkJCR9++KGynzNIrkqqzoyRYWSOBvDaAQAASGAwQcrIyDh//vy1a9f0LohbSAtpTFQvUVwf\nhDYAAAC4wwCC1N7ebrVa9+3bN3/+/I0bN2p4Zs0TfotmiWaMDIu0BJJeu6V5FdpeEQAAwDQYILlq\neXn55MmTtTpbpCWQXiw4GfA9Y2QY/vDkyH6C2QReOwAAAHcYQJDCwsLWrl2LPx8/fjwjI0P+bxkL\ngMhfhz+QFhL22o29+ScAAABAwACCFBwcnJSUhD/7+XFdYNL2EnQIe+0ErdpXWmsgQbLampfurogM\nCYQdNAAAoA3X/TttrLYmpJ02iCaQyIi+SEsgKrt5UePsVFtYVjfu3aP4c+aR6orFo/QtDwAA5sbk\ngnS+6jy6pTd5JDKkG6U5JJcTSDcvyiKHnubM3Pat8Nlqa5657dT7yXE6lgcAAGXExsbqXQRZGCDK\nTg3ht4ULn8nd+WhDKhBCaEw0MY1kkLiGce8eE2lnZnE1hK0DgBEpIdC7LFKYXJBYQnbfLPdKp0Fh\nWZ3LjHykzQQAAKAtXi1IlZp6zyRMH5HLjnOvndXWLEwdOf9p5rZTjMsDAICX4HWCRC97ECkzIveg\n6KKcxzWIJKcgZTgZuQ6OOwAAKOF1gsQG0RwSYpJGTxMyi6tJZ93YqJCxUb1EsQzguAMAgAYmFyRm\nzjHRhZzlJzLkR5uJ591jyVjBSEtgQcqd+AOpSZCUDwAAGphckKTR0HXW6amMYiGR5tGSxB8Xw84Y\nGdZxKw3Y2wkAAI0x2Dqk8ePHKw5bxG40Z2eaJlhtP0Y0uNQe8rqFZXVLaBRCNaKAC1KBEEJjo3oJ\ncuUyBs8cWG3N1trmfaW1hWV1kZbAJ0eGGSizBgAYGoMJkiEgvXMCZCA4t0uRRDIjUtYx0SFo941s\n5ZwHCioD50nqkHGjDGUWV7+fHCda6QwAAA2822WnXa9a2SHEzpWFxHBZrmIkkk0gJ+PSfNNI49Yf\ndbnN/NLdsGkIALDAqwVJQzqfQzLCUiSX2coFRHs78RyaoYDM4mp3jQKrrwCADV4nSJRyKHRYhORq\nmsoQS5EkJpBuHvxRpTKPmCquQWQGRVoCySYrLKvlcwwBAGbC6wSJEuS0kDvN4zzQTiJbuUCHpHwm\n6qBF5lFByvCKxaMK5g4XjuDpJT2KBgBehFcLEqXgAnfTRZwvRZKeQMKYdRopq/ii8BmvBUYIRVoC\nyXqAFBUAQBszC1Jrx40nMJSCCzp12SHuLSTpCSSMKaeRRJlknxzZT/hMrsRCkKICAChjZkHiDdFS\nJB1L4hKJXHwk5ptGEi3yJa0i5xQVLsPwAIANmcXV5n4CvUiQ6BkoInVxdyGelyKJnnKJpaDkNJI5\nIO89PVG8U/vYqJCOKSouIgDQg5nbTgn/9C4LLbxIkOjRaZqGm3/idykSuY5KehGoKH6dQ1PPI5bm\ndQhVeNJ59ZUlcEnSjypl4hQVAkvzKpbmVQxcfnDmtlNmCl0xOsLIycTTmV6dqYHGy+YyTcPNP/G7\nkbn8MHQ8jST0y/tKaw2dWYf0Os4YGdZp2ieEUGFZnaFvWZrM4ur0m/GEmbbqwrLasVEhkD9Jd8zt\nqRPwOguJRi67TtM0uISrMQ75uLuLaBAg+yau7sJTRNHezuYRRhTKITKqzITzEmA8bTbu3aMmvmtD\nIPIVmyOeyBmvEyQayDQvRGst+UEkKi6XxJJ0WI3E2WSYR4gi3SWMADL0ztC3LI3E5ET67gpDDz6M\njshXrKwtBi4/OHD5wYqEVwcuP8ibkwbj7YKkeatIW2B8LkUi58CQDCPPNNNIciLdMaRIG/qWJVia\nVyHam9H5C2xLBNzA2V+nYFSEM5ZZbc2tt/S22pr5HBybWZDaurlYh0QDUtUopSaiipwlsSSmWY3U\naaokAdPcsjsKy+rSiVQUY6NCClLurFg8igw7LCyrNaUS849zbKdZR0VmFiSX0BgXyB+t8Dn7It9Q\nEODzRjxCTqokEhPcsgSiNb84sBBvB0XWjPk2ZrTamnFUod4FkcJlbKfIsWEOvE6QaCBzSSni1X6S\nbygImGAayVO70AS37A5ReHd64kBBfSMtgUsSTRv1brU1D1x+MH13RfruioHLD+pdHNe4i6/zdEkc\nnwmdRXiLILkb/2reSJ3MIRFyxUmnJnNVrwieQ9hl4qldaJqZM2dEU0fkuivkNH/GuTHhEWQQBxYn\nHQvjjn1unjROOhBt8RZBooeoO5bu0Dnsx2Wu6hUh+qYRe2dP7ULRNJJpnFeFZXVkVZCpkjDiJLOm\nyBeFnII4EK8bX5EWklY7F/MZ0YC8U5C0bQw1NhYP/Xil5/465y8bzp3t6QQSy9W8eAAAIABJREFU\npuM0kkmcV1kyqoL02pkjoZ8oiEMgs7iaKxNQVNWihuChA9EWbxQkbfHIwuBwKZImTstKPqw9+Xg6\ngYQhV86apjsgldXttiOipcEG3xrKamse9+5Rd39N381RjIPoQTV9tKe3C5K2Q3uJvEEuv8ODYaEg\nxA5j6KgzZXdtyu5ATq4KdDPuTviJ4VqcROSXS08cKHJUZh5xu589Y5xTqHTcyNiDVuCht+kUMwuS\ny/2QkDzZkI+yvEHOv9UL+SGCIgwddaYgsPDmlw0sw87Iz/IeGWKS/EmZxdXOQRwzRoald/SG8TBH\n6NQ6IUg0D220965TzCxIbPDU5cVzj+ZRoj8OAzRkomwCCWNoGXaGdL51kuXdEmiO/EmkEyzSEliQ\ncif+jGVJ+BMP76YoBz9+UEV+Y2Vn1nZQriEgSGqRs1cst6jpmkXw8ALLRNkEEsZkwd8y/XUYc+RP\nIp950Y7AHfp6DhRXzmDXoK3gDm8RJHcjAvVOM/LBlbPuldshtkeeK+Q0oWIIDzWG9Nh4OoYw0zSS\nfH8dRhSSY8R7d+kEE+BttOFysGumJ9AZbxEkNvC8BZ9LRO4LNafiYT5MJuR7rmADXJ6drh5BLvWX\naSmSAzsj3rtz0Br5V976eneDXWVPoCHeUG8UJI3XIXnosuN26kWBv9GIXbOnZoEz3Nq4nqIg1LBD\nrJ0B771T45h8HnRfAuwu9Eb9E8jb4hMBbxQkDVGvKPpqEvl+KsizZ8SuWc0EEobbIYVHiIRZtoXE\nl1PLUzqdMzPEUjNRKxj3IXTG2wVJ5bJQ0c/ljDtE39E346Hi6GeMEbtmNRNIbk7IY5/VKcqEmTen\nlkfICeERv576zYxKJJnkqg/RFm8UJA1j4TRJBKcXHmXhk4Mhumb1m1cZulMW8GjTehIj+mkxMjW4\nY8ZCzzJqa4h032KCJ9AlZhakVlYb9GHkh/Zzsm+s+mxshgu0E2mwAqPw5g+N2iljlPnrMEb002Jk\najAn82TSSSYVDB8NYUiZWZDkoNLRpCZNA1dokkWY/zAekQYrbjLjdsrOeBzub8xpJOmAbxJOblBa\nPzrma+D9vZOPH/tL7tmzJycnp7S0NDg4OCYmZvbs2REREdI/yc3NLSgoEB3s3bv3okWLZF6Ukloo\nexTGRvUSekYdHyZ3+6x4BHkvhWV1S9SfkSaipe+Kz2PEyTMSNeH+2CwWGn1faa2CSEX2SAd8k3By\ng/LDdxU8gdwu4WctSMuXL8/Ozg4ICIiPj29oaPjoo48++eSTd95551e/+pXEr3bt2pWfnx8UFEQe\n7N+/v7IyUNq2VVkb69iddViOo/SVGxMdgm6mn+HfViDlX8N3srCszhCdsoDijLoYY41CMB45qHm4\nQekV93zuPa0epoJ06NCh7OzsAQMGZGZmhoeHI4Ty8vLmz5+/ePHivLy8wEC3VVxSUjJ8+PAtW7Yw\nLKwslM2Qc/IwaaIfxrIVVIa5C3AyiFaMyuhKY41CMB4lSeLhBqXbiMO9pzWB6RzSxx9/jBBauHAh\nViOEUFJS0n333Xfx4sUDBw64+1VjY2NVVVVcnHgjSx5QuiqNi4dJZa/kEs5nFDS8ZePGNahPYMjJ\nLIt8PE6SZLQblAP/40XEWJAOHTrk6+s7ZswY8uD48eMRQkVFRe5+debMGYfDMWjQIBpFUqkHivdu\ncHkGlkiscvAIAwXaaZhJFhk5rkGDpcFGC3z39JZ1366+09fTWJ4J+bATpNbW1pqamvDwcJFrLioq\nCiFUVVXl7oclJSUIoX79+q1aterJJ5985pln3nzzze+//15xSShlnJM/JyH6pi7Pk7IVVJ3Cf6Ad\nRr1FaNzV8posDTaWgahgzkzf7eo9fT09ffw4mTVwhp0gNTQ02O32nj17io7jI/X19e5+iAVp3rx5\n2dnZNputqKhow4YNkyZNOnjwoPQVW7tZVJeaCjwstJZe5eARRumbNMwk63wGA4Xeql8ajIxmICpw\n1SrIuqshnb6exn38pGFqISGE/PzEYRT4CP6rS0pKSnx8fObMmXPkyJFPP/3066+/njdvXkNDw6JF\ni3744QeZV9++fbvSgrtFTaYD3RctafgEG6VvUhla5oyx3FYY0YhB8VjEQLMsyly1+t6gnNfToz4k\nNjb28OFDKkrECHaC5Ovri1wJDz6C/+qSNWvW5OfnP/PMMwEBAfibc+fOffDBBy9dupSXlyfz6lOn\nTnV5XI2nRU2fTiZr0GXeRcPe2SjubM3LZhTTkETkC9Jq7pBnPVa2+ExN5ahH820/S0pKbr1jmPrz\n0IadIHXr1g0h1NQk7nzxEfxXl4SGhjovOUpMTEQInT17VkFJaCwKU/Ps6jLvoj4cwx18ds1aJQ0i\nMYppSEL6LbUKrUS8NjpG8cBRxxRfcrb91H1QSwN2ghQcHNy9e/dz5861t7eTx61WK0IoLMztyMVu\ntzscDuezIYSuXLki8+o0REjNQ6C7y45EZeXoO5aUiVZJgzqcxIBxDRpaxqSByLMeK75lHS1gTye9\nPB3UcruVKNOw77i4uOvXr588eZI8ePToUYRQfHy8y59Yrda4uLgFCxaIjuNIh+joaPWl0qQfkZ9Z\n9eb39UxFpW0ANOImXawEWiUNIjHixLKGK7H0nfaXj+Jb5twC5n8UqACmgoT9bKtWrRKOVFdXb926\n1dfXNyEhAR+x2Wx5eXl5eXl4bikiIsJiseTn558+fVr4VV1d3aZNm/z8/CZOnMiy/CK0yqyq78ha\nE78Neft89suUSmWUeRSMtgMRQ8Q1qLllvSZHZa4RNGV+VaaClJycHB0dffjw4alTp27atOm11157\n+OGHGxsbZ82aJcwSnT17NjU1NTU1tbGxESHk4+OTnp7e3NycnJz8+uuv79ixY/369ZMnT66pqXn2\n2Wdvv/12BcXgYWQR0aEHZz340jxJOf9xDZqH2GGMFdegfkksiSHiGjS0jJm1r4I1gnJeOj5fTBFM\nc9n5+/tnZWUtW7Zsz5492HEXFBT0/PPPz549W+JXSUlJ69evX7169ebNmxFCPj4+AwYMWLdunb7m\nEVIZZaerD1fz8RQPub+koZEnCRnhxkk03y2348l5zLKq7iXtkLHQamtCTDIWylwjqCoZI2T7xvTp\n02ft2rUSX/jFL36B54dIEhISBJ+efNpukbVBn7W2WYtJFA9z+OtqUmiS55uEcwtJqzxJRkeTJbEk\nHTdS4VGPNbSMmUXDyhRRTlJiaou3b9CnBg1fb8adOO3HlzfnFaU8Scgg8ygYrZbEkvAf16DSMtbF\nJav5IiQD4aWCpEmvpKZb1zdAS3P/FeeR35QW3yDub5xEtEqBxtwhb3qsPohDl0A7OYuQEPduCWV4\niyBxmExQr46M0rPLc+Q3pYgGDM83TqJtRAPGEHENGGUDEV2WmikbL0qXzSiK5S2CRBsFQQp6LbQW\nWWNa6aJRDAXNMUqgHVVVvnkJvm5ffTpdIy41MzReKkiaiIGGrl6W2YMoTaiQNcBbx0QpxA7D+fJJ\nAUpjZJ7zNWiiwYxNQPkBOGrEktvho5cKEg/o9UxouPEECbf9suZpKUQYwpVPI5Ufhue4Bk0GIowt\nYI/Gi9zqimJAkBSiZu+JGz/RyaRg4Hbgtl/Wau5EAt6sQwyNVH43TsVrXAONgQiDkZZH40Xz5Vc1\nsyC1dpO1DkmZu8y43mRKQaU8bIPrkn2Uu0hDbOJOI5UfxhC3r8YiZGwCKu5YpPsxo/RXZhYkztHL\nxyUzqNRTRAHQ/LwADCbzSfjcxJ1Zc/Bz+1ptEKyjCdjpeBFcdiZBfUPSW2tJG3o7IRklAFpz+A+0\no6rKfN6+VrfMeKkZWWyPxovyxxz+jXJ37WGPlwqStni698TNX+k/Ge4N68CphthhuA3oEPDCOEMN\nb5nlSMuj8SIPfYi2eIsgSTStMm+G5q4JNs8T1axuHI6U2WSx47xfYBxnyEPTa9vuej3YGo4X+Zzb\nc8ZbBEmEtpaBgdbcUfU0cjhSJu+XknnkDA89sjtoVAKHcQ3aPufMHmxPY3f13cWGBl4qSOrRRD/0\nnXxS5mk0FlrNbEvDYY9MwqYSBHiIa9A2eyEzC1gUnd8pyuaA/ZtgDolj1D9hiu0t9ssIqPZNuiT+\nkg+zCTMeemQSBnGGHHprBWjs/KTtCV3CzKDnChAkhWje27LvxTR/UTlM/KU4ZslTeO6RGYwMePPW\natvuzCxgT7dy5nzyUgFmFqRW9xv0ads3KT4b+4U7tDto3jxXDELsMLz1yAL0kgaR8GYcU213emNH\nlZ2A7tWuHjMLElU06XT0HeDQ3kZdd8+V+vRO8uGtRxaglzRI4rT6Gsc0QivZW8ByHBjyq133l1Em\nIEgaQLtn1xDam1Fy5bnydIpYDVz1yCQ0tkFyCT8bIxl30boCB4aCG/SDhbG6I9H5KrN1NOnZGft5\nGFgMXMWh0kvg5hJ+emR3UA3rIB8nfcciNDYI5tYla7L8qt4iSCIMZNNoCAOLgauKZWym8DkYZxbW\nwWfWD61KxcYly2zKk1u8VJBUotXjqOMckjdMbjNOq8rnFoVeGNZBQ4MZuGTVvywSE0X8+JClAUHS\nAK2GxrS7b0+DShXA1VQKYznkp0cWoJ00qMPJuUkgREmDaUeQil4WmY3Fp12uGBAkJd2WVv0s4+6b\njTzwFvmNYeAA4co6dIZ2DTDOiu0ONhpMI2hNWSxGh6fOIGaQBF4qSCo9yxqG8eiV2Z6BCwvpGmzK\n0jhweQkeegfGSYN4239EWw1mGUGqLK2XzDGQf5NNwcnZYFpBYjY+VZkRzmRBMoizyG8Ms/lh3gLt\nGA9BeGh6ehpMO4JUmUedaqAKe0wrSEaEtj3hVTE8jI0DxheSCb3NGF3CwywavZV2tCtQmUlNloqT\nmUs1gCAh5Lk5pWF0ALPsQcxMRh56JRHMIpK5CrQTLztjG5at1xQa+chpaz3QniMkT6i5OcvbjKY7\nvEWQRLKhUkU0VA5mkd/KYngUwEnCR2brb0i4EmM2SYM6XKKj5ukiyfTcALTnCJU9MApeN9h+AnAL\ns75Sl63qkH6Ggi7+Sa4C7RgnqkAc7AtFO5KF6t1p4l81iiXkDhAkhNQNdlR6Qpi5gBksQhJOrvtU\nCpudy53hKtBO9zA/fRN60h6FUL07+b0KV4+cekCQlEAOQ4wS5cLySdU9/FfHeEV+lmGxj/JHegfa\n0Y5koXd3tPNM6u5Alomf3gWQy549e3JyckpLS4ODg2NiYmbPnh0REaHmhJGWQMXmrYaty24OieaU\nqYhISyAqo3qFTmDvreq0GOzRxWk5JjoE7a64UQDmnSDtZPb07k5Nnkk1XRlvGMNCWr58+bPPPpuf\nnx8cHNzQ0PDRRx89+OCDX375pcRPmBkE2gaD0nuwWPYOugeb6ei44GEtjvOl2TktdZ1FoxdixxJP\nRw9mWstoAEE6dOhQdnb2gAEDPvvss23btuXm5r711lutra2LFy9ubtYptFS714yZC5jleFn3YDNd\nvFU3Lqf3vd+4tE4dk75TGrQfcnpyq9UUrxyj3B/2Q1LDxx9/jBBauHBheHg4PpKUlHTfffddvHjx\nwIEDik+r1bBCvWdA9xAA88F4QSgJJ1HvOjot9UpXwSBZlChgR0O5VRVXZaIOxACCdOjQIV9f3zFj\nxpAHx48fjxAqKirSqVBawsDiZuzA4Sr6Wd99enSLetfPacnDTn30fAAMAnY8fWLNlF+Vd0FqbW2t\nqakJDw8PDOzQSFFRUQihqqoqOSfRtv+lGg9DaRqc8abO+vpt2KdVFV2Oh0A7HZ2Weo0A2CSLonRm\nrdZxuxv8GSXqgXdBamhosNvtPXv2FB3HR+rr6zW5ikcyoHn3yix70I3LqcsGKxMeOmXEQco+vQLt\ndMxbqNcsGhsNZhCw46mT2VMB84NMDYppbW1FCPn5icPT8RH81045X3U+NjZW87IhjYZLDGYddMk0\nKsC4U9b3ZhEHgXZ6hdjduBw3O/XRgJLcqglYN1N+Vd4FydfXF7kSHnwE/7VTwm8LLykpER1U/JZq\nPt5nHJ/KxqOie6eM0cV9pHugHWMPrQi9UnWwMQppDB+N4k9jAO+C1K1bN4RQU5P4rcZH8F9dwsZN\npIn7i8EAx4ipJRSjS1pVd+jS15A2oi5OS/apOvSaONRksKVh4mOjaxvvghQcHNy9e/dz5861t7eT\nx61WK0IoLEybeFaPZm6YJYXTEPbjdB2tBN23feIh6bWALjaivvYx1UanYf+ptGg7bWIDqRTvgoQQ\niouLu379+smTJ8mDR48eRQjFx8crPq2+0cAkDOaQ2PfRei3H4eHd0z29rO42IvvhCMuJQ6r2n3qn\nS6evACyMVUViYiJCaNWqVcKR6urqrVu3+vr6JiQksC8P7UA4qht/6QWzYTL7TYBcom96Wd1tRBI2\njx/LMHfNw2JVqqlR/DRyMIAgJScnR0dHHz58eOrUqZs2bXrttdcefvjhxsbGWbNm9e/fX84Z6AU6\na2Jm0d/4i9HWfKKr6P6e6NgX6+iz0jfE7sZFeXJaag5V6199l2LotbEGECR/f/+srKykpKTTp0+v\nXLkyKyurqanp+eefnz9/vlaX8OipohEgQLXX0CvmShcrQfeYb0xEh0E04yk0fXZiJGG/OtjQqRrV\ndynSj7qBJMoY20/06dNn7dq1epfiBjT6l8iQbsJDabU1IWo+B5Y9FLkJBbNXguUuGxIwTqBHwokk\nk9BeiMY4xE5zC8noi4c0xAAWEiV0Dwh2h+Zv7z6dHCa6xDVw8m7ruDiU9p5AMtHLacneKFR/d+rN\nO9PsQOG9gqQJWg2EmWUPYtlD6RL5zcl8vo5TaJzsCcSy9RkbhTzMj0ogPaLlueTIxILkYXo6D14Y\nGiNQqsaE7kHAiF2olf7z+T9eXadAO14kmWHGd/aZZDVsXK4eWt0xrSB1io5efme49R+qQaTWDDRJ\n35Q5InTxWemb6Vzi0lTtfvYLGzRsXE08bJ0ENRjHiee9gqQMSntPUM0epNeQmf0mFJV8GAcYHQPt\nMLrXgC479el+12pQXHjTbInkFYKk4TjRcI1NdfemTmEd+8tT6+hignMVYsdmZlQXo1DDGTLNm0za\nXmSz+4xivEKQOkX3XAb05pD46aMZbEKh4650zugSaMdzDTC4IjPziNKt8ZPPTC+8V5CUtT2bWQoN\nH3F9p1UYz6N0iDfRe45Ql1gsrmqATaCdLkah6EJqnm1NYo5MMwntvYKkHg2NX3rTLfpOq+j4nvAw\n2GQcaCd2z+pdA2wC7fQyCmm4oxWPIaQnofXatlgBIEi8QGlwp6/LjuVelvwEmAkwNhA5SSxLFoDB\nNJLu/nakrsfnZCEzJ4Ag3UDmY01vMyQGo2n2jzvLxSgknIRaMV4aTD6cM0Zqs1WYShjbiCzbXZPR\nBo2Yo06CGvQepkhjWkGiNRzjJkZAJvquimUZ+c1VgBmG8aw+hw8nbRuRQ7PYI0RGrWJMY1qZVpA6\nRf2zq+1DwMC5ocssN+PIbwyf7ydtrx1XIXYYlouxGJvFmpu/GpafBzemMrxXkEToPrqkNJrmykNN\ndXKVhwxJIhjvwsBJ0iAS2mMgHc1iTdzRWpXfNNtPeIUgadgR09gMicbZMPquisUwm9jnsDsWQVmP\neUyJRnsxlo5GoebuaC27KeMokAivECR3KHhpSdtc29EfjYA08nXVq49mM7HP7VzCkqSBwufMI9US\n31RJFlEDnEQ0ICcbUdu4BpFdwv4JV39rHQQ1WlX5ZT7zurtJpPFqQSKR6U6h55yl4bKjFxMoHzaB\ndhwGmGGYxTVw5ZsloWci6x7mrv7WNGw1c2yJ5NWCpHJlK9XXXqOdKPW33NkE2vFwp3Kg57TkZBsk\nZ+iZyOQEjC6jEJUhG/Q86iLnsIFiHLxakEh0X8xMo+PmZD9vBhP7HAaYYZjFNXA7hcYmp58uRqFK\np71WMd8YftzUajCtINEYFNCOEdD8hJzs501CSfi5yuEmAaXb53YKDdHM6ad7XKWG7mhmYwiurGdn\nTCtIcvD0PaHtF9J8WTsn3TTtQDvecriJYJxAiKspNAylfA26G4UqvRraxqybY0skrxAkOYMC4zah\nTHTspmkH2uk+uS0NgzjDfax2pFUGDUnmxCjUKoZQ29fTQJNGIrxCkNzh6UNAeysH8r1VL5B8Lkyh\nDVfTJxgGcYbcTqFhaEgyt3GV8mHmcuTQde8OrxYkEk+7CdobL6rvtvTdCYmEdo/MYRY7EgZxhro7\nr6ShEdfAiUtDK+NPvUed88khmYAgeQDtSDx6vh199y2m3SNzbh8g7Rw7LuHfFKYR18BJo6t5Z7Vd\nOiZzZT3PIT/IywVJzZiCw9deBFd2A9UemX+oZs41xCpIqgE7HBqFncL/MEIXvFqQSOSMbuhH2dFa\n1a974BnV943/volqvgbd14fKQdu4Bk4iGpAKbyTVYQQENXCHHIFRY70y6OJVvre6r9IgIatL29Bn\nfvomCZjt1Kf7yMMd9GpA3yGIJgufNbkFiaY3kD6ZVpBoQLtdtXW1c7IICcOmR+bTPEKUsxVwNfJw\nh7Y1wJU7mkT+NDPtWzCQCJGAIN3A0/aj9OZTypCo+8CZns+K276JhGoCIf49lkjrwRYnEQ0Y9d5I\nTV5PyPZtGNzZBx6vQ2Ibzq8mqI/zKVMNTQRO8vV5hIbhmobwWGI0jGvgSoOVWf+07VpOwuI9xSsE\nSSYejdwpOcG0CsfiZxGSUAZKxTDKoj9KCYQMEdGA0aoGeNNg9cvstOpMdK8K9fixv+SePXtycnJK\nS0uDg4NjYmJmz54dEREh/ZPc3NyCggLRwd69ey9atIhaMV3AwC1Lw7Wl+xASExnSTbijfaW1YzWy\nZrgaLEswJjoE7a7AnzUUUXKgzbk3hkYN8NDizsvs5AgDs/2rjDWZxFqQli9fnp2dHRAQEB8f39DQ\n8NFHH33yySfvvPPOr371K4lf7dq1Kz8/PygoiDzYv39/lYVRtYk9nWdI5Q4rAhwmN4u0BKKyG5+1\n8idw7pkkcZ7V10SSye6G24gGjFY1wOGs4dioEGFkIGewRem5Jcd8VlsTMogHm4SpIB06dCg7O3vA\ngAGZmZnh4eEIoby8vPnz5y9evDgvLy8w0G2rlJSUDB8+fMuWLVSLJ3NoQxUankBOBs40jD/SM8nD\nYFkCHNcgdFua9Beifo3/Goi0BKpveq4iGjBjo3oJpSosq1vS2fcZrGV2N0+pexcnDdM5pI8//hgh\ntHDhQqxGCKGkpKT77rvv4sWLBw4ccPerxsbGqqqquLg4j64l87mX3zy0N0O6cVqNem0OQ4FpRH5z\nOFiWiSZxDaJ+jf8a0CSugUMnrZpnW8Nb4P8B6BSmgnTo0CFfX98xY8aQB8ePH48QKioqcverM2fO\nOByOQYMGUS+fJLpErSie++VqERKGarYCxI0hKAHpyck8Ui3xTZkYKKIBoz6ugX8nrZxn27gDKdqw\nE6TW1taamprw8HCRay4qKgohVFVV5e6HJSUlCKF+/fqtWrXqySeffOaZZ958883vv/9ek1IpW/dD\n7xmiEY3GZ0+tSaQZh4agBOQ4WhM4dF5Jo95K5i18FCN6xTx6tjV8PV3u0Wes+G92gtTQ0GC323v2\n7Ck6jo/U19e7+yEWpHnz5mVnZ9tstqKiog0bNkyaNOngwYMyL61JkzNLYal+bSyfo0gaWsuh90YC\nzfM1GOv2kRY1QNoW/Nyypwufych1SgMpYwXXCTC1kBBCfn7iMAp8BP/VJSUlJT4+PnPmzDly5Min\nn3769ddfz5s3r6GhYdGiRT/88IOcSyeMT4iNje30a/Ld+sx2c1A208DnKBJpnfKZT92VQNRtqawB\n3pbjyEF9DRjCKJR+bUVSwVJW/RuvyOkJdYRKlN3ChQtbWlrII8uWLfP19UWuhAcfwX91yZo1a1pb\nW4Ugb19f37lz55aXl+/cuTMvL2/atGmdlid/b76715WMRZamssOsDMWXX314dCWvA2dPg5Gk4VZ3\nZaK+BgS4amVpyGcg80j1kqSB8n9bWFbHrVEo/9kmNRVp+uh2amyF3xZeUlLCsyZREaS9e/c2NjaS\nR1588cUePXoghJqaxMYsPtKtm1ubIzQ01PlgYmLizp07z549q0FxOUP95D+3XmNtl0by6b2Rhuy2\nVNaAQSfGOzwDtmarzYO1FjxHFcp/tumForjco88Q22UJUBGkY8eOuTzevXv3c+fOtbe3k/aQ1WpF\nCIWFuW0Yu93u4+Pj4+NDHgwODkYIXblyRWVRXU4DuoRZL6/V2lgMVxENzklWtFsSyNFtSkB2Wyox\nhPPKGVFLebT4j+eoQvnPtkEbjg1Mw77j4uKuX79+8uRJ8uDRo0cRQvHx8S5/YrVa4+LiFixYIDqO\nIx2io6PplLQTqHZ/6qO0uY0903Yvc25vUwIN4xq4dV5JI5pGWprngTyT02a8deWimB2JZ5tNuAEE\nNXROYmIiQmjVqlXCkerq6q1bt/r6+iYkJOAjNpstLy8vLy8Pzy1FRERYLJb8/PzTp08Lv6qrq9u0\naZOfn9/EiRNdXkihp8uTqUh6aOCy47ir0nBWn+fbdIdWcQ1GjGgQIFcjyfcBiG6ZwxaXE7MjGoJo\n7LLrbKDMLBpLMUwFKTk5OTo6+vDhw1OnTt20adNrr7328MMPNzY2zpo1S4hZOHv2bGpqampqKp6F\n8vHxSU9Pb25uTk5Ofv3113fs2LF+/frJkyfX1NQ8++yzt99+u8oiKRtZsxyPa7uZm+6QXaeZUj4r\nQ5PG5bBrluZJoheWbyaK/HUctricZb8sI3GMaCQxzWXn7++flZW1bNmyPXv2YMddUFDQ888/P3v2\nbIlfJSUlrV+/fvXq1Zs3b0YI+fj4DBgwYN26de7MI8VID9aYbXOgMuUX58HQHSfttMkey9t0gjSa\nxDUYNKIBI3rCZaZ+FwWncYicuAaqkTiGexKcYZ3tu0+fPmvXrpX4wi9DAfBKAAAXw0lEQVR+8Qs8\nP0SSkJAg+PQUINFOyiZsaCfjUZO1l/NgaA1n9QWMEtGA0aQGjD4xTj7hMsPfySHak1wOQeQ421k2\nHA4Y0XA3SAbABn1yYbZ/iQhPnyduFyFhtJrVN2JEA0aTGjDi/BkJufxIjpnoNIFkAA122bK0G056\nDMrhCFUECBKPqNk3lttFSBjRrL5ijNsjq49rMMH8maeqbAgPbacty7jhjLUCCePtgiQzpI3x9KCa\nQDu9LDllKAszM3qPrDLttyF6Z2k8VWWjuCilW5b0XlBqOP7j6KTxdkFSBu0eUM3aWPL7fPqy1O9B\nYPQeuUPSa8+9djwvx5GPR4+BUQxi6QBCxt4LY80eYcwpSIob3p05wvhJUhM0wf+rq+1Offxbgc6I\npkA8cq3wvxxHJvIfAwMZxCLLL6tjyRnYec6Vw7kPX4Q5BUk+Mh9uxqFrtPey0xfnJCuensG4EQ0C\nHbuti/J/yP9yHJmIHgOJlA3GMog7Wn4dXJH8DxZ1x9sFSQSfown5Lh3OFyFh1CcQMsGL7WmYmQDZ\nxxnROhQQGRMSc2nGclG689qxsfOkk3Py/8CAIHn8ZDCYNlS8lx3ni5AE1ISZGciBI4Hi4G/+l+PI\n5/3kOOGz1dYsalnMzG2nyP/yP/5wl6yPvZ1nROeK+QVJkw6L2WZIP15F0b6xnC9CElAT12AsB447\nlAV/m0OMBcR9t9N6YZFKGcVF+eTIfsJnYhsII4W/6oX5BalT1G8ZThv50TJ8uhyd0WqLDUO/2ApU\n2RxiTNLBdelkKYokakmiB7v56UjH7cybM4url+ZVsJn4dD6zsewkECRZsO/ola2NNcpyjQ47iXn4\nwpggogFDhpnJTNRmrNkUOYyN6kV23zO3fSt8LiyrM6J5hJwsv5nbTqV3VFZ63guXe/QJ8P++gCB1\nQI4twmZUrizQroNbgHLCPTWI6tAjr50JIhownlaCaQK+RbgzkkRxd+SEE/9IbM3+fnKcUZSVPSBI\nsiA7QTajDAVOLZFu8ezOEg0h5XtKzTSJ4ulWdaYJ+BbhbCRlFldnFleTVmO6QZx1Au7evoKU4cx8\nrcZy1mHMKUgeTQXx+WIrsG9Ebh8+78sl8ifJGCRfYYnL2W938L//gmI61IOteea2U6LgOgmDg08i\nLYGi5zPSEliQMpx2WlhnIWS2b44mmFOQFONutoZ9oypYPWqszlpZPjeygXg2AWUimv2W9tqZKeBb\nxNioEAkPpLGcdQJkCEakJbBi8Sj2ScoNZySBIHncr7GZm1GwetQoIXYYZfncSJcd/zO0nSKdaYZE\n5NAzxP4L8om0BJJGkuhP/I+uXIJFKD1xYHriwIrFo5hdlM2FKGF+QfJoHavbXHZ6rCHwdJ2KUULs\nMAryuZlyVl+UacbdE0iGaRm0g5ZmxsiwisWjZowMI62l9MSBBXOH61swNURaApckDdTR3ygapPIc\n6IRhvWMsIJ9ISyAqu/FZzq6ahgs/GxsVIohoVvHFTvtZU87qPzkyTBAbq6156e4KZw+VyDwymb9O\nINISKNy71dZsjvZlD7k9vOEwv4XUKZ16fvRq3Y5pqTqd8TZAFjsRnuZzM5YJKBORS6qwrNbZeyky\nj0zmr3OJIR5gPiF9Qsr2G9MREKROlpK5+j6jV4WcZekUo2SxI/Eon1thWZ3hTECZkLPfVluzaCZJ\nHG9mtABogDGiDCDGspZAkDpHr2ABj/pr0p1llM7ao3xuokkmo4iuHCItgeQ6m8ziaqGtDZrMDdAR\niY2m+A9MBUHqHL2MD8U5v/l/7ATkB38vNfWsvmhaSJg0AvMI8BRD76ZmTkHyaO9enttPvjvYoBne\n5Ad/m3gVDuo4n48QKiyrHbj84MDlB0XZCsA8AsyNOQWJHgw2QyKRnxDaoPMrMoO/nQK+TTirPzYq\npENGXafV0IbLVgDogsgTbizML0ieDiqdjST2myEJSLiDSQyd4U3OZt7m23bBmUhLoIRHznDJ3AAd\nMVYPQGJ+QeoUnhtPgTvRcIMjOcHf5tt2wSVjo0JcCk96op6LKwHD4W4Wmee+DgMLY8VYa8Ur8vhJ\nyVNYVufSW0UaEPw/cyKcgwlF9yiy/8xqIaGbC/ufHBkmTB1FWrpFhiiMbQG8ljHRIchp+11DAILk\nGYwD2LA7WOie9pXWuhQkQ68Y7fQevcFfRxJpCZxhMf9tAvQwUJytCHDZIdSZVcF+MyRvQzr429By\nCwDsMW5cAwiSGOdAL303FJETaGfQEDsBUfA3mbptaV6F0e8OANhjUDevOQXJ01kf6WBuffcF7zTQ\nztAhdhjRnqHpuyuwJmUWV5Np3ESB0QAAuMOgXjtzCpIapBfVsm9mj3bqM64BIUpxnXmkemlehdG3\nDQUAvXDOhGmIwRwIUieIBIB9o4oSCDkbf4YOsRMQpSqw2prTO4YJMdj+GQBMA1hILKivr09ISPjw\nww/l/0ROw0j04zzEfEsnEDLNnP+MkWHuln+mJw4ENQIA+Rg0rsFggpSRkXH+/Plr167Ru4RIgXjY\n1kE6rsFMc/5PjgxzDuweGxUCzjoA8AYMIEjt7e1Wq3Xfvn3z58/fuHEjjUtIWFGVHHT3EnENJoho\nIMHpc8i7GBsVUpByp45FAgCDInIqMM7DqQwDLIwtLy+fPHkys8uJJo14cNmRiIpnvk0ZIi2BBXOH\nL91dUVhWO2NEGNhGAKAMI+ZrMIAghYWFrV27Fn8+fvx4RkaGXiXRa55QdF0huU5mcbUpN2UQBTgA\nAKAAI8Y1GECQgoODk5KS8Gc/PyoFlsi/wEOaBlFynaziaixIZG7ssVEhMO0PAICAqN8wBAaYQ1KA\nmn32RJM0+qZpECCnVYQtrslH7cmR/XQoFgAAgHaYU5BI3l63LjY2NjY2VuI7EvkXOIlhE+2UM3Pb\nt6JFo+aYQAIAQEPIcerRwtxOe0Ld4chlt3DhwpaWFvLIsmXLQkLUasDqVxbOGPmGypPoDp5WEUTI\namvOtP0YXwdqBACAM2OjQgSvXf1HixFajBDiWZM4EqS9e/c2NjaSR1588UX1guQppEkkWvSjb1D1\njJFhWcUXXXqEJXYaBQDAa8EDWQOtBuFIkI4dO6bXpd2Fozhn/taX95PjBi4/KDo4Y2SYgR44AABY\nYqzOwfxzSApwGRPBQxIEl/HQhk4XBAAAIACCJEVlh40nuBhoYKcweQQmkAAAMAcgSAi5Fxve0jQg\nJ6cwLCAFAMA0cDSHpCEqFw9Za5udJYqfZc84uU5hWS1sWAcAgJkwmCCNHz++pKSExpkjLYHOU0dk\nSJteaRpcEmkJnGEBTx0AAKYCXHYAAAAAF5hfkCSyMHT4GpGbXYj25iRNAwAAgDdgfkECAAAADAEI\nkgtwtDdXaRoAAABMDwjSDZz1hofNywEAALwHECRZGGL3XwAAAEMDgnQDcpkRXg/LYZoGAAAAE2NO\nQVKzQZ/wcw7TNAAAAJgYcwqSAsZE/xjVjRM9kKrGT5oGAAAAs2J+QZKpJR1cdrbmwrI6Mv8QV2ka\nAAAATIn5BUkmkZZAcunr0ryKjn+FoAYAAAC6gCD9CLn/vLW2CVx2AAAALAFB+hHSQlIZFgEAAAB4\nCgjSj4i8dqI/MS4MAACAtwGC1IGxrrYDBzUCAABgAAhSB8jgbwFI0wAAAMAAEwqS1dY8NipkbFSI\nf+MVT42bsVG9XOwVCxYSAAAAfQy2Y6wcIi2BBSl3IoRiY5MVbC87Niok01ZNoVwAAACAFCa0kFTy\n5Ejx1uBjXE0sAQAAANoCgiQGlhwBAADoAgiSGOfgb0jTAAAAwAAQJBeQKRsQ2EwAAABMAEFygZOF\nBIIEAABAHRAkF0RaAt9PjsOfC1KG61sYAAAAL8GEYd+aMGNk2NioELCNAAAAmAEWkltAjQAAAFgC\nggQAAABwAQgSAAAAwAUgSAAAAAAXgCABAAAAXACCBAAAAHABCBIAAADABSBIAAAAABfotjC2vr5+\nypQps2bNmj59eqdfzs3NLSgoEB3s3bv3okWL6JQOAAAAYI1ugpSRkXH+/Plr167J+fKuXbvy8/OD\ngoLIg/3796dTNAAAAEAHmApSe3v7uXPnKisrd+zYkZubK/+HJSUlw4cP37JlC72yAQAAAPrCdA6p\nvLw8KSnpqaee8kiNGhsbq6qq4uLi6BXMxMTGxupdBP2BSkBQCQghqATuYWohhYWFrV27Fn8+fvx4\nRkaGnF+dOXPG4XAMGjSIZtEAAAAAnWEqSMHBwUlJSTcu7Cf30iUlJQihfv36rVq16sSJE926dYuN\njX388cf79OlDq6AAAAAAcwyw/QQWpHnz5rW0tERGRp47d66goGDr1q1vvfXWqFGj9C4dAAAAoA3G\nECQfH585c+bMnDkzICCgvb1948aNa9asWbRoUW5ubnBwsMRvwWWMoBIQQlAJCCGoBIQQVALfUBGk\nhQsXtrS0kEeWLVsWEhLi7vvSrFmzprW1VQjy9vX1nTt3bnl5+c6dO/Py8qZNm+buh9i0AgAAAAwB\nFUHau3dvY2MjeeTFF19ULEihoaHOBxMTE3fu3Hn27Fll5wQAAAB4g4ogHTt2TMOz2e12Hx8fHx8f\n8iD21F25ckXDCwEAAAA6wnsuO6vVGhcXt2DBAtFx7I6Ljo7Wo1AAAACA9nAnSDabLS8vLy8vr7W1\nFSEUERFhsVjy8/NPnz4tfKeurm7Tpk1+fn4TJ07Ur6QAAACAlnAXZXf27NnU1FSE0OHDh3v27Onj\n45Oenp6ampqcnPzYY4/FxsZWV1dv2bLl8uXLaWlpt99+u97lBQAAALSBO0FyJikpaf369atXr968\neTNCyMfHZ8CAAevWrQPzCAAAwEz4OBwOvcsAAAAAAPzNIQEAAADeCQgSAAAAwAUgSAAAAAAXgCAB\nAAAAXGCAKDtP2bNnT05OTmlpaXBwcExMzOzZsyMiIvQuFBWU3emOHTv27dt35syZHj16DBkyZM6c\nOX379mVQWkqobO6mpqbHH3/cYrFs3LiRXiFpo6wSLly4sGfPnqKiooqKip/85CcPP/xwYmIig9JS\nQkElOByO3bt3//vf/66srOzVq1d0dPSMGTOioqLYFFgX6uvrp0yZMmvWrOnTp+tdFheYLcpu+fLl\n2dnZAQEB8fHxDQ0N5eXlAQEB77zzzq9+9Su9i6YxCu60vb3997///d69e4ODgwcPHnzx4sXKysrA\nwMBNmzaNGDGCZeG1Qn1zv/zyyx999NHgwYO3b99Otaj0UFYJ33333WOPPfb999/37du3d+/e3377\nLUIoJSUlLS2NVcG1RFkl/OEPf/jkk0+CgoKGDBlSWVl58eJFPz+/NWvWmHhJyRtvvLFx48YFCxY8\n/fTTepfFFQ4TUVRUFBMTM2HChKqqKnzk888/j4uLu/fee5uamvQtm7You9MtW7bgkaPwnY8++igm\nJuaee+65fv06i3Jrivrm3rt3b3x8/LBhw6ZMmUKzpBRRVgk//PDDuHHjhgwZ8uWXX9rtdofDUV5e\nPmLEiEGDBgnnMRDKKmHv3r0xMTEPPfRQbW2tw+Gw2+25ubkxMTF33313W1sbo6Izoa2traKiorCw\ncN68eTExMTExMRs2bNC7UK4x1RzSxx9/jBBauHBheHg4PpKUlHTfffddvHjxwIEDuhZNY5TdaWZm\npp+f32uvvRYYGIiPPPLII/fcc8+lS5fOnDnDoNjaorK5v//++8WLF8+dO9disdAtKE2UVcLu3bvP\nnz//wgsvjB49GqctHjhw4OzZs0NCQr755hs2JdcQZZXw9ddfI4Rmz57dq1cvhJCPj8+kSZOio6Ov\nXLlSXl7OpOCMKC8vT0pKeuqpp3Jzc/UuSyeYSpAOHTrk6+s7ZswY8uD48eMRQkVFRToVigoK7tTh\ncFy4cCE2NvbWW28ljw8cOBAhVFVVRa2wtFDZ3C+99FJYWNjcuXNplY8JyiohJyfH19d36tSp5MFn\nnnnm4MGD9913H6Wi0kNZJbgciLS3t3fp0qVPnz6aF1JHwsLC1t5kzpw5ehdHCvMENbS2ttbU1Nx2\n223C8B+DpyiN2OG6Q9mdtre3L1261Dl+obS0FCEUGRlJpazUUNnc27ZtO3DgwPbt2319fSmWkjKK\nK+Hw4cPDhg0LDg4+e/bssWPHzp8/Hx0dfffddxuxI1ZcCRMmTHjjjTcyMzMTEhKCgoIQQrt3766o\nqPjlL3+pePM2PgkODk5KSsKf/fy47vO5LpxHNDQ02O32nj17io7jI/X19XoUigrK7tTPz895d92D\nBw9+9dVX0dHRhtvIQ01zW63W119/PS0t7Y477qBYRPooq4SrV69ev349LCwsIyNj1apVwvEePXr8\n5S9/mTBhAr0C00DxkxAREZGZmZmampqQkHDXXXeVl5dXVFSMHj16zZo1dEsMuMc8Lju8XYWz/uMj\n+K/mQKs7/eKLL1JSUgICAv785z8bzlBQXAnt7e0vvPDCoEGDZs2aRbWEDFBWCXhby6Kiorfffvvl\nl1/+8ssvv/zyy8WLFzc3Ny9YsOD8+fOUS60xal6Hq1ev2u32urq6wsLCiooKhFBjY2NtbS21wgKd\nYB5Bwl2q8/OHjxiuw5VA/Z02NDS8+OKLzz33XFBQ0KZNm4YPH06jnFRRXAnvvPNOaWnpihUrunQx\n/MOvrBLwX2022wsvvDB9+vTQ0NDQ0NAnnngiLS2tpaUlIyODcqk1RvGT8OGHH6akpISGhmZnZx8/\nfry4uHjx4sXffvvtQw89hPf/BNhj+HdSoFu3bgihpqYm0XF8BP/VHKi80/z8/Pvvv3/79u0PPPBA\nTk6OQVcgKauEb775ZsOGDc8884zFYrl6E7vdbrfbr1692tjYSLvY2qKsEnBQGUJoypQp5PHJkycj\nhE6dOqV5Oami+HX4+OOP/f39N27c+Itf/MLPz69Hjx5PPPHEggULmpqavvjiC6plBtxhHkEKDg7u\n3r37uXPn2tvbyeNWqxUhFBYWpk+xKKDmTt9///25c+d27do1Kyvrr3/9q3Enb5VVwokTJ9rb2//6\n17+OILh48eKpU6dGjBjx2GOPMSi5hiirhF69evn7+99yyy14Jl8APww2m41WcemgrBJqa2tPnToV\nExMjRIpjcK4KkwXlGgjzBDUghOLi4g4fPnzy5Mlhw4YJB48ePYoQio+P169c2qPsTvPz81esWDFi\nxIj33nsvODiYRUFpoqASBg8enJKSIjqYnZ3drVu3Rx55xIgxZgoqwd/f/5e//OV//vOf8+fPk90x\njrc0Yp4tBZUQFBTk6+t79epV0fG6ujp0U5sB9pjHQkI3Rzdk4FB1dfXWrVt9fX0TEhL0K5f2yLlT\nm82Wl5eXl5eHnent7e0rVqzo3r27OdQIKaqEYcOGpTnRo0ePW2+9NS0t7X//9391uRE1KKgE4Vdv\nvfWW42bmMLvd/s4776CbjjtjoaASunbtGh8f/9133+Xl5Qm/cjgcGzZsQAjdddddTG8AuImpLKTk\n5ORt27YdPnx46tSpDzzwwKVLl3bt2tXY2Pi73/2uf//+epdOS+Tc6dmzZ1NTUxFChw8f7tmzZ2lp\naWVl5YABA958803nE86cOfO2225jeg+qUVAJupaXCsoqYdq0aYWFhZ988snFixdxb/75558XFxff\nddddv/71r3W8HWUoq4T09PTk5OTU1NQhQ4ZMmjSpa9euubm5x44dGzp0qBGHJubAVILk7++flZW1\nbNmyPXv2nDx5EiEUFBT0/PPPz549W++iaYyCOz1+/DhC6Lvvvvvwww+d/zp58mTDCZL3NLcEyirB\n19f3zTff/Mtf/rJ///4///nPCCGLxfLEE0/84Q9/MGLwobJKGDx48IcffvjGG28cOnToxIkTCKGu\nXbs+8cQTzz33nL+/P6OiAx0xW7ZvTGtrq9VqDQoKCgsLw6m6zIr33KkEUAlIRSVcuHDBx8fHHFE/\nyiqhqampqqoqICAgPDzcTOtDjIg5BQkAAAAwHMYzzwEAAABTAoIEAAAAcAEIEgAAAMAFIEgAAAAA\nF4AgAQAAAFwAggQAAABwAQgSAAAAwAUgSAAAAAAXgCABAAAAXACCBAAAAHABCBIAAADABSBIAAAA\nABeAIAEAAABcAIIEAAAAcAEIEgAAAMAFIEgAAAAAF/w/g1t3KAgIkIcAAAAASUVORK5CYII=\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"m = 100; n = 15;\n",
"t = (0:m-1)'/(m-1);\n",
"A = t.^0;\n",
"for j = 2:15, A(:,j)=t.^(j-1); end\n",
"b = exp(sin(4*t));\n",
"x = A\\b;\n",
"plot(t,b-A*x)\n",
"title('Residual of least squares fit')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see that the fit is quite good. As in the text, we will renormalize the vector $b$ to make the last coefficient of the solution equal to 1. (See the section below for more on this constant.) "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"b = b/2006.787453104851834;\n",
"x = A\\b; y = A*x;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As in the text, we compute all the numbers associated with the conditioning of the least squares problem, as described by Theorem 18.1. Since the residual is small, the conditioning of $x$ will be close to $\\kappa(A)$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"kappa =\n",
" 2.2718e+10\n",
"theta =\n",
" 3.7461e-06\n",
"eta =\n",
" 2.1036e+05\n"
]
}
],
"source": [
"format compact, format short e\n",
"kappa = cond(A)\n",
"theta = asin(norm(y-b)/norm(b))\n",
"eta = norm(A)*norm(x)/norm(y)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"K =\n",
" 1.0000e+00 1.0800e+05\n",
" 2.2718e+10 3.1909e+10\n"
]
}
],
"source": [
"K(1,1) = sec(theta); K(2,1) = kappa*K(1,1);\n",
"K(1,2) = K(2,1)/eta; K(2,2) = kappa + kappa^2*tan(theta)/eta"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Backslash\n",
"\n",
"The internal backslash does Householder QR, but includes column swapping to improve accuracy. Here's how it does."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x(15) = 0.999999588983230\n",
"errBS =\n",
" 4.1102e-07\n"
]
}
],
"source": [
"x = A\\b;\n",
"fprintf('x(15) = %.15f\\n',x(15))\n",
"errBS = abs(x(15)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We've lost about 9 digits of accuracy, which is in line with the conditioning of $x$ with respect to $A$ (the (2,2) element of $K$ above). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Householder QR solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, QR without column pivoting happens to do a little better."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x(15) = 0.999999919426074\n",
"errHH =\n",
" 8.0574e-08\n"
]
}
],
"source": [
"[Q,R] = qr(A,0);\n",
"x = R\\(Q'*b);\n",
"fprintf('x(15) = %.15f\\n',x(15))\n",
"errHH = abs(x(15)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modified GS QR solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Classical Gram-Schmidt is too unstable even to bother with. The modified form is better, but problems remain:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x(15) = 0.939892735084297\n",
"err =\n",
" 6.0107e-02\n"
]
}
],
"source": [
"[Q,R] = mgs(A);\n",
"x = R\\(Q'*b);\n",
"fprintf('x(15) = %.15f\\n',x(15))\n",
"err = abs(x(15)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Crikey! That's much worse than is accounted for by the problem conditioning. The culprit is our nonorthogonal $Q$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
" 3.5947e-07\n"
]
}
],
"source": [
"norm(Q'*Q-eye(15))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, we can imitate Gaussian elimination and use an augmented matrix in order to get an \"implicit Q\" result. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x(15) = 1.000000091910366\n",
"errMGS =\n",
" 9.1910e-08\n"
]
}
],
"source": [
"[Q,R] = mgs([A b]);\n",
"x = R(1:15,1:15) \\ R(1:15,16);\n",
"fprintf('x(15) = %.15f\\n',x(15))\n",
"errMGS = abs(x(15)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Normal equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This problem is a perfect storm for the normal equations: a very good fit (small residual) and a badly conditioned matrix. Squaring that condition number is catastrophic."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[\b> In pymat_eval (line 31)\n",
" In matlabserver (line 24)]\b \n",
"[\bWarning: Matrix is close to singular or badly scaled. Results may be\n",
"inaccurate. RCOND = 4.833016e-19.]\b \n",
"x(15) = -0.094071919024642\n"
]
}
],
"source": [
"x = (A'*A)\\(A'*b);\n",
"fprintf('x(15) = %.15f\\n',x(15))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SVD"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The SVD provides the gold standard of stability in least squares. It is preferred to QR when the matrix is nearly or numerically rank deficient."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x(15) = 0.999999919426927\n",
"errSVD =\n",
" 8.0573e-08\n"
]
}
],
"source": [
"[U,S,V] = svd(A,0);\n",
"x = V*(S\\(U'*b));\n",
"fprintf('x(15) = %.15f\\n',x(15))\n",
"errSVD = abs(x(15)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Extended precision\n",
"\n",
"If you compare my results to those in the text, you will find that they are different in some respects. For one thing, I use a different normalization constant for $b$. This constant is supposed to come from a higher-precision computation of the problem. \n",
"\n",
"There's some ambiguity here. Should all computations be evaluated in extended precision from the start of the fitting problem, or only once the data $A$ and $b$ have been formed in double precision? Here's the second interpretation."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
" 2006.7875278311014035999679357743\n",
" 2006.7875278311014035999679357742\n",
" 2006.7875278311014035969129213443\n",
" 2006.7875278311014035999679357743\n"
]
}
],
"source": [
"b = vpa( exp(sin(4*t)), 64 );\n",
"A = vpa(A,64);\n",
"[Q,R] = qr(A,0);\n",
"x1 = R\\ (Q'*b);\n",
"[Q,R] = mgs([A b]);\n",
"x2 = R(1:15,1:15) \\ R(1:15,16);\n",
"x3 = (A'*A)\\(A'*b);\n",
"[U,S,V] = svd(A,0);\n",
"x4 = V*(S\\(U'*b));\n",
"\n",
"[x1(end);x2(end);x3(end);x4(end)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice how the results not from the normal equations agree up to the last digit. The disadvantage of this interpretation is that it could be dependent on the precise double precision values found for the data, before extended precision is used. By putting extended precision in from the start, we get results that will hopefully be machine- and language-independent. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
" 2006.7874531048518338761038143559\n",
" 2006.7874531048518338761038143553\n",
" 2006.7874531048518338766907539159\n",
" 2006.7874531048518338761038143555\n"
]
}
],
"source": [
"t = vpa(0:m-1,64)'/vpa(m-1,64);\n",
"A = t.^0;\n",
"for j = 1:14, A=[A,t.*A(:,j)]; end\n",
"b = exp(sin(4*t));\n",
"\n",
"[Q,R] = qr(A,0);\n",
"x1 = R\\ (Q'*b);\n",
"[Q,R] = mgs([A b]);\n",
"x2 = R(1:15,1:15) \\ R(1:15,16);\n",
"x3 = (A'*A)\\(A'*b);\n",
"[U,S,V] = svd(A,0);\n",
"x4 = V*(S\\(U'*b));\n",
"\n",
"[x1(end);x2(end);x3(end);x4(end)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Unfortunately, neither of my results matches the T&B number of 2006.787453080206, which I can't explain. "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Matlab",
"language": "matlab",
"name": "matlab"
},
"language_info": {
"codemirror_mode": "octave",
"file_extension": ".m",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "matlab",
"version": "0.11.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment