{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Randomized SVD" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import numpy.linalg as la\n", "import matplotlib.pyplot as pt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Obtaining a low-rank matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's begin by making a low-rank matrix:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.0999639133084558e-13\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VdW9//H3NyNkgCSQIJCEEEAGkTFEJhVQKeCAOIJWcbo429pah2trrb1t9dpqa8UBFcHW4gQKKg5UUWQmjDLPkJBAwhTGAEnW748c+8vFBEKmnZzzeT3PeXL2Omvv/d1snu/ZZ+211zLnHCIiEjiCvA5ARERqlxK/iEiAUeIXEQkwSvwiIgFGiV9EJMAo8YuIBJjTJn4zSzKzmWa2xsxWmdnPfOVxZjbDzDb4/saWs/5oX50NZja6ug9ARETOjJ2uH7+ZNQeaO+eWmFk0sBi4ErgF2Ouce9rMHgVinXOPnLRuHJABpAHOt25P59y+aj8SERGpkNNe8TvncpxzS3zvDwJrgJbAcGCir9pESr4MTvYTYIZzbq8v2c8AhlRH4CIiUjkhZ1LZzFKA7sACoJlzLgdKvhzMLKGMVVoCmaWWs3xlp9S0aVOXkpJyJqGJiAS0xYsX73bOxVekboUTv5lFAZOBnzvnDphZhVYro6zMtiUzGwOMAUhOTiYjI6OioYmIBDwz21bRuhXq1WNmoZQk/bedc1N8xbt87f8/3AfILWPVLCCp1HIikF3WPpxz45xzac65tPj4Cn1piYhIJVSkV48BbwBrnHPPlfpoGvBDL53RwNQyVv8CGGxmsb5eP4N9ZSIi4pGKXPH3A24CBpnZMt9rGPA0cImZbQAu8S1jZmlm9jqAc24v8Htgke/1lK9MREQ8ctrunF5IS0tzauMXEak4M1vsnEurSF09uSsiEmCU+EVEAowSv4hIgPGbxF9U7Bg7cyMrsvZ7HYqISJ3mN4n/0LFC3p6/jfsnLeVgwQmvwxERqbP8JvE3bhjKC6O6k7XvKI9/uJK62FtJRKQu8JvED5CWEseDF7dj2vJs3l+c5XU4IiJ1kl8lfoC7B7Slb5sm/HbqKjbmHvQ6HBGROsfvEn9wkPH89d1oGBbMrz9a6XU4IiJ1jt8lfoBmjRpw38C2zN+8l7mbdnsdjohIneKXiR/ghvOSadYonL/O2KAbvSIipfht4m8QGsx9A9uycOteZm/UVb+IyA/8NvEDXNcriRaNG/DcjPW66hcR8fHrxB8eEsx9g9qxdPt+vlmf53U4IiJ1gl8nfoBreiaSFNeQJ6etYv+R416HIyLiOb9P/GEhQfz1+u7k7C/g/klLKSwq9jokERFP+X3iB+jZKpb/ubIz323YzdOfrfU6HBERT4WcroKZjQcuA3Kdc519Ze8C7X1VYoD9zrluZay7FTgIFAGFFZ0dpiZc1yuJ1TkHeH32FpLiIripdyuCgsyrcEREPFORK/4JwJDSBc65651z3XzJfjIw5RTrD/TV9Szp/+DxSztyfrum/HbaKi5/cTbfrMtVbx8RCTinTfzOuVlAmROkm5kB1wGTqjmuGhEaHMSEW9N57rqu5B89wS1vLuKBd5Yp+YtIQKlqG//5wC7n3IZyPnfAl2a22MzGVHFf1SI4yLiqRyJf/3IAYy5I5ePl2Xy3QQ94iUjgqGriH8Wpr/b7Oed6AEOBe83sgvIqmtkYM8sws4y8vJrvcx8WEsQvB59NYmxDnvl8LcXFuuoXkcBQ6cRvZiHAVcC75dVxzmX7/uYCHwLpp6g7zjmX5pxLi4+Pr2xYZyQ8JJiHBrdnVfYBPl6RXSv7FBHxWlWu+C8G1jrnypzxxMwizSz6h/fAYKDOjZN8RdcWdGzeiD9/uY7jherjLyL+77SJ38wmAfOA9maWZWa3+z4ayUnNPGbWwsym+xabAbPNbDmwEPjUOfd59YVePYKCjEeGtCdz71H+tWCb1+GIiNS40/bjd86NKqf8ljLKsoFhvvebga5VjK9WXHh2PH1Sm/DM5+socnBL3xSC1cdfRPxUQDy5ezpmxl9HdqN3ahy//2Q1V700hzU5B7wOS0SkRijx+zRr1IDxt/TibyO7kbnvKJe+8B0PvruMzXmHvA5NRKRanbapJ5CYGcO7teSCdvG88u0mJs7bytRlO7imZyJPXnEOEWH65xKR+k9X/GWIjQzjsWEd+e7hQdzarzUfLM7ixtcXsO+whnUWkfpPif8U4qPD+c1lnXjpxh6syj7Ata/OI3v/Ua/DEhGpEiX+ChjSuTlv3ZbOrvwCrn55Ltv2HPY6JBGRSlPir6DeqU14984+FJwo4obXFrBDV/4iUk8p8Z+BTi0a8Y/bz+NAwQlufG0+uQcKvA5JROSMKfGfoc4tGzPxtnTyDh5j1Gvz+XLVTg31ICL1ihJ/JfRIjuWNW3qRf7SQMf9YzHl//DdPTF3JXvX6EZF6QIm/knqnNmHeY4N485Ze9GvblHcWZjJ87GzW7TzodWgiIqekxF8FocFBDOyQwIs39OC9u/pQcKKYq16aw1drdnkdmohIuZT4q0m3pBim3deP1Pgo7ngrg2c+X8uxwiKvwxIR+REl/mrUvHFD3ruzD9enJfHyN5u44u9zWLkj3+uwRET+DyX+atYwLJinr+7Cm7f0Yt+R41w5dg5frtrpdVgiIv+hxF9DBnZI4MsHL6Bds2ienLaKI8cLvQ5JRARQ4q9RMRFhPDX8HLLzC3jlm01ehyMiAlRs6sXxZpZrZitLlT1pZjvMbJnvNaycdYeY2Toz22hmj1Zn4PVFr5Q4rujagldnbSZz7xGvwxERqdAV/wRgSBnlzzvnuvle00/+0MyCgbHAUKATMMrMOlUl2PrqsWEdCDLjj9PXeB2KiMjpE79zbhawtxLbTgc2Ouc2O+eOA+8AwyuxnXqveeOG3DuwDZ+t3MmEOVs4elzdPEXEO1Vp47/PzFb4moJiy/i8JZBZajnLVxaQ7jg/lW5JMTz58WrO++O/eXLaKpZu30dRsfM6NBEJMJWdS/Bl4PeA8/39C3DbSXWsjPXKzXJmNgYYA5CcnFzJsOquBqHBfHhPXxZu2cvbC7bzrwXbmTB3KzERofRv2/Q/XwwiIjWtUonfOfefMQnM7DXgkzKqZQFJpZYTgexTbHMcMA4gLS3NLy+DzYzzUptwXmoT9h85zqwNu5m1Po+v1+YyZ+NuPv/5BTRr1MDrMEXEz1WqqcfMmpdaHAGsLKPaIqCdmbU2szBgJDCtMvvzRzERYVzRtQV/vrYr793Zh6Mnivjle8spVtOPiNSwinTnnATMA9qbWZaZ3Q78r5l9b2YrgIHAg766LcxsOoBzrhC4D/gCWAO855xbVUPHUa+1TYjit5efw+yNu3l99mavwxERP2fO1b0rzLS0NJeRkeF1GLXKOcdd/1zM12tz+fCefnRu2djrkESkHjGzxc65tIrU1ZO7dYSZ8fRVXWgSGc4Nr81nypIs6uKXsojUf0r8dUhsZBjvjOlNu2bR/OK95dwxMYNdmtdXRKqZEn8dk9I0kvfu7MOvL+3I7I27ufrlueQfOeF1WCLiR5T466DgIOOO81OZNKY3uw4U8NAHy9XsIyLVRom/DuuRHMujQzsyY/Uu3pi9xetwRMRPKPHXcbf1S2Fwp2Y8/dlalmzf53U4IuIHKjtkg9QSM+PZa7py6d+/4+qX55IUG0H7s6LpkRzLyF5JxEaGeR2iiNQz6sdfT2TuPcKUJTtYv+sga3ceYFPeYRqGBjMyPYk7zk+lZUxDr0MUEQ+dST9+Jf56at3Og7z67SamLs+mQUgQH9/fn9T4KK/DEhGP6AGuAND+rGieu74b//7FhYSGBPHgu8s4UVTsdVgiUg8o8ddzrZtG8scR57I8K58XvtrgdTgiUg8o8fuBYec255qeiYyduZFFWyszWZqIBBL16vETT15xDgu37OXufy6mW1IMUeEhJDRqwH+dn0p8dLjX4YlIHaIrfj8RFR7CSzf2oMNZjcjJL2Bp5n7enLOFoX+bxcy1uV6HJyJ1iHr1+LH1uw7ywKSlrN15kNF9WvH4pZ0IC9F3vYg/Uq8eAeDsZtF8dG8/bu/fmonztvHo5BUa80dE1Mbv7xqEBvObyzrRqEEoz/97PYlxEfzikrO9DktEPHTaxG9m44HLgFznXGdf2bPA5cBxYBNwq3NufxnrbgUOAkVAYUV/hkj1e+CitmTuO8ILX20gMbYh16UleR2SiHikIk09E4AhJ5XNADo757oA64HHTrH+QOdcNyV9b5kZf7rqXPq3bcp/T/mem8cv5E/T1/DR0h0cOV7odXgiUotOe8XvnJtlZiknlX1ZanE+cE31hiU1ITQ4iJd+2oNnPlvLssz9vDlnD8eLikmMbcgfRpzLhWfHex2iiNSC6mjjvw14t5zPHPClmTngVefcuGrYn1RBowah/GHEuQCcKCpm0Za9/GbqSkaPX8iV3Vrwuys60zgi1OMoRaQmValXj5k9DhQCb5dTpZ9zrgcwFLjXzC44xbbGmFmGmWXk5eVVJSypoNDgIPq2bcr0n53PAxe149Pvc/j5u0vV80fEz1U68ZvZaEpu+t7oyskUzrls399c4EMgvbztOefGOefSnHNp8fFqcqhN4SHB/OKSs/nvYR2ZuS6Pdxdleh2SiNSgSiV+MxsCPAJc4Zw7Uk6dSDOL/uE9MBhYWdlApeaN7pNCn9Qm/P6T1WTuLfO0iogfOG3iN7NJwDygvZllmdntwItANDDDzJaZ2Su+ui3MbLpv1WbAbDNbDiwEPnXOfV4jRyHVIijIePbaLpgZD72/nOJiNfmI+CMN2SA/8t6iTB6evIIR3Vvy4MVnk9wkwuuQROQ0zmTIBj25Kz9ybVoim3Yf4s3ZW5m2PJsrurbg7gFtOLtZtNehiUg10BW/lGvXgQJem7WZtxds5+iJIi7u2Ix7BrahR3Ks16GJyEk0565Uq32HjzNh7lYmzN1K/tET9GwVy/VpSVzapTmR4frRKFIXKPFLjTh8rJB3FmXyrwXb2JR3mMiwYC7r0oJr0xLp2SoWM/M6RJGApcQvNco5x5Lt+3h3USafrsjh8PEiUptGMio9mVHnJROlXwEitU6JX2rN4WOFTP8+h/cyMlm0dR+NG4Yyum8Kt/ZNITYyzOvwRAKGEr94Ylnmfl6auZEvV++iaVQY74zpQ9uEKK/DEgkImoFLPNEtKYZxN6fxyf39AeOnry/QE8AidZASv1S7zi0b88870jl6oogbXp/PzvwCr0MSkVKU+KVGdDirEW/dls6+wycY9dp8NuUd8jokEfFR4pca0zUphom39SL/6AmGvziHL1bt9DokEUGJX2pYz1ZxfHx/f9rER3LnPxbzp+lryD2gph8RL6lXj9SKghNFPDltFe8sysQMerduwhXdWjCsc3PN+CVSDdSdU+qsjbmHmLY8m4+XZ7Nl92HCgoMY1CGBq3smcnHHBD39K1JJSvxS5znn+H5HPh8tzWba8mx2HzrGo0M7cNeFbbwOTaReUj9+qfPMjC6JMTxxeSfmPzaIS7s055nP1/LNulyvQxPxe0r84rmQ4CCevaYLHc5qxP2TlrJl92GvQxLxaxVK/GY23sxyzWxlqbI4M5thZht8f8scpN3MRvvqbPBN0C7yIxFhIYy7qSchQcZ/vZVB9v6jXock4rcqesU/ARhyUtmjwFfOuXbAV77l/8PM4oDfAucB6cBvy/uCEEmKi2DsjT3YvvcIA579hiemrmSXun6KVLsKJX7n3Cxg70nFw4GJvvcTgSvLWPUnwAzn3F7n3D5gBj/+AhH5j75tmjLzoQFc3TORfy3Yzvn/O5P7Jy3l36t3cbyw2OvwRPxCVQZOb+acywFwzuWYWUIZdVoCmaWWs3xlIuVqGdOQP111LvcMaMNr323mY1/3z5iIUB4b2oHreyV7HaJIvVbTN3fL6pRdZv9RMxtjZhlmlpGXl1fDYUl9kBQXwVPDO7Pw8YsZf0saHc9qxCOTv+fJaas4UaSrf5HKqkri32VmzQF8f8vqh5cFJJVaTgSyy9qYc26ccy7NOZcWHx9fhbDE34QGBzGoQzP+cXs6/3V+aybM3cro8QvZfeiY16GJ1EtVSfzTgB966YwGppZR5wtgsJnF+m7qDvaViZyxkOAgHr+0E3+5tisZW/cx4NlvGDtzI0ePF3kdmki9UtHunJOAeUB7M8sys9uBp4FLzGwDcIlvGTNLM7PXAZxze4HfA4t8r6d8ZSKVdnXPRKb/rD992jTh2S/WMeDPM5m1Xs2DIhWlIRukXlu4ZS+//uh7cvYX8MkD/WnVJNLrkEQ8oSEbJGCkt45j/C29CAoy7nl7CQUn1OwjcjpK/FLvJcZG8Nx1XVmVfYDff7La63BE6ryq9OMXqTMu6tiMOy9M5dVvN2MGZzeLJiG6AT1axZAQ3cDr8ETqFCV+8RsPDW7PlrzDvLMwk8LikntXMRGhvPLTnvRObeJxdCJ1h27uit8pLnbsOXycrXsO88jkFWTuPcIfRpzLdWlJp19ZpJ7SzV0JaEFBRnx0OL1S4vjw7n6c17oJD3+wgsemfK8hn0VQ4hc/1zgilDdv7cWt/VJ4PyOTgX/+hp++voDPvs/RoG8SsNTUIwEj90AB7y7K5F8Lt5OTX0BcZBhXdmvJqPQk2jWL9jo8kSrRnLsip1BYVMx3G3bz/uJMZqzeBcAfrjyX63rpHoDUX2eS+NWrRwJOSHAQAzskMLBDAnsOHePn7y7j4ckr2Jh3iEeGdCA4qKxBZUX8h9r4JaA1iQrnzVt6cXOfVoybtZmbxy9g8uIscvI19aP4L13xS8ALCQ7iqeGdadcsmudnrOeXG5cD0CY+kku7tOCKri1omxDlcZQi1Udt/CKlFBc71u48yNxNu/l6bS7zNu/BOeiS2Jjnr+9Gm3h9AUjdpJu7ItUk90ABn6zI4aVvNlLsYOKt6Zyb2NjrsER+RA9wiVSThEYNuK1/a96/qy8NQ4MZ9dp85m3a43VYIlWixC9SAa2bRjL57r60iGnA6PELuW3CIl7/bjOrsw9w6FghdfGXs0h5dHNXpILOatyA9+7sw/Mz1jNrQ8k9gB80CA2iWaMGPHn5OQzskOBhlCKnV+k2fjNrD7xbqigVeMI599dSdQZQMhfvFl/RFOfcU6fbttr4pT7Ysf8oGVv3sjO/gD2Hj/P12lx25hcw9b5+ugksta7Wb+6aWTCwAzjPObetVPkA4CHn3GVnsj0lfqmPduw/yuV/n01cZBgf3duPqHD9oJba48XN3YuATaWTvkigaRnTkBdHdWdz3iF+9f5ytftLnVVdiX8kMKmcz/qY2XIz+8zMzqmm/YnUSX3bNuXRoR34bOVObnpjoUYBlTqpyk09ZhYGZAPnOOd2nfRZI6DYOXfIzIYBf3POtStnO2OAMQDJyck9t23Tjwepn5xzvPLtZt6at5Wc/AKaRIZxSadm9GnThD5tmmgqSKkRtdrGb2bDgXudc4MrUHcrkOac232qemrjF39QVOyYtSGP9zMy+W7Dbg4WFAIwqEMCTw0/h8TYCI8jFH9S26NzjqKcZh4zOwvY5ZxzZpZOSdOSnn6RgBAcZAxsn8DA9gkUFTtWZeczc20er87axCXPzeKXg8/m1n6tNRqo1LoqJX4ziwAuAe4sVXYXgHPuFeAa4G4zKwSOAiOd7nhJAAoOMrokxtAlMYZr0hL5zUcr+Z9P1zBn425eH91LyV9qlcbqEfGAc44Jc7fyu49Xc+cFqTw2rKPXIUk9p4lYROo4M+PWfq3ZnHeYV2dtpmPzRlzZvaXXYUmA0Fg9Ih564vJOpLeO45HJK1iWud/rcCRAKPGLeCg0OIiXbuxB06hwRrw0hxtem897izLJO3hMD4BJjVFTj4jHmkaFM+WevkxauJ2py7J5ePIKAMJDgmgR05COzaO584I2dE2K8ThS8Re6uStShzjnWJ6Vz5Jt+8jJP0p2fgFzNu5m/5ETDGgfzwMXtaNHcqzXYUodpJu7IvWUmdEtKYZupa7uDx0r5K15W3lt1mauemku/do24d6BbemT2gQzdQOVM6crfpF64vCxQt5esI3XvttC3sFjpKfEMfbGHsRHh3sdmtQBmnpRxA9Fhocw5oI2fPfwQJ4afg4rduxn1GvzyT1Y4HVoUs8o8YvUMw1Cg7m5TwoTbk1nx76jjBo3n9wDSv5ScUr8IvVU79QmTLwtnZz8Aq57dR7vZ2Ry+Fih12FJPaDEL1KPpbeO463b0gH41QcrSPuff/PL95az+9AxjyOTuky9ekTqubSUOGY+NIDF2/YxeUkWU5bsYOHWPbx5Sy/aJkR7HZ7UQbriF/EDZkZaShx/uqoL797Zh6PHixnx0lzmbDzl1BcSoNSdU8QPZe07wu0TMtiQe5B2CdG0TYiiXbMoru+VRPPGDb0OT2qAunOKBLjE2Ag+uLsP9w5sS2JsQ1Zl5/PCVxsY/NwsJi3crnGAApyu+EUCxPY9R3hk8grmbd5D3zZNeHRoB7okavwff1Grc+7WBCV+kZrhnGPSwkz+NH0NB48Vck6LRoxKT2ZE95ZEhquvR31W25OtbwUOAkVA4ck7tpLBRP4GDAOOALc455acaptK/CI160DBCaYuy+ZfC7azJucACdHhPDykA1d1b0mQpoGsl7xo4x/onOtWzk6HAu18rzHAy9W0TxGppEYNQrmpdyumP9Cf9+7sQ/OYhjz0/nKufGkOK3fkex2e1LDauLk7HHjLlZgPxJhZ81rYr4ichpmR3jqOD+/uy/PXd2XXgQKueWUuM1bv8jo0qUHVkfgd8KWZLTazMWV83hLILLWc5SsTkToiKMgY0T2RT+4/n/bNornzHxlMnLvV67CkhlRH4u/nnOtBSZPOvWZ2wUmfl9Vg+KMbC2Y2xswyzCwjLy+vGsISkTMVHx3OpDG9GdShGb+dtoonp63ieGGx12FJNaty4nfOZfv+5gIfAuknVckCkkotJwLZZWxnnHMuzTmXFh8fX9WwRKSSIsJCePWmntzWrzUT5m5l5Lh55OQf9TosqUZVSvxmFmlm0T+8BwYDK0+qNg242Ur0BvKdczlV2a+I1KzgIOOJyzvx4g3dWbfzIJe+MJtv1+uXuL+o6hV/M2C2mS0HFgKfOuc+N7O7zOwuX53pwGZgI/AacE8V9ykiteSyLi2Ydn9/4qPCGT1+IU9OW8XR40VehyVVpAe4ROS0Ck4U8czna3lzzlZS4yN5dEgH2iRE0TKmIQ1Cg70OT9CTuyJSQ2Zv2M2vPlhOTv7/n/ErNT6Sm3q34tq0JKL09K9nlPhFpMYcPlbIyh35ZOcfZce+o3y9Npcl2/cTFR7Cjb2TefDis/UrwANnkvj19SwiZyQyPITzUpv8Z/m+Qe1Ylrmf8bO38Oq3m5m9YTdjb+hBStNID6OUU9GwzCJSZd2SYnhhVHfeGJ1G1r6jXP732Xy8PFvDP9dRSvwiUm0u6tiMTx/oT2pCFPdPWsqIl+Yyc12uvgDqGCV+EalWibERfHBXH/444lzyDh7j1jcXMeKluWzfc8Tr0MRHiV9Eql1ocBA3nJfMzIcG8PRV57J1z2GuGDubuZs0B3BdoMQvIjUmLCSIkenJTL23H02jwrnpjYW8NW8rxcVq+vGSEr+I1LhWTSL58J6+XHh2PE9MXcXFz33LhDlbOFhwwuvQApL68YtIrSkqdnyyIpsJc7eydPt+IsOCubRLc67qkUh6Spxm/6oCPcAlInXe8sz9/HP+NqZ/n8Ph40W0jGnI+e2a0qNVLD2SY0ltGqkvgjOgxC8i9caR44V8uWoXHy/PJmPbPvKPljT/RDcIoUtiY7okxnB1j0TaJkR5HGndpsQvIvVScbFj8+5DLN62j+VZ+azI2s/anIMA3NI3hQcubkejBqEeR1k3KfGLiN/Yc+gYf/5yHe8syqRJZDiPDGnP1T0S1Qx0kjNJ/OrVIyJ1WpOocP50VRem3duf5LiG/OqDFVzzylxW7sj3OrR6S1f8IlJvFBc7Ji/J4unP1rLvyHGGdD6LoZ2bM7BDQsAPCa3ROUXELwUFGdemJTH4nLN48esNfLh0B9O/30lYSBCXnducXww+m8TYCK/DrPMqfcVvZknAW8BZQDEwzjn3t5PqDACmAlt8RVOcc0+dbtu64heRiigqdizZvo9PlmfzzqJMHHBr3xTuHtCGmIgwr8OrVbVyc9fMmgPNnXNLfBOuLwaudM6tLlVnAPCQc+6yM9m2Er+InKns/Uf5y5frmbI0iyAzuiXF0L9tUy7v2py2CdFeh1fjauXmrnMuxzm3xPf+ILAGaFnZ7YmIVEWLmIb85bqufPaz87nrwlQKi4p54esNDH9xDhtzD3odXp1SLTd3zSwFmAV0ds4dKFU+AJgMZAHZlFz9rypnG2OAMQDJyck9t23bVuW4RCSwZe49wpVj5xAbGcbUe/sR6cc3gGu1O6eZRVGS3H9eOun7LAFaOee6An8HPipvO865cc65NOdcWnx8fFXDEhEhKS6Cv4/qzua8QzwyeYUmhPGpUuI3s1BKkv7bzrkpJ3/unDvgnDvkez8dCDWzplXZp4jImejbtikP/aQ9n6zI4c05W70Op06odOI3MwPeANY4554rp85ZvnqYWbpvf3squ08Rkcq4+8I2XNKpGX+YvoY5GzUZTFWu+PsBNwGDzGyZ7zXMzO4ys7t8da4BVprZcuAFYKTTby0RqWVmxnPXdaVNfCT3vL2EzXmHvA7JU3pyV0QCRubeIwwfO4eYhqF8eE8/Gkf4z4BvGqtHRKQMSXERvPLTnmTuO8LtExcxd+PugJwGUolfRAJKeus4nr2mK2t3HuSG1xfQ75mvefqztazbGTh9/dXUIyIB6ejxIv69ZhcfLt3Bt+vzKCp2dGreiGvTErm5TwrB9WzYZw3SJiJyGg3Dgrm8awsu79qC3YeO8emKHKYs3cHvPl7N4WOF3Deondch1hg19YhIwGsaFc7ovil8dE9fhndrwXMz1rNwy16vw6oxSvwiIj5mxh9GnEtyXAQPTFrK3sPHvQ6pRijxi4iUEhUewos39GDv4eP84r1l5B4s8Dqkaqc2fhGRk3Ru2ZhfX9aRJ6auIv0PX9GicQO6J8feZC1hAAAHAklEQVRyx/mt6Z4c63V4VabELyJShpv7pNAlMYaMrXtZnpXP3I27+fT7HK7pmcjDQ9qTEN3A6xArTYlfRKQc3ZJi6JYUA8ChY4W8+PVG3pi9mc9X7mRUehIj05NpEx/lcZRnTv34RUTOwJbdh/nzl+v4YuVOCosd57WO47+HdaSr7wvCKxqyQUSkhrRuGsnYG3ow97FBPDykPdv2HOG6V+fx8fJsr0OrMCV+EZFKSIhuwD0D2vLpA/3pktiY+yct5fkZ6+vF2D9q6hERqaJjhUX8+sOVvL84i8iwYJLiImjVJIKLOjbj6h6JtTL8w5k09Sjxi4hUA+ccH6/IYcm2fWTuPcKG3ENs33uEdglRPDKkAxd1TMA3L1WNUOIXEfGYc47PV+7k2S/WsXn3Ydo3i+bSLs0Zdu5ZtE2Irvb9KfGLiNQRJ4qKmbw4i8lLssjYtg/noEdyDI9f2pGereKqbT+1lvjNbAjwNyAYeN059/RJn4cDbwE9KZlr93rn3NbTbVeJX0T80a4DBXyyIodXv91E7sFjDDv3LB7+SQdSmkZWedu1kvjNLBhYD1wCZAGLgFHOudWl6twDdHHO3WVmI4ERzrnrT7dtJX4R8WdHjhcybtZmXv12M8cKixjauTljLkit0rMAtdWPPx3Y6Jzb7Jw7DrwDDD+pznBgou/9B8BFVpN3N0RE6oGIsBB+fvHZfPurAYy5oA2zNuQxfOwcrn91HgUnimp8/1UZsqElkFlqOQs4r7w6zrlCM8sHmgC7T96YmY0BxgAkJydXISwRkfohoVEDHh3agfsGteWdhdvZmHuIBqHBNb7fqiT+sq7cT243qkidkkLnxgHjoKSppwpxiYjUK1HhIdxxfmqt7a8qTT1ZQFKp5UTg5GeW/1PHzEKAxoD/TmsjIlIPVCXxLwLamVlrMwsDRgLTTqozDRjte38N8LWri/1HRUQCSKWbenxt9vcBX1DSnXO8c26VmT0FZDjnpgFvAP8ws42UXOmPrI6gRUSk8qo0Hr9zbjow/aSyJ0q9LwCurco+RESkeml0ThGRAKPELyISYJT4RUQCjBK/iEiAqZOjc5pZHrCtkqs3pYwng/1cIB4zBOZxB+IxQ2Ae95kecyvnXHxFKtbJxF8VZpZR0YGK/EUgHjME5nEH4jFDYB53TR6zmnpERAKMEr+ISIDxx8Q/zusAPBCIxwyBedyBeMwQmMddY8fsd238IiJyav54xS8iIqfgN4nfzIaY2Toz22hmj3odT00xsyQzm2lma8xslZn9zFceZ2YzzGyD72+s17FWNzMLNrOlZvaJb7m1mS3wHfO7vlFi/YqZxZjZB2a21nfO+/j7uTazB33/t1ea2SQza+CP59rMxptZrpmtLFVW5rm1Ei/48tsKM+tRlX37ReL3zf87FhgKdAJGmVknb6OqMYXAL51zHYHewL2+Y30U+Mo51w74yrfsb34GrCm1/AzwvO+Y9wG3exJVzfob8LlzrgPQlZLj99tzbWYtgQeANOdcZ0pG/h2Jf57rCcCQk8rKO7dDgXa+1xjg5ars2C8SPxWb/9cvOOdynHNLfO8PUpIIWvJ/5zeeCFzpTYQ1w8wSgUuB133LBgyiZC5n8M9jbgRcQMnw5jjnjjvn9uPn55qSUYMb+iZvigBy8MNz7ZybxY8npirv3A4H3nIl5gMxZta8svv2l8Rf1vy/LT2KpdaYWQrQHVgANHPO5UDJlwOQ4F1kNeKvwMNAsW+5CbDfOVfoW/bHc54K5AFv+pq4XjezSPz4XDvndgB/BrZTkvDzgcX4/7n+QXnntlpznL8k/grP7esvzCwKmAz83Dl3wOt4apKZXQbkOucWly4uo6q/nfMQoAfwsnOuO3AYP2rWKYuvTXs40BpoAURS0sxxMn8716dTrf/f/SXxV2T+X79hZqGUJP23nXNTfMW7fvjp5/ub61V8NaAfcIWZbaWkGW8QJb8AYnzNAeCf5zwLyHLOLfAtf0DJF4E/n+uLgS3OuTzn3AlgCtAX/z/XPyjv3FZrjvOXxF+R+X/9gq9t+w1gjXPuuVIflZ7feDQwtbZjqynOucecc4nOuRRKzu3XzrkbgZmUzOUMfnbMAM65nUCmmbX3FV0ErMaPzzUlTTy9zSzC93/9h2P263NdSnnndhpws693T28g/4cmoUpxzvnFCxgGrAc2AY97HU8NHmd/Sn7irQCW+V7DKGnz/grY4Psb53WsNXT8A4BPfO9TgYXARuB9INzr+GrgeLsBGb7z/REQ6+/nGvgdsBZYCfwDCPfHcw1MouQ+xglKruhvL+/cUtLUM9aX376npNdTpfetJ3dFRAKMvzT1iIhIBSnxi4gEGCV+EZEAo8QvIhJglPhFRAKMEr+ISIBR4hcRCTBK/CIiAeb/AfP6kltVyn+vAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "n = 100\n", "A0 = np.random.randn(n, n)\n", "U0, sigma0, VT0 = la.svd(A0)\n", "print(la.norm((U0*sigma0).dot(VT0) - A0))\n", "\n", "pt.plot(sigma0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A randomly drawn Gaussian matrix: Emphatically not low-rank. Let's swap out the singular values with cleanly exponentially decaying ones:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD+CAYAAAA6c3LAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4lfX9//HnO4swZA+ZMmQY2YQNiVpGEARUUNBSRBRBEUw6xNa22tpabRuWiIIiOJEismRbJWxIUGbYMiIriIICAoHP74+c/ppvSiSQcSfnvB7XlUvuD2e87uvG88o59513zDmHiIhIkNcBRESkYFAhiIgIoEIQEREfFYKIiAAqBBER8VEhiIgIoEIQEREfFYKIiAAeF4KZ9TazyWY2x8y6eJlFRCTQ5XohmNkUMztuZlszrceY2U4z22NmowCcc7Odc48CDwH353YWERHJvrx4hzAViMm4YGbBwASgGxAB9DeziAw3edb39yIi4pGQ3H5A51yCmdXMtNwK2OOc2wdgZtOBXmaWDPwNWOic23i1xy5fvryrWTPzQ4uIyE9JSko64ZyrcLXb5XohZKEqcCjDdgrQGngS6ASUMrObnXOvZb6jmQ0BhgDUqFGDxMTEfIgrIuI/zOxAdm6XX4VgV1hzzrlxwLifuqNzbhIwCSAyMlKjWUVE8kh+XWWUAlTPsF0NOJxPzy0iItmQX4WwAahrZrXMLAzoB8zNp+cWEZFsyIvLTj8A1gD1zSzFzAY759KA4cBiIBmY4ZzbltvPLSIi1y8vrjLqn8X6AmBBbj+fiIjkDo2uEBERQIUgIiI+AVEIFy9d5vl52zj83Tmvo4iIFFgBUQjbDp9m+vpDdBmdwLtrD3D5sn6cQUQks4AohKbVS7P4qSiaVC/Fs7O30m/yWr46ccbrWCIiBUpAFAJAjXLFeHdwa16+tzHJR04TMyaB15fvJe3SZa+jiYgUCAFTCABmxn0tq7MsLproehV4ceEO7n51NclHTnsdTUTEcwFVCP9RqWQ4rw9owYQHmnPk1DnuGr+Sfy7Zyfm0S15HExHxTEAWAqS/W+jeuDJLY6Pp2bQK4/+9h+7jVpJ04Fuvo4mIeCJgC+E/yhQPI/6+pkwd1JKz59Po89pqnp+3jbMX0ryOJiKSrwK+EP7jtvoVWRIXzYA2N/HWqv10GZ3Ayt0nvI4lIpJvVAgZlCgSwp96NWTGY20JCw7i52+u4zczN3Hq7EWvo4mI5DkVwhW0qlWWBSM7Muy2Ony08Ws6jV7Ooq1HvY4lIpKnVAhZCA8N5umYBsx5oj0VShRh6LtJPPHeRlK/P+91NBGRPKFCuIqGVUsxZ3h7ft21PkuTj9EpfjkfJaXgnMZfiIh/USFkQ2hwEE/cfjMLRnSkbsUS/PJfmxj41gZSvj3rdTQRkVyjQrgGN1cswYzH2vJ8z1tJ3H+SrqMTeHvNfg3LExG/oEK4RkFBxsB2NVkSG0Xzm8rwhznbuH/SGvam/uB1NBGRHPG8EMysuJklmVkPr7Nci2plivH2w634R98m7Dr2A93GruDVz/doWJ6IFFrXXQhmNsXMjpvZ1kzrMWa208z2mNmobDzU08CM683hJTOjT4tqLI2L4o76FXl50U56v7qKbYdPeR1NROSa5eQdwlQgJuOCmQUDE4BuQATQ38wizKyRmc3P9FXRzDoB24FjOcjhuYo3hPPagBZMfLA5R0+dp+crq/j74h38eFHD8kSk8Ai53js65xLMrGam5VbAHufcPgAzmw70cs69CPzPR0JmdjtQnPTyOGdmC5xzhfYzl26NKtO2Tjle+CSZCZ/tZeHWo7x8b2Mia5b1OpqIyFXl9jmEqsChDNspvrUrcs79zjn3FPA+MPlKZWBmQ8ws0cwSU1NTczlu7itdLIx/9G3C2w+34vzFy/R9fQ3Pzd3GmfMaliciBVtuF4JdYe2q12Q656Y65+Zn8XeTnHORzrnIChUq5DhgfomqV4ElsVEMbFuTaWvSh+Ut31XwC01EAlduF0IKUD3DdjXgcC4/R6FRvEgIz/W8lX891pbw0CAGTlnPL2ds4ruzF7yOJiLyP3K7EDYAdc2slpmFAf2Aubn8HIVOZM2yfDKiI0/cXofZX35Np/gEFm454nUsEZH/IyeXnX4ArAHqm1mKmQ12zqUBw4HFQDIwwzm3LXeiFm7hocH8umsD5g5vT6WSRRj23kaGvpPE8dM/eh1NRAQAK0xD2iIjI11iYqLXMXIs7dJlJq/4itHLdhEeEsTve0TQp0U1zK50CkZEJGfMLMk5F3m123n+k8qBKCQ4iGG31WHRyI40uLEkv565mV9MWc+hkxqWJyLeUSF4qHaFEkwf0oY/927IxgPf0nVMAm+t+krD8kTEEyoEjwUFGQPa3MSSuGha1izL8/O20/f1New5rmF5IpK/VAgFRNXSRZk6qCXx9zVhb+oP3Dl2BRM+28NFDcsTkXyiQihAzIx7mldjaWw0nW+txN8X76TnK6vY+rWG5YlI3lMhFEAVbijChAea8/qAFnzzw3l6TVjFS4s0LE9E8pYKoQDreuuNLI2Npk/zakz8fC93jl3B+q9Oeh1LRPyUCqGAK1UslJf6NObdwa25ePky972+ht/P3soPGpYnIrlMhVBIdKhbnsVPRfFw+1q8u+4AXeKX89nO417HEhE/okIoRIqFhfCHuyL4aFg7ihcJYdBbG4j78Eu+PaNheSKScyqEQqh5jTLMH9GBEXfczNxNh+kUv5z5mw9TmMaQiEjBo0IopIqEBBPXpT7znuxA1TJFGf7+Fzz2ThLHNCxPRK6TCqGQu6VySWYNa8dv72zA8l2pdIpfzocbDurdgohcMxWCHwgJDmJIVB0WPxVFROWSPP3RFh58Yx0Hv9GwPBHJPhWCH6lZvjgfPNqGv9zdkM0pp+g6JoE3VuzjkobliUg2qBD8TFCQ8WDrm1gaF0XbOuV44ZNk7p24ml3Hvvc6mogUcCoEP1W5VFHeHBjJ2H5NOfDNGbqPW8HYZbu5kKZheSJyZflWCGZW28zeNLOZGdaKm9k0M5tsZg/mV5ZAYWb0alqVZXHRxDSszOhlu+j5yko2HfrO62giUgBlqxDMbIqZHTezrZnWY8xsp5ntMbNRP/UYzrl9zrnBmZbvAWY65x4Fel5Tcsm2ciWKML5/Myb/IpJvz17g7ldX8dcFyZy7oGF5IvJf2X2HMBWIybhgZsHABKAbEAH0N7MIM2tkZvMzfVXM4nGrAYd8f9arUx7rHFGJpXHR3N+yOpMS9tFtbAJr9n7jdSwRKSCyVQjOuQQg85jNVsAe33f+F4DpQC/n3BbnXI9MX1kN3UkhvRSynUVypmR4KC/e05j3H2nNZQf9J6/ltx9v4fSPF72OJiIey8mLcFX++909pL+4V83qxmZWzsxeA5qZ2TO+5VnAvWY2EZiXxf2GmFmimSWmpqbmIK5k1O7m9GF5j3SoxfT1B+kSn8Cnyce8jiUiHgrJwX3tCmtZXvDunPsGGJpp7Qww6KeexDk3CZgEEBkZqQvqc1HRsGCe7RFBjyZVeHrmZgZPS6Rnkyr88a4IypUo4nU8EclnOXmHkAJUz7BdDTicszjihabVSzPvyQ6M/FldFm49QufRCcz58muNvxAJMDkphA1AXTOrZWZhQD9gbu7EkvwWFhJEbOd6zH+yI9XLFmPk9C95ZFoiR06d8zqaiOST7F52+gGwBqhvZilmNtg5lwYMBxYDycAM59y2vIsq+aH+jTcwa1g7nu1+C6v2nqBLfALvrzvIZY2/EPF7Vpg+FoiMjHSJiYlexwgYB785y6hZm1m99xva1C7L3+5pTM3yxb2OJSLXyMySnHORV7udLvWULNUoV4z3HmnNS/c2Ytvh08SMTWBygoblifgrFYL8JDPj/pY1WBobTYebK/CXBcnc8+oqdh7VsDwRf6NCkGy5sVQ4k3/RgvH9m3Ho23P0GL+C0Ut3aVieiB9RIUi2mRl3NanCsrhoujeqzNhPd9Nj/Aq+OPit19FEJBeoEOSalS0exph+zZjyUCTf/5jGvRNX88L87Zy9kOZ1NBHJARWCXLc7GlRiSWwUD7SuwRsrvyJmzApW7znhdSwRuU4qBMmRG8JDeaF3I6YPaUOQwQNvrGPUR5s5dU7D8kQKGxWC5Io2tcux6KkoHouuzYzEQ3QZvZyl2zUsT6QwUSFIrgkPDeaZbrcw+4n2lCkWxqNvJzL8/Y2c+OG819FEJBtUCJLrGldLH5b3y871WLLtGJ3il/PxFykalidSwKkQJE+EBgfx5M/q8smIDtQqX5zYDzfx8NQNHP5Ow/JECioVguSpupVuYObQdvzxrgjW7jtJl9EJvLP2gIbliRRAKgTJc8FBxqD2tVgSG0XT6qX5/eyt9Ju0ln2pP3gdTUQyUCFIvqlethjvDG7Fy30as+PoabqNXcFry/eSdknjL0QKAhWC5Csz477I6iyLi+a2+hX428Id3P3qarYfPu11NJGAp0IQT1QsGc5rP2/Bqw8258ipc/R8ZSX/XLKT82mXvI4mErBUCOIZM+PORpVZGhtNzyZVGP/vPXQft5KkAxqWJ+IFFYJ4rkzxMOLvb8rUQS05d+ESfV5bzfPztnHmvIblieSnfCsEM6ttZm+a2cwMa0Fm9hczG29mA/MrixRMt9WvyOLYKAa0uYm3Vu2n65gEVuxO9TqWSMDIViGY2RQzO25mWzOtx5jZTjPbY2ajfuoxnHP7nHODMy33AqoCF4GUawku/qlEkRD+1KshMx5rS1hwEAPeXM9vZm7i1FkNyxPJa9l9hzAViMm4YGbBwASgGxAB9DezCDNrZGbzM31VzOJx6wNrnHNxwLDr2wXxR61qlWXByI4Mu60OH238mk6jl7No61GvY4n4tWwVgnMuATiZabkVsMf3nf8FYDrQyzm3xTnXI9PX8SweOgX4zxlEXV4i/0d4aDBPxzRgzhPtqVCiCEPfTeLx95I4/v2PXkcT8Us5OYdQFTiUYTvFt3ZFZlbOzF4DmpnZM77lWUBXMxsPJGRxvyFmlmhmiamp+jw5EDWsWoo5w9vz6671WZZ8nM7xCXyUpGF5IrktJAf3tSusZfl/qHPuG2BoprWzQObzCpnvNwmYBBAZGalXgAAVGhzEE7ffTNdbb+Tpjzbzy39tYs6mw/z17oZUK1PM63gifiEn7xBSgOoZtqsBh3MWR+Sn3VyxBP96rC3P97yVxP0n6To6gbfX7NewPJFckJNC2ADUNbNaZhYG9APm5k4skawFBRkD29VkSWwULWqW5Q9ztnHf62vYq2F5IjmS3ctOPwDWAPXNLMXMBjvn0oDhwGIgGZjhnNuWd1FF/q9qZYoxbVBL/tm3CbuP/0C3sSuY8NkeLmpYnsh1scJ0Yi4yMtIlJiZ6HUMKoOPf/8gf52xj4daj3FqlJC/d25iGVUt5HUukQDCzJOdc5NVup9EV4hcq3hDOxJ+3YOKDzTl2+jy9Jqzi5UU7+PGirmYWyS4VgviVbo0qsywuirubVeXVz/dy57gVJO7P/CM0InIlKgTxO6WLhfGPvk14++FWnL94mb6vr+G5uRqWJ3I1KgTxW1H1KrAkNoqBbWsybc1+uoxOYPku/XCjSFZUCOLXihcJ4bmetzJzaFvCQ4MYOGU9v5yxie/OXvA6mkiBo0KQgNDiprJ8MqIjw2+/mdlffk2n+AQWbjnidSyRAkWFIAEjPDSYX3Wtz9zh7bmxVBGGvbeRoe8kcfy0huWJgApBAtCtVUox+/H2PB3TgH/vPE6n+OXMSDykYXkS8FQIEpBCgoMYdlsdFo3sSIMbS/KbmZv5xZT1HDp51utoIp5RIUhAq12hBNOHtOHPvRuy8cC3dB2TwFurvuKShuVJAFIhSMALCjIGtLmJJXHRtKpVlufnbafva6vZfex7r6OJ5CsVgohP1dJFeeuhloy+vwn7Tpyh+7iVjP90t4blScBQIYhkYGbc3away+Ki6XxrJf65dBd3jV/JlpRTXkcTyXMqBJErKF+iCBMeaM7rA1pw8swFer+6ihcXJmtYnvg1FYLIT+h6640sjYumT/NqvL58H93GrmDdvm+8jiWSJ1QIIldRqmgoL/VpzLuDW5N2+TL3T1rLs7O38P2PF72OJpKrVAgi2dShbnkWPxXFw+1r8d66g3QdncBnO457HUsk1+RrIZhZbzObbGZzzKxLVmsiBVWxsBD+cFcEHw1rR/EiIQyauoHYD7/k5BkNy5PCL9uFYGZTzOy4mW3NtB5jZjvNbI+Zjfqpx3DOzXbOPQo8BNyf1ZpIQde8Rhnmj+jAiJ/VZd6mw3SOX878zYc1/kIKtWt5hzAViMm4YGbBwASgGxAB9DezCDNrZGbzM31VzHDXZ3334yprIgVWkZBg4jrXY96THahapijD3/+CIe8kcUzD8qSQynYhOOcSgMy/i7AVsMc5t885dwGYDvRyzm1xzvXI9HXc0r0ELHTObQS40ppIYXJL5ZLMGtaO397ZgIRdqXSKX8709Qf1bkEKnZyeQ6gKHMqwneJby8qTQCegj5kN/Ym1/8/MhphZopklpqbqt11JwRQSHMSQqDosfiqKiMolGTVrCw++sY4D35zxOppIttm1fBdjZjWB+c65hr7tvkBX59wjvu0BQCvn3JO5HxUiIyNdYmJiXjy0SK65fNnxwYaD/G3BDi5evsyvutRnUPtaBAeZ19EkQJlZknMu8mq3y+k7hBSgeobtasDhHD6mSKEWFGQ82PomlsRF0b5OeV74JJl7J65m51ENy5OCLaeFsAGoa2a1zCwM6AfMzXkskcKvcqmivDEwknH9m3Hw5Fl6jF/BmGW7uJCmYXlSMF3LZacfAGuA+maWYmaDnXNpwHBgMZAMzHDObcubqCKFj5nRs0kVlsZG0a1hZcYs281d41ey6dB3XkcT+R/XdA7BazqHIIXdsu3H+N3sLaR+f57BHWoR17k+RcOCvY4lfi6/ziGIyDXoFFGJpXHR3N+yBpNXfEXM2ATW7NWwPCkYVAgi+axkeCgv3tOI9x9tDUD/yWv57cdbOK1heeIxFYKIR9rVKc+ikVEMiarN9PUH6RKfwKfJx7yOJQFMhSDioaJhwfz2zluY9Xh7ShUNZfC0REZ88AXf/HDe62gSgFQIIgVA0+qlmfdkB2I71WPh1iN0Hp3AnC+/1vgLyVcqBJECIiwkiJGd6vLJiI7UKFuMkdO/5JFpiRw5dc7raBIgVAgiBUy9Sjfw0bB2PNv9FlbtPUHn+ATeW3eAy5f1bkHylgpBpAAKDjIe6VibJU9F07haKX738VYeeGMt+09oWJ7kHRWCSAFWo1wx3nukNS/d24hth0/TdUwCkxL2knZJ4y8k96kQRAo4M+P+ljVYFhdNVL0K/HXBDu6ZuJrkI6e9jiZ+RoUgUkhUKhnOpAEteOWBZnz97TnuGr+S+KW7OJ92yeto4idUCCKFiJnRo3EVlsVFc1eTKoz7NH1Y3hcHv/U6mvgBFYJIIVSmeBij72/KWw+15Psf07hn4mr+PH87Zy+keR1NCjEVgkghdnuDiiyJjeLB1jV4c+VXxIxZweo9J7yOJYWUCkGkkLshPJQXejdi+pA2BBk88MY6Rn20mVPnNCxPro0KQcRPtKldjkVPRfFYdG1mJB6ic/xylmw76nUsKURUCCJ+JDw0mGe63cLsJ9pTtngYQ95JYvj7GzmhYXmSDSoEET/UuFr6sLxfdq7Hkm3H6BS/nI+/SNGwPPlJ+VYIZtbbzCab2Rwz6+Jbq2Fmc81sipmNyq8sIoEgNDiIJ39Wl09GdKB2+eLEfriJh6du4OvvNCxPrixbheB7wT5uZlszrceY2U4z23O1F3Tn3Gzn3KPAQ8D9vuV6wCfOuYeBiGuPLyJXU7fSDfxraDv+0COCtftO0iV+Oe+s1bA8+V/ZfYcwFYjJuGBmwcAEoBvpL+b9zSzCzBqZ2fxMXxUz3PVZ3/0AvgD6mdm/gc9ysiMikrXgIOPhDrVYEhtFsxpl+P3srfSbtJZ9qT94HU0KkGwVgnMuATiZabkVsMc5t885dwGYDvRyzm1xzvXI9HXc0r0ELHTObfQ9xiDgj865O4DuubRPIpKF6mWL8c7gVrzcpzE7jp4mZuwKJn6uYXmSLifnEKoChzJsp/jWsvIk0AnoY2ZDfWuLgBFm9hqw/0p3MrMhZpZoZompqak5iCsikD7+4r7I6iyLi+b2+hV4adEOer+6iu2HNSwv0Fl2rzows5rAfOdcQ992X6Crc+4R3/YAoJVz7sm8iQqRkZEuMTExrx5eJCAt2HKEP8zZyndnLzI0ug7D77iZ8NBgr2NJLjKzJOdc5NVul5N3CClA9Qzb1YDDOXg8EfHAnY0qszQ2mp5Nq/DKZ3voPm4FSQcyf0IsgSAnhbABqGtmtcwsDOgHzM2dWCKSn8oUDyP+vqZMHdSSHy9eps9ra3hu7jbOnNewvECS3ctOPwDWAPXNLMXMBjvn0oDhwGIgGZjhnNuWd1FFJK/dVr8ii2OjGNDmJqau3k/XMQms2K1zd4Ei2+cQCgKdQxDJP+u/Osmojzaz78QZ+raoxrPdIyhVLNTrWHId8uMcgoj4sVa1yrJgZEcev60Os774mk6jl7Noq4bl+TMVgohkKTw0mN/ENGDOE+2pUKIIQ99N4vH3kjj+/Y9eR5M8oEIQkatqWLUUc4a359dd67Ms+Tid4xOYmaRhef5GhSAi2RIaHMQTt9/MghEdqVuxBL/61yYGvrWBlG/Peh1NcokKQUSuyc0VSzDjsbY83/NWEvefpMvoBKau+krD8vyACkFErllQkDGwXU2WxEYRWbMsz83bzn2vr2HPcQ3LK8xUCCJy3aqVKca0QS35R98m7D7+A3eOXcGEz/ZwUcPyCiUVgojkiJnRp0U1lsZF0SmiIn9fvJNer6xi69envI4m10iFICK5ouIN4bz6YAte+3lzUn84T68Jq3hp0Q5+vHjJ62iSTSoEEclVMQ0rsyw2mnuaVWXi53u5c+wKNuzXsLzCQIUgIrmuVLFQ/t63CW8/3IrzaZfp+9oa/jBnKz9oWF6BpkIQkTwTVa8CS2KjeKhdTd5Ze4CuoxNYvkvD8goqFYKI5KniRUJ4ruetzBzalvDQIAZOWU/cjC/57uwFr6NJJioEEckXLW4qyycjOjL89puZ++VhOsUvZ8GWIxp/UYCoEEQk34SHBvOrrvWZO7wDlUsV5fH3NjL03SSOn9awvIJAhSAi+S6iSkk+frwdo7o14POdqXSKX86MxEN6t+AxFYKIeCIkOIih0XVYOLIjDSqX5DczNzPgzfUcOqlheV5RIYiIp2pXKMH0R9vw594N+fLQd3QZncCUlV9xScPy8l2+FYKZ3WJmr5nZTDMblmG9uJklmVmP/MoiIgVLUJAxoM1NLImNonXtsvxp/nb6vLaa3ce+9zpaQMlWIZjZFDM7bmZbM63HmNlOM9tjZqN+6jGcc8nOuaHAfUDG3+35NDDjWoOLiP+pUroobz3UkjH3N2X/iTN0H7eS8Z/u1rC8fJLddwhTgZiMC2YWDEwAugERQH8zizCzRmY2P9NXRd99egIrgU99252A7cCxXNkbESn0zIzezaqyNC6arg1v5J9Ld3HX+JVsTvnO62h+z7J7Vt/MagLznXMNfdttgeecc119288AOOdezMZjfeKc625mfwGKk14o54C7nXNZfisQGRnpEhMTs5VXRPzD0u3HeHb2FlK/P8+jHWsT27ke4aHBXscqVMwsyTkXebXbheTgOaoChzJspwCtfyLQbcA9QBFgAYBz7ne+v3sIOHGlMjCzIcAQgBo1auQgrogURp0jKtGqVln+tjCZ1xP2sXjbUf52b2Pa1C7ndTS/k5OTynaFtSzfbjjnPnfOjXDOPeacm5Dp76Y65+Zncb9JzrlI51xkhQoVchBXRAqrUkVDefGexrz/SGsuO+g3aS2/+3gL3/940etofiUnhZACVM+wXQ04nLM4IiJZa3dzeRY91ZHBHWrx/vqDdBmdwL936BRkbslJIWwA6ppZLTMLA/oBc3MnlojIlRULC+H3PSKYNawdN4SH8PDURJ6a/gUnz2hYXk5l97LTD4A1QH0zSzGzwc65NGA4sBhIBmY457blXVQRkf9qVqMM857swMif1WX+5iN0jl/OvE2HNf4iB7J9lVFBoKuMRORKdhw9zdMzN7Mp5RSdbqnEX+5uSKWS4V7HKjCye5WRRleISKHX4MaSfDSsHb+78xZW7E4fljd9/UG9W7hGKgQR8QshwUE8GlWbxU9FcWuVkoyatYUH31jHgW/OeB2t0FAhiIhfqVm+OO8/0oa/3t2ILSmn6DomgTdW7NOwvGxQIYiI3wkKMh5oXYMlcVG0r1OeFz5J5p6Jq9l5VMPyfooKQUT8VuVSRXljYCRj+zXl0Mmz9Bi/gjHLdnEhTcPyrkSFICJ+zczo1bQqS2OjuLNRZcYs281d41ey6ZCG5WWmQhCRgFCuRBHG9mvGmwMjOXXuIne/uoq/fLKdcxcueR2twFAhiEhA+dktlVgSF0W/VjWYvOIrYsYmsGbvN17HKhBUCCIScEqGh/LXuxvxwaNtAOg/eS3PzNrC6QAflqdCEJGA1bZOORaNjGJIVG0+3HCQLvEJfJocuMPyVAgiEtCKhgXz2ztv4ePH21O6WCiDpyUy4oMv+OaH815Hy3cqBBERoEn10swd3oHYTvVYuPUInUcnMOfLrwNq/IUKQUTEJywkiJGd6vLJiI7UKFuMkdO/5JFpiRw5dc7raPlChSAikkm9Sjfw0bB2PNv9FlbtPUHn+ATeW3eAy34+/kKFICJyBcFBxiMda7PkqWgaVyvF7z7eygNvrGX/Cf8dlqdCEBH5CTXKFeO9R1rzt3sase3wabqOSWBSwl7SLvnf+AsVgojIVZgZ/VrVYFlcNFH1KvDXBTu4Z+Jqko+c9jparsrXQjCz4maWZGY9MmxPM7PJZvZgfmYREblWlUqGM2lAC8b3b8bX357jrvEriV+6i/Np/jH+Iru/U3mKmR03s62Z1mPMbKeZ7TGzUdl4qKeBGRm27wFmOuceBXpmO7WIiEfMjLuaVGFZXDQ9m1Rh3Ke76TFuJRsPfuutoNKcAAAHfUlEQVR1tBzL7juEqUBMxgUzCwYmAN2ACKC/mUWYWSMzm5/pq6KZdQK2Axl/DLAacMj3Z/+oWBEJCGWKhxF/f1PeGtSSM+fTuHfiav40bztnL6R5He26hWTnRs65BDOrmWm5FbDHObcPwMymA72ccy8CPTI/hpndDhQnvTzOmdkCIIX0UvgSnc8QkULo9voVWRwbxcuLdjJl1VcsTT7Ki3c3pkPd8l5Hu2Y5eRGuyn+/u4f0F/eqWd3YOfc759xTwPvAZOfcZWAWcK+ZTQTmXel+ZjbEzBLNLDE1NTUHcUVE8sYN4aH8uXdDPhzShpCgIH7+5jqenrmZU+cK17C8bL1DyIJdYe2qP7XhnJua4c9ngEFXuf0kYBJAZGSkf/9UiIgUaq1rl2PhyI6MXraLyQn7+GzncV7o3ZAut97odbRsyck7hBSgeobtasDhnMURESncwkODeabbLcx+oj1li4cx5J0knnh/I6nfF/xheTkphA1AXTOrZWZhQD9gbu7EEhEp3BpXK828Jzvwqy71WLrtGJ1HL2fWxpQCPSwvu5edfgCsAeqbWYqZDXbOpQHDgcVAMjDDObct76KKiBQuocFBDL+jLgtGdqB2+eLEzdjEoKkb+Pq7gjkszwpyW2UWGRnpEhMTvY4hInLNLl12vL1mPy8v2kmQwahuDXiw9U0EBV3pdGzuMrMk51zk1W6nSz1FRPJBcJAxqH0tlsRG0fymMvx+zjb6TVrLvtQfvI72/6kQRETyUfWyxXj74Vb8vU9jdhw9TczYFbz6+Z4CMSxPhSAiks/MjL6R1VkWF83t9Svw8qKd9H51FdsOn/I0lwpBRMQjFUuG8/qASCY+2Jyjp87T85VV/H3xDn686M0kHxWCiIjHujWqzLK4KHo3rcqEz/bSfdwKEvefzPccKgQRkQKgdLEw/nlfE6Y93IofL16m7+treG7uNs6cz79heSoEEZECJLpeBRbHRvGLNjcxbc1+uoxOIGFX/sxxUyGIiBQwJYqE8Hyvhsx4rC1FQoP4xZT1PDt7S54/b06G24mISB5qWbMsC0Z0ZNynuylbPCzPn0+FICJSgIWHBvObmAb58lz6yEhERAAVgoiI+KgQREQEUCGIiIiPCkFERAAVgoiI+KgQREQEUCGIiIhPofoVmmaWChzIwUOUB07kUpzCIhD3GQJzvwNxnyEw9/ta9/km51yFq92oUBVCTplZYnZ+r6g/CcR9hsDc70DcZwjM/c6rfdZHRiIiAqgQRETEJ9AKYZLXATwQiPsMgbnfgbjPEJj7nSf7HFDnEEREJGuB9g5BRESyEBCFYGYxZrbTzPaY2Siv8+QVM6tuZp+ZWbKZbTOzkb71sma21Mx2+/5bxuusuc3Mgs3sCzOb79uuZWbrfPv8oZnl/W8XyUdmVtrMZprZDt/xbhsgxznW9297q5l9YGbh/niszWyKmR03s60Z1q54fC3dON/r22Yza369z+v3hWBmwcAEoBsQAfQ3swhvU+WZNOCXzrlbgDbAE759HQV86pyrC3zq2/Y3I4HkDNsvAaN9+/wtMNiTVHlnLLDIOdcAaEL6vvv1cTazqsAIINI51xAIBvrhn8d6KhCTaS2r49sNqOv7GgJMvN4n9ftCAFoBe5xz+5xzF4DpQC+PM+UJ59wR59xG35+/J/1Foirp+zvNd7NpQG9vEuYNM6sGdAfe8G0bcAcw03cTv9pnMysJRAFvAjjnLjjnvsPPj7NPCFDUzEKAYsAR/PBYO+cSgJOZlrM6vr2At126tUBpM6t8Pc8bCIVQFTiUYTvFt+bXzKwm0AxYB1Ryzh2B9NIAKnqXLE+MAX4DXPZtlwO+c86l+bb97ZjXBlKBt3wfk71hZsXx8+PsnPsa+AdwkPQiOAUk4d/HOqOsjm+uvcYFQiHYFdb8+tIqMysBfAQ85Zw77XWevGRmPYDjzrmkjMtXuKk/HfMQoDkw0TnXDDiDn308dCW+z8x7AbWAKkBx0j8uycyfjnV25Nq/90AohBSgeobtasBhj7LkOTMLJb0M3nPOzfItH/vPW0jff497lS8PtAd6mtl+0j8OvIP0dwylfR8rgP8d8xQgxTm3zrc9k/SC8OfjDNAJ+Mo5l+qcuwjMAtrh38c6o6yOb669xgVCIWwA6vquRAgj/STUXI8z5QnfZ+dvAsnOufgMfzUXGOj780BgTn5nyyvOuWecc9WcczVJP7b/ds49CHwG9PHdzN/2+ShwyMzq+5Z+BmzHj4+zz0GgjZkV8/1b/89+++2xziSr4zsX+IXvaqM2wKn/fLR0rQLiB9PM7E7Sv2sMBqY45/7icaQ8YWYdgBXAFv77efpvST+PMAOoQfr/VH2dc5lPWBV6ZnYb8CvnXA8zq036O4aywBfAz51z573Ml5vMrCnpJ9HDgH3AINK/wfPr42xmzwP3k35F3RfAI6R/Xu5Xx9rMPgBuI32q6THgj8BsrnB8feX4CulXJZ0FBjnnEq/reQOhEERE5OoC4SMjERHJBhWCiIgAKgQREfFRIYiICKBCEBERHxWCiIgAKgQREfFRIYiICAD/Dy+bKdvLO3rJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sigma = np.exp(-np.arange(n))\n", "pt.semilogy(sigma)\n", "\n", "A = (U0 * sigma).dot(VT0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Find the approximate range" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's fix parameters first. What accuracy should we obtain for the values of $k$ and $p$ below? (where `p` represents the 'extra' dimensions)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "k = 10\n", "p = 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Draw a random Gaussian matrix `Omega` and obtain orthogonal columns in a matrix `Q` spanning the range:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Omega = np.random.randn(n, k+p)\n", "\n", "Y = A @ Omega\n", "\n", "Q, _ = la.qr(Y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an alternative to the above, use a few iterations of the power method on `Omega`:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Omega = np.random.randn(n, k+p)\n", "\n", "Y = A @ Omega\n", "\n", "Y = A @ A.T @ Y\n", "\n", "Q, _ = la.qr(Y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Observations about associativity?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Convert to factorization form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reconstruct $C$ in the factorization-form LRA $A\\approx BC$:\n", "\n", "(Recall $A\\approx QQ^TA$, $B=Q$)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(15, 100)\n" ] } ], "source": [ "C = Q.T @ A\n", "\n", "print(C.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sanity-check that form:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3.7208685050055835e-06" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Aapprox1 = Q.dot(C)\n", "\n", "la.norm(A - Aapprox1, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute the approximate SVD" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the SVD of $C=U_C \\Sigma_C V_C^T$:\n", "\n", "(Make sure to pass `full_matrices=False` to the SVD.)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "UC, sigmaC, VTC = la.svd(C, full_matrices=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reconstruct the SVD of $A$: $A\\approx QU_C \\Sigma_C V_C^T$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "UAapprox = Q @ UC\n", "\n", "sigmaAapprox = sigmaC\n", "\n", "VTAapprox = VTC" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the 2-norm of the reconstructed $A$ with the original $A$:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3.7208685050056457e-06" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Aapprox = (UAapprox*sigmaAapprox) @ VTAapprox\n", "\n", "la.norm(A - Aapprox, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimate the error\n", "\n", "Compute an a-posteriori estimate of approximation error in the spectral norm:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.4702805615731663e-07" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "omega = np.random.randn(n)\n", "\n", "Aomega = A @ omega\n", "\n", "err = Aomega - Q @ Q.T @ Aomega\n", "\n", "la.norm(err, 2) / la.norm(omega, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Is this the right direction for the error estimator?\n", "* Is the estimator supposed to be conservative? Can it be?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }