{
  "metadata": {
    "name": "",
    "signature": "sha256:aac1bb3788907f361665f9126a7c1933c6af9a1e35dceededf8617edb19beb19"
  },
  "nbformat": 3,
  "nbformat_minor": 0,
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "# Randomized SVD"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "import numpy as np\n",
            "import numpy.linalg as la\n",
            "import matplotlib.pyplot as pt"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 2
        },
        {
          "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",
          "collapsed": false,
          "input": [
            "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)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "output_type": "stream",
              "stream": "stdout",
              "text": [
                "2.67019320578e-13\n"
              ]
            },
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 3,
              "text": [
                "[<matplotlib.lines.Line2D at 0x7fd35095af98>]"
              ]
            },
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGu5JREFUeJzt3Xu81VP+x/HXp3tSU9QQZZIuP6UQiZhxptvJJdU0iVya\nCo0ImQn9hMK4DRpEv9+PNLnkmlK5dNOZMpTRXSIxIXIqkkshtX5/rN10pE7n7LP3Xnt/v+/n47Ef\nj305e+/P+T7q0+qz1vosc84hIiLRUC50ACIikjpK6iIiEaKkLiISIUrqIiIRoqQuIhIhSuoiIhFS\nbFI3s/pmNtvMlpvZW2Z2WeL5/cxshpmtNLPpZlYzM+GKiEhxrLh16mZ2IHCgc26xme0LLAC6AX2B\nDc65O8zsaqCWc+6ajEQsIiJ7VOxI3Tn3mXNuceL+N8AK4GDgDGBc4sfG4RO9iIgEVuKaupk1AI4G\n5gMHOOcKEy8VAgekPDIRESm1EiX1ROllAnC5c+7roq85X79RrwERkSxQYW8/YGYV8Qn9UefcpMTT\nhWZ2oHPuMzOrC6zbzfuU6EVEkuCcs2Tfu7fVLwaMAd52zv2tyEuTgT6J+32ASbu+NxGYbs5xww03\nBI8hW266FroWuhbF38pqbyP1E4FzgaVmtijx3FDgNuBpM+sPrAbOLHMkIiJSZsUmdefcq+x5NN8h\n9eGIiEhZaEdpBuTl5YUOIWvoWuyka7GTrkXqFLv5qEwfbObS9dkiIlFlZrh0TZSKiEhuUVIXEYkQ\nJXURkQhRUhcRiZC0JvWXXoJt29L5DSIiUlRak/oNN8Ahh8DQobBxYzq/SUREIM1J/Y03YPp0WLMG\nunSBLVvS+W0iIpKRderbt0Pv3r4U89RTUE6VfBGR3cqJderlysHf/w6FhTBkSCa+UUQknjI2Zq5S\nBSZNghdfhAceyNS3iojES8bbBLz/Phx/PLz6KjRtmpavFhHJWTlRfinqsMPg+uuhXz8tdxQRSbUg\nU5aXXALly8N994X4dhGR6ArWpfG99+CEE2DePGjUKC0hiIjknJwrv+zQuDFcey307QtffRUqChGR\naAm6Yvyyy6BZMzjiCN9SQEREyiYrDsmYNQsuvBBOOgnuvRdq1kxLSCIiWS9nyy9FtW8Py5ZB5cqQ\nn69yjIhIsrJipL6DczBwICxf7ssx1aqlJTQRkaxV1pF6ViV18H1i+vWDTz6BKVP8TlQRkbiIXFIH\nvynpnHN8gn/qKbCkfz0RkdwSiZr6rsqXh3HjYMUKeOyx0NGIiOSOrByp77BokZ84XbAA6tdPUWAi\nIlkskiP1HY4+2q9l79/fT6KKiEjxsjqpA1xzDWzaBKNHh45ERCT7ZXX5ZYd334UTT/SblI48MiUf\nKSKSlSJdftmhaVMYNQrOOMOfniQiIruXE0kd4Kyz4A9/gO7d4bvvQkcjIpKdcqL8ssP27dCrF1St\n6pc8av26iERNLMovO5Qr55P522/7CVStiBER+amcSuoA++wDL7/sJ00HDvSjdxER8XIuqQPUrg2v\nvOJ3nJ53HmzdGjoiEZHskJNJHaBGDd/JcdMmvyrmiy9CRyQiEl7OJnXwE6YTJ0Lz5tCqlT/vVEQk\nznJq9UtxJk/2pydddRUMHuwnVUVEck0kW+8ma/VqOPtsX5oZOxYOOiijXy8iUmaxWtK4Nw0awNy5\n0Latbwb27LOhIxIRyaxIjdSLmj8fzj0Xzj8frrsuWBgiIqWi8ksxCguhdWsYORJ69AgaiohIiSip\n78WCBdC5s9+s1LJl6GhERIqnmvpeHHMM3HsvdO0KGzaEjkZEJL0in9TBr4jp1csn9k2bQkcjIpI+\nsUjqALfc4g/Y+O1vYd260NGIiKRHbJJ6uXJw//3QpQucdJJf0y4iEjWxSerg+6+PGAGDBsGvfw1P\nP60ujyISLXtN6mb2sJkVmtmyIs8NN7M1ZrYoceuc3jBTa9Ag+Pvf4c47/SalyZPVm11EomGvSxrN\n7NfAN8AjzrkWieduAL52zt1dzPuyYkljcZyDKVP85qTKleHGGyE/XycqiUg4aV/S6JybC2zc3Xcn\n+6XZwsy37V20CIYMgSuv9PX2uXNDRyYikpyy1NQHmdkSMxtjZjVTFlEA5cpBz56wbBlcconfffra\na6GjEhEpvQpJvm80cGPi/k3AXUD/XX9o+PDh/7mfl5dHXl5ekl+XGeXLQ+/eUKuWT/Lz50O9eqGj\nEpEoKygooKCgIGWfV6I2AWbWAJiyo6ZektdyoaZenNtv910e58zxh3GIiGRCkDYBZla3yMPuwLI9\n/WyuuuoqaNTIH7yRw/82iUjMlGT1yxPAyUBtoBC4AcgDjgIc8G9ggHOucJf35fRIHWDzZujYEapU\ngf/5H2jcOHREIhJ16tKYZj/+6BuC3XILXHGFH8FXqhQ6KhGJKnVpTLMKFfxSxwUL/MRpkybwwAOw\nZUvoyEREfk4j9VJ6/XW49Vb417/88sf8fL8rtUKy64hERIpQ+SWQpUvhwQehoAA+/thvWrr/fvjV\nr0JHJiK5TEk9C2zY4PvIzJ/vT1gqp6KWiCRJNfUsULs2/OUv8P33vt4uIhKKRuoptHIltG0L8+b5\nNe4iIqWlkXoWadIEhg2DP/wBtm0LHY2IxJGSeopddpnvIXPxxbBmTehoRCRulNRTrFw5f6JS1arQ\nsqVvELZwYeioRCQulNTT4IAD4J574IMPoFUrOPVUGD06dFQiEgeaKM2A99/3ib1rV7jtNi15FJE9\n0zr1HPH559CtGxx0EIwZA/vuGzoiEclGWv2SI/bfH2bMgGrV/CqZ0aNh69bQUYlI1CipZ1CVKvDw\nw/6w64kToVkzmDw5dFQiEiUqvwQ0cyYMGADt2sHIkSrJiIjKLzmtQwdYvNj3bG/Vynd+FBEpCyX1\nwKpXh7Fjfe+Y007za9xFRJKlLuBZomdPaNrUL31cv973ahcRKS0l9SzSsiXMnQudOsG6dTB8OFjS\nlTURiSOVX7LMoYfCq6/6VTG33RY6GhHJNVr9kqU++QRat4Zx46Bjx9DRiEimaPVLRB18MIwfD+ed\nBx9+GDoaEckVSupZLC8PhgyBHj3gu+9CRyMiuUBJPctdeSU0bAinnw6vvAKqaIlIcVRTzwHffefX\nst9/vz9RaeBAf7pS9eqhIxORVFOXxhhxzi95vO8+mD0bLrgABg3y9XcRiQZNlMaIGfzmN/DMM76l\nwJYt0KIFjBgB27eHjk5EsoFG6jlu7Vro1Qtq1IBHH4VatUJHJCJloZF6zNWtC7NmQaNGcOyxsHRp\n6IhEJCSN1CNk/HgYPNj3az/uuNDRiEgyNFEqP/HCC9Cvn28z0KZN6GhEpLRUfpGfOO00v/yxSxeY\nNy90NCKSaUrqEXTqqfDIIz6xDxvmW/mKSDwoqUdU584wfz58/rnv03755UruInGgpB5hDRvC6NGw\nfDn88IM/C/Xzz0NHJSLppInSmHAOrr7a70SdNcuvaxeR7KPVL1Jizvlj8pYtg5dfhmrVQkckIrvS\n6hcpMTMYNcqXZX77W5g2TV0fRaJGI/UY2rbNb1S6/XaoVAmuugq6doWqVUNHJiIqv0jStm+HF1+E\nu++GN9/0h3Kcfjr07KkeMiKhKKlLSnzxhS/HTJwIr7/uz0Zt1y50VCLxo6QuKTd9OvTtC717w803\nQ+XKoSMSiQ9NlErKdeoES5bAqlVw1FHw4IOweXPoqESkJDRSlz1yzq9pv/deX5Lp18+3HmjdWqN3\nkXTRSF3Sxgw6dPAdH+fN84+vuAJq1/b19mnTQkcoIrvSSF1KbdMmmDEDLr4Ypk5Vi1+RVNJEqQQz\ndSpceCG8+iocdljoaESiIe3lFzN72MwKzWxZkef2M7MZZrbSzKabWc1kA5DcdfrpcN11vtWvGoWJ\nZIeS1NTHAp13ee4aYIZzrgkwK/FYYmjgQL8b9ZRT4MMPQ0cjIntN6s65ucDGXZ4+AxiXuD8O6Jbi\nuCSH3HYb9OjhV8U8+WToaETirUQ1dTNrAExxzrVIPN7onKuVuG/AFzseF3mPauoxs2CB37DUpo3v\nK1O3buiIRHJPWWvqFcoagHPOmdlus/fw4cP/cz8vL4+8vLyyfp1ksWOOgYUL/RF6zZr5HjJ//jM0\naRI6MpHsVVBQQEFBQco+L9mR+jtAnnPuMzOrC8x2zv3XLu/RSD3G1q/3bX4feMAn95EjtWFJpCRC\nbT6aDPRJ3O8DTEo2AImmOnVgxAh4/31Yu9ZvVlq7NnRUItFXkiWNTwCvAU3N7GMz6wvcBnQ0s5VA\nu8RjkZ+pUQMmTPAHYbdu7dsNiEj6aPORZMzkyX6zUvPmcNFF0L27SjIiu9KOUskp338Pzz8P//d/\nsHQpPPEEtG8fOiqR7KGkLjmroADOPBNeeMGXZkREXRolh+XlwUMP+Xa+77wTOhqRaCjzOnWRsjjj\nDPjyS8jPh0mT4PDDoUqV0FGJ5C4ldQnu/PPh6699q4FPPoGaNaFFCxgzBn71q9DRieQW1dQlq2zf\nDoWFMH483HUXPPccHH986KhEMkcTpRJZU6f6A7BHjYJevUJHI5IZSuoSaUuW+Lp7w4Z+XXv37lC/\nfuioRNJHq18k0o480q+MGTwYFi2Co4/2h3J88knoyESyk5K6ZL2qVf1ofexY3z/m+OOhVSt45pnQ\nkYlkH5VfJCfNnw/nnecT/D33QK1ae3+PSC5Q+UViqU0bX46pUQOOOMKvcRcRjdQlAubMgQsu8PX2\nv/1NJy5JbtNIXWLvN7/xq2QaNvSj9uuug6++Ch2VSBhK6hIJVavCrbf6kszHH/sj9EaP9puZROJE\n5ReJpKVL4Y9/hEqVfNOwRo1CRyRSMiq/iOxGy5Ywdy507epXyIwcCVu3ho5KJP2U1CWyypf3m5Ze\nf933bG/Rwq+S0X8gJcpUfpFYcA6mTYOrroJf/MKvlmnb1pdlLOn/6Iqknnq/iJTCtm3+CL2pU+G1\n12DLFr9b9c47tYFJsoNq6iKlUL48nHsuPPkkfPSRXy1TrZqvwc+cGTo6kbLTSF0EmD4d+vWDbt3g\n4ouhWTOVZSQMjdRFUqBTJ78MskIF3wWyYUO49FJYuTJ0ZCKlo6QukrDffr7NwOrVMGUK1KkDJ53k\n74vkCpVfRIoxbx707AkXXgjDhkE5DYMkzbT6RSTN1q6F3//ej+QfftiP4EXSRTV1kTSrWxdmz/aT\np0ce6TcyiWQrjdRFSmHOHOjTBzp08BuYWrWCihVDRyVRopG6SAbtaPNbq5avs++/P+Tnw1//CsuW\nqQWBhKeRukgZbNjgG4fNnAkvvQTffw9nngk33+w3NYmUliZKRbKEc7BqFdx0EyxYAE8/Dc2bh45K\nco3KLyJZwgwaN4ZHHoEhQyAvD8aMCR2VxI1G6iJp8vbb0L07DBgAV14ZOhrJFWUdqVdIZTAislOz\nZjBjBpx4ItSr52vtIummpC6SRocc4tv8duwIBx7oV8+IpJNq6iJpduSRMH68bzewZEnoaCTqlNRF\nMqBDBxg1Ctq3h9tvhx9/DB2RRJUmSkUyaPVq6N8fvv3W95Fp1ix0RJJttKRRJIc0aOAnT/v08fX1\nDh3g0Ud9khdJBSV1kQwrV86frrRmjV/u+NRTUL++1rRLaqj8IpIFVqzw7X3btPG19332CR2RhKLy\ni0gEHH44vPEG/PADnHCC7yWzbp0ahEnpaaQukkWcg4cegnHj/I5UM2jb1j93wAGho5NMUEMvkYhy\nzo/WR43y69xfegmaNAkdlaSbkrpIDIwZA9deC88950fuEl2qqYvEQP/+MHYsdO0Kd9zh+7aL7I6S\nukiOOOUUeO01+Oc//aaliRM1kSo/V6byi5mtBr4CtgFbnXPHFXlN5ReRNJk507fzPeQQeOwxqFkz\ndESSKqHLLw7Ic84dXTShi0h6dejgT1c67DBo3RqWLw8dkWSLVJRfkv4XRUSSV7Ei3HMPDBvmT1ka\nPx62bAkdlYRW1vLLB8AmfPnlf51zDxZ5TeUXkQx580249FJYtsyfi3rCCX5ytWXL0JFJaYU++ehE\n59xaM6sDzDCzd5xzc3e8OHz48P/8YF5eHnl5eWX8OhHZnWOPhXnzYPNmX5b5xz8gPx9OPhlGjICm\nTUNHKHtSUFBAQUFByj4vZevUzewG4Bvn3F2JxxqpiwT0zTdw331w991w0kn+9KV27XyCNxVNs1aw\niVIz28fMqifuVwM6AcuS/TwRSa1994WhQ+G996BHD1+iyc+Hhg1h9uzQ0Um6JD1SN7NDgYmJhxWA\nx51ztxZ5XSN1kSzjHEybBn37+va/114L5cuHjkqKUpsAESm1Tz+Fc87xCb1PH6ha1bf7bdbMH+Qh\n4Sipi0hStm2DkSNh8WK/FPLbb/0k6333wVlnhY4uvpTURSRlFi+G3/0OunXzB2RXrBg6ovgJvaNU\nRCLkqKP8hOqKFX7X6qefho5ISktJXUR+Yr/9YOpUaN8ejjnGT6xK7lD5RUT26B//gHPPhd694cYb\noXLl0BFFn8ovIpI2J58MCxfCypVw4IF+xcyECX5SVbKTRuoiUiJr18Lzz/uk/u9/Q0EB1KsXOqro\n0eoXEcm4v/7VH4ZdUAB164aOJlpCN/QSkRgaMsQfqde+vU/sv/xl6IhkByV1EUnKsGHwww++l/vA\ngb7+3rw5lNNMXVAqv4hI0pyDZ5/1yx7nzIHPP4cBA+D666FKldDR5SbV1EUka3z0EfzpT/6wjjFj\n4MQTQ0eUe5TURSTrTJgAgwb5dgMjRkCdOqEjyh1apy4iWadHD3jrLahQAQ4/HP7yF38qk6SfRuoi\nklarVvm+7QUFcMghUKmS35nauzdccEHo6LKPyi8ikhM++MBPpH7/PXz1FVx+OfTs6UfxOl5vJyV1\nEclJ69dDly7QuLGfVK1UKXRE2UE1dRHJSXXqwCuv+FH7CSfAHXfAkiV+maQkTyN1EQlq2zZ44QW/\n1n36dJ/k69WDGjX87ayz4OyzQ0eZOSq/iEikrFkDhYU+ua9fD3/+M1x5JVxxRejIMkNJXUQi7cMP\noWNHP2IfMSL6k6pK6iISeevWQX4+tG4Nd97pyzJRpYlSEYm8X/7Sr3P/4Qdo2hQefNDX4uXnNFIX\nkZzy5psweDBs2uRPYsrLg1atoGLF0JGlhsovIhI7zsHLL/tbQYE/iemcc/yyyOrVQ0dXNkrqIhJ7\nGzbA0KEwaxaMHet7u+cqJXURkYQXXoCLLoJeveDmm2GffUJHVHqaKBURSTjtNFi61K9zb9ECZs4M\nHVHmaaQuIpH04otw8cXQrp1vGnbQQaEjKhmN1EVEduPUU31P9/33hyOOgD/+0XeKjDoldRGJrOrV\n/Wald9+F2rXhuON8y4Eor3FX+UVEYmPjRn8qU61a8NhjULVq6Ih+TuUXEZESqlULXnrJ927v1Am+\n+CJ0RKmnpC4isVK5Mjz+uC/FtG0LCxaEjii1VH4Rkdh6/HFfY7/oIrjuuuw4fUmbj0REymDtWhgw\nAFav9q0GGjeGRo1847DKlTMfj5K6iEgZOQfPPw+vvgrvvedvX34J118P/ftntlmYkrqISBosWABX\nXw0ffww33QTdumWmPKOkLiKSRjNm+KT+1lvQpQv8/vdwyilQoUJ6vk9JXUQkAz79FCZMgPHj4euv\n4d57fQuCVFNSFxHJIOdg4kS/aqZ1a7jlFj+5mipK6iIiAWzZ4g/luP9+34LglFP8AdkNGvjj92rV\nSu6QbCV1EZGAtm+HhQt9V8jZs32ZZt062LzZr3+/667STbAqqYuIZKGNG6FvX5/gn3kGDj64ZO9T\n7xcRkSxUqxY895w/uKN1a5gzJzPfq5G6iEiaTZsG558PI0b4vu7FCTZSN7POZvaOmb1nZlcn+zki\nIlGXnw///KdfBjlwIGzdmr7vSiqpm1l5YBTQGWgGnG1mh6cysCgpKCgIHULW0LXYSddipzhci0aN\nYN48v0O1Q4f0ncKU7Ej9OGCVc261c24r8CTQNXVhRUsc/sCWlK7FTroWO8XlWtSoAZMm+aP2jjvO\n95bZvDm135FsUj8Y+LjI4zWJ50REpBjly/ueMosX+8ZhzZrB3Lmp+/xkuxdoBlREpAzq1YMnnoCC\nAjjggNR9blKrX8zseGC4c65z4vFQYLtz7vYiP6PELyKShIxvPjKzCsC7QHvgU+AN4Gzn3IpkAxER\nkbJLqvzinPvRzC4FpgHlgTFK6CIi4aVt85GIiGReWtoExHljkpnVN7PZZrbczN4ys8sSz+9nZjPM\nbKWZTTezmqFjzRQzK29mi8xsSuJxLK+FmdU0s2fNbIWZvW1mbWJ8LYYm/o4sM7PxZlY5LtfCzB42\ns0IzW1bkuT3+7olr9V4ip3ba2+enPKlrYxJbgcHOuebA8cAlid//GmCGc64JMCvxOC4uB95m56qp\nuF6Le4AXnXOHAy2Bd4jhtTCzBsCFQCvnXAt8Cfcs4nMtxuLzY1G7/d3NrBnQC59LOwMPmFmxeTsd\nI/VYb0xyzn3mnFucuP8NsAK/hv8MYFzix8YB3cJEmFlmVg84FXgI2DGjH7trYWa/AH7tnHsY/LyU\nc24TMbwWwFf4wc8+iUUX++AXXMTiWjjn5gIbd3l6T797V+AJ59xW59xqYBU+x+5ROpK6NiYlJEYk\nRwPzgQOcc4WJlwqBFK5MzWojgSHA9iLPxfFaHAqsN7OxZrbQzB40s2rE8Fo4574A7gI+wifzL51z\nM4jhtShiT7/7QfgcusNe82k6krpmXgEz2xeYAFzunPu66GuJ9pWRv05mdjqwzjm3iJ2j9J+Iy7XA\nrzRrBTzgnGsFfMsu5YW4XAszOwy4AmiAT1r7mtm5RX8mLtdid0rwuxd7XdKR1D8B6hd5XJ+f/ksT\neWZWEZ/QH3XOTUo8XWhmByZerwusCxVfBrUFzjCzfwNPAO3M7FHieS3WAGucc/9KPH4Wn+Q/i+G1\nOBZ4zTn3uXPuR+A54ATieS122NPfiV3zab3Ec3uUjqT+JtDYzBqYWSV8kX9yGr4nK5mZAWOAt51z\nfyvy0mSgT+J+H2DSru+NGufcfzvn6jvnDsVPhL3inDuPeF6Lz4CPzaxJ4qkOwHJgCjG7FvgJ4uPN\nrGri70sH/ER6HK/FDnv6OzEZOMvMKpnZoUBj/GbPPXPOpfwGnILfcboKGJqO78jWG3ASvn68GFiU\nuHUG9gNmAiuB6UDN0LFm+LqcDExO3I/ltQCOBP4FLMGPTn8R42txFf4ftWX4icGKcbkW+P+1fgr8\ngJ9/7Fvc7w78dyKXvgPk7+3ztflIRCRCdEapiEiEKKmLiESIkrqISIQoqYuIRIiSuohIhCipi4hE\niJK6iEiEKKmLiETI/wNH69Pu2CCf4AAAAABJRU5ErkJggg==\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd352980358>"
              ]
            }
          ],
          "prompt_number": 3
        },
        {
          "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",
          "collapsed": false,
          "input": [
            "sigma = np.exp(-np.arange(n))\n",
            "pt.semilogy(sigma)\n",
            "\n",
            "A = (U0 * sigma).dot(VT0)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "display_data",
              "png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD7CAYAAACIYvgKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+81/P9//Hbo4MwP6KRlcM7KUQ/UB016Y3MySSEZPWt\nqBnLj/34qtmGbebnLKY1oZL2UbHYataSH+/EGSqVM5XFvKlsMeMzzNavx+eP5xtvRyfvzvvH6/3j\nfr1culy8Xr1fr9fjvHTO8zxez+fr8TB3R0REKluzqAMQEZHoaTAQERENBiIiosFARETQYCAiImgw\nEBERYIeoA9gaM9N6VxGRJnB3a8pxRZsZzJ3r7L+/c/HFznvvOe6V+efqq6+OPIZi+aN7oXuhe7Ht\nP9ko2sGgthbq6+GDD6BLF3jyyagjEhEpX0U7GAC0aAH33AO33gqDB8Pll8O//x11VCIi5aeoB4OP\n9O8PL7wAb70FXbtCXV3UERVOPB6POoSioXvxCd2LT+he5IZl+5wpH8zMG4vrwQfhm9+EIUPgJz+B\nnXcucHAiIkXKzPBym0BuzJlnhiwhmYQjj4Tnnos6IhGR0ldymUG6mTPh0kth5Ei46ipo3rwAwYmI\nFKmKygzSDRoEy5fDiy9Ct27w/PNRRyQiUppKejAA2G8/eOghGDMG+vWDa66BDRuijkpEpLQUfDAw\ns0PN7Fdmdr+ZXZCbc4YJ5aVLYdEiqKkJ8woiIpKZyOYMzKwZMMPdz9nK32U0Z7A17uHdhCuuCO8l\njBkDOxRl0Q0RkdyKZM7AzCab2Xozq2+wv9bMVpnZajMb08ix/YGHgRlNvX7jccGIEWH+4MknoWfP\nMKcgIiKNy+Yx0RSgNn2HmVUB41P7OwKDzewwMxtqZuPMrDWAu89x937AsCyuv03V1fDHP8KoURCP\nw403wubN+bqaiEhpy+oxkZnFgDnu3im13RO42t1rU9tjAdz9hrRj+gBnAjsDK9391q2ct8mPibYm\nmYQLLgilLKZMgUMPzdmpRUSKRjaPiXL9NL0NsCZtey1Qk/4Bd18ALPi8E8XjcWKxGLFYjHg8ntUr\n57EYzJ8Pd9wBxx4LV14Jl10GVVVNPqWISOQSiQSJRIJkMkkymczqXLnODAYCte4+KrU9BKhx90u2\n87w5zQzSvfJKmFNwD1nCwQfn5TIiIgVXTC+drQOq07arCdlB0WjXDhIJGDgQjjkGbr8dtmyJOioR\nkWjlejBYDLQ3s5iZ7QQMAmbn+BpZa9YsLDutq4P77oMTT4RXX406KhGR6GSztHQ6UAd0MLM1ZjbC\n3TcBo4F5wApgpruvzE2oudehAzz1FJxyCnTvHuYUirBUk4hI3pV0obpcWrEChg+HPfeESZPggAMK\nenkRkawV05xByerYMTw2Ov54OProMCAU4TgpIpIXygy24oUXYNgwaN0a7rwT2rSJLBQRkYwpM8ix\nzp1D05wePUIDnWnTlCWISHlTZvA5li4NWcJBB4UJ5v32izoiEZGtU2aQR0ceGcpiH3EEdOkCM2Yo\nSxCR8qPMYDssWhSyhMMPhwkTYJ99oo5IROQTygwKpHv3UBq7bdswr/Dgg1FHJCKSG8oMmqiuLtQ4\nOvroUNKiZcuoIxKRSqfMIAK9eoXJ5VatoFMnmF10RTdERDKnzCAHnnwyZAnHHgu33gp77RV1RCJS\niUoqMzCzY83sV2Z2l5k9Xejr58Nxx4UX1XbfPcwlzJ0bdUQiItsnsszAzAYA+7r7XVv5u5LKDNI9\n/njoqnbiiXDLLaHWkYhIIUSSGZjZZDNbb2b1DfbXmtkqM1ttZmO2cYrzgPuaev1idcIJIUuoqgpZ\nwqOPRh2RiMjny+Yx0RRC4/uPmVkVMD61vyMw2MwOM7OhZjbOzFqnPncA8L/u/kEW1y9au+8OEyfC\nXXfB+efDRRfB++9HHZWISOOaPBi4+0LgnQa7ewAvu3vS3TcCM4AB7j7N3b/l7m+kPnc+MLmp1y4V\nX/lKyBL++9+QJSQSUUckIrJ1O+T4fG2ANWnba4Gahh9y92s+70TxeJxYLEYsFiMejxOPx3MWZCG1\naAGTJ8PDD8OQIXDmmXD99fCFL0QdmYiUukQiQSKRIJlMkkwmszpXVhPIZhYD5rh7p9T2QKDW3Uel\ntocANe5+yXaet2QnkLflnXfg0kvhmWdgypSwFFVEJFeKaWnpOqA6bbuakB0I4f2DadPg5pvhnHPg\nO9+BDz+MOioRkdwPBouB9mYWM7OdgEGA3s1t4PTTw1zCunXQtWvIFEREopTN0tLpQB3QwczWmNkI\nd98EjAbmASuAme6+MjehlpcvfjGUw7722jA4jB0L//lP1FGJSKVSOYoisH59WH760kswdSp06xZ1\nRCJSioppzkCaoFUrmDULvv99OOUU+OEPYcOGqKMSkUqiwaBImMF558Hy5bBsWeidsGxZ1FGJSKXQ\nYFBkvvSlUA7729+Gk06CH/8YNm6MOioRKXeaMyhia9fCyJHw1lthLuGII6KOSESKmeYMytT++4dy\n2BddBMcfH95c3rQp6qhEpBwpMygRr70WSmO/9x7ccw8cdljUEYlIsVFmUAEOPBDmz4fhw6F379Ar\nYfPmqKMSkXKhzKAE/fWvoc3m5s2hxlH79lFHJCLFQJlBhTnoIHjiiVDfqGdP+MUvYMuWqKMSkVKm\nzKDErV4dHh3tuGMolX3QQVFHJCJRKdrMwMzamtndZvbA1rYle+3bw5NPQv/+UFMDv/qVsgQR2X4F\nyQzM7AF3P7ux7a18XplBE6xaBcOGhbabkyaFSWcRqRx5zwzMbLKZrTez+gb7a81slZmtNrMxTQlA\ncufQQ+Hpp6Fv31Ds7u67QWOqiGQi08dEUwhN7j9mZlXA+NT+jsBgMzvMzIaa2Tgza53bUCUTO+wQ\nymE/8UR4ZHTKKeFNZhGRbcloMHD3hcA7DXb3AF5296S7bwRmAAPcfZq7f8vd3zCzvc3sDqCrmY1p\nuJ3Tr0Q+5YgjQtOcnj3hqKNCOQtlCSLSmB2yOLYNsCZtey1Qk/4Bd/8n8I0GxzXc3qp4PE4sFiMW\nixGPx4nH41mEWpl23BGuugpOOy3MJcyaBRMnhmJ4IlL6EokEiUSCZDJJMpnM6lzZDAZ5/T0zkUjk\n8/QVpWtXWLQIfvKT8N/jxsHgwaFstoiUroa/KFsW39TZLC1dR2h4/5FqQnYgRWinncJg8PDD8NOf\nwllnwZtvRh2ViBSLbAaDxUB7M4uZ2U7AIGB2bsKSfOnWDZYsgQ4doHNneEBvfIgIGb5nYGbTgT5A\nS+BN4Cp3n2Jm/YBbgSpgkrtfn5Og9J5BQTzzTHh7uWtXGD8evvjFqCMSkWxk856BylFUuA8/hB/8\nAKZPhwkT4PTTo45IRJpKg4FkbeHCUAn1o8J3e+0VdUQisr2KtjaRlI7evWH58jAIdOoUJppFpHIo\nM5DPSCTg/PMhHg/LUPfcM+qIRCQTygwkp+JxeOEFaN48ZAnz5kUdkYjkmzID2ab582HkSDj55NBq\nc/fdo45IRBqjzEDy5qSTQpawZUvIEh5/POqIRCQflBlIxubOhVGjYMAAuPFG2G23qCMSkXTKDKQg\n+vWD+np4/33o0iUsRxWR8qDMQJpk9my46CI45xy47jrYZZeoIxKRos0MttIDOW5mC83sV2bWJ5/X\nlvw67bQwl7B+fShn8ac/RR2RiGQjr4OBu7/q7iPTdm0B3gOaowqnJa9lS7jvvpAZnHEGXHEF/Oc/\nUUclIk1R6B7IC939FGAs8KMmxCtFaODAkCX89a+hq9qiRVFHJCLbq6A9kNMmAt4lZAdSJvbdN5TD\nvuoqOPXUUPxuw4aooxKRTBWyB/JYMzsjtX0vcHsuvxCJnhmce26ocVRfH3onLF0adVQikokoeiA/\nlMnJ1QO5dO23H/z2t/DrX4c3l7/5TbjyytCTWURyJ5c9kDNeWmpmMWCOu3dKbQ8Eat19VGp7CFDj\n7pdkFRFaWlpO1q2Dr38d/vY3mDo1vMUsIvkR1dJS9UCWz9WmDfz+9zB6NJxwQlh5tGlT1FGJSEPq\ngSx5ZxZKYi9ZEspj9+oFK1ZEHZWIpMt0ael0oA7oYGZrzGyEu28CRgPzgBXATHdfmb9QpdQdcEAo\nh33BBdCnD9x8M2zeHHVUIgIqRyERefXVkC38979wzz3QoUPUEYmUvqItRyHSmLZt4bHH4LzzwmOj\nW28NZbJFJBrKDCRyL78Mw4dDs2YwZQq0axd1RCKlSZmBlLSDD4YFC0J9o5oa+OUvlSWIFJoyAykq\nL70Ew4bBrrvC5MkQi0UdkUjpUGYgZeOQQ+Dpp6G2Frp3hzvvBP1eIJJ/ygykaL34YsgSWraEu++G\n6urPP0akkikzkLJ0+OGhaU7v3qE09pQpyhJE8kWZgZSE5ctDltCmDdx1F7T+TIF0EVFmIGWvSxd4\n7rlQFrtr11ARVb8viOSOMgMpOUuWhPcS2rWDiROhVauoIxIpDkWbGZhZWzO728weSG13NLOZZjYh\nVQJbZLsdfTQsXgwdO4aM4f77o45IpPQVJDMwswfc/Wwz+zbwnLs/ZWa/c/cBjXxemYFk5LnnwlxC\np07hZbV99ok6IpHo5D0zMLPJZrbezOob7K81s1VmttrMxmRwqmnAuWZ2E9CyCfGKfEqPHvD883Dg\ngdC5MzyUUR89EWkoo8zAzHoD7wP3pnU6qwJeAvoSGt0sAgYD3YCjgJvd/Y3UZx9w97PTzlcFzHL3\n0xu5njID2W51dWEuoXt3uP122HvvqCMSKay8ZwbuvhB4p8HuHsDL7p50943ADGCAu09z92+5+xtm\ntreZ3QF0NbMxZnagmU0EpgI3NSVgkcb06gXLloVHRZ06wZw5UUckUjp2yOLYNsCatO21QE36B9z9\nn8A3Ghx3YSYnj8fjxGIxYrEY8XiceDyeRahSKXbdNZTDPuOM0C9h1qyw3aJF1JGJ5F4ikSCRSJBM\nJkkmk1mdK+MJZDOLAXPSHhMNBGrdfVRqewhQ4+6XZBURekwkufH++3DFFSFDuOuuUO9IpJxFtbR0\nHZBeLaaakB2IFIXddoMJE0IZiwsvhFGj4F//ijoqkeKUzWCwGGhvZjEz2wkYBMzOTVgiudO3L9TX\ng1lYcfTYY1FHJFJ8Ml1aOh2oAzqY2RozG+Hum4DRwDxgBTDT3VfmL1SRpttjj1AOe+LEsOLo4ovD\nYyQRCVSOQirOu+/C5ZfDwoWhgU6fPlFHJJIb2cwZaDCQivX734e5hLPPhuuuCyuRREpZ0dYmEilm\np54a5hL+8Y9QCbWuLuqIRKKjzECEUMbim9+Er30Nfvxj2GWXqCMS2X7KDESydMYZoYHOa6+FrmrP\nPht1RCKFpcxApIH774dLLw1vMF99NTRvHnVEIplRZiCSQ+ecE7KElStDZ7UlS6KOSCT/NBiIbEWr\nVvDggzBmDPTrFzKEDRuijkokfzQYiDTCDIYMCZVQFy8OvROWL486KpH80GAg8jlatw7vJFx+eSht\nce21sHFj1FGJ5JYmkEW2w5o1MHIkvP023Htv6MMsUiyKdgLZzAaY2Z1mNsPMTjKzQ83sV2Z2v5ld\nkM9ri+RDdTX88Y/w9a+HMhY33QSbN0cdlUj2CpIZmFkL4GfuPjK13QyY4e7nNPJ5ZQZS9JLJsPz0\nww/hnnvgkEOijkgqXd4zAzObbGbrzay+wf5aM1tlZqvNbMw2TvEDYHzqmP7Aw4Q2mSIlKxaDRx8N\nk8xf/jL8/OfKEqR0ZZQZmFlv4H3g3rROZ1XAS0BfQqObRcBgoBtwFHAz8DfgBuARd3+swTl/5+4D\nGrmeMgMpKa+8AiNGgHtopnPwwVFHJJUo75mBuy8E3mmwuwfwsrsn3X0j4Tf9Ae4+zd2/5e5vAJcA\nJwJnmdmFZtbHzG4zs4nAE00JWKQYtWsHiQQMHAjHHAO33w5btkQdlUjmdsji2DbAmrTttUBN+gfc\n/RfALxoctyCTk8fjcWKxGLFYjHg8TjwezyJUkfxr1iwsPz3llNBA58EHQ7+Etm2jjkzKVSKRIJFI\nkEwmSSaTWZ0rm8Egr89xEolEPk8vkjcdOoTGOePGhRfVrr02rD6yJiXvIo1r+IuyZfGPLJulpeuA\n6rTtakJ2IFLxqqrgu9+FBQtg0iT4ylfg9dejjkqkcdkMBouB9mYWM7OdgEHA7NyEJVIeOnYMTXOO\nPx6OPjo8NtLaCClGma4mmg70AVoCbwJXufsUM+sH3ApUAZPc/fqcBKXVRFKG6uth2DDYbz+46y5o\n0ybqiKTcqAeySInYuBGuvx7Gj4dbbgnvKGguQXJFg4FIiVm6NGQJbdvCxIkhWxDJVtHWJhKRrTvy\nyFAWu1Mn6NIFZszQXIJES5mBSMQWLQpZQseOMGEC7Ltv1BFJqVJmIFLCuneH558PbzF36QKzZkUd\nkVQiZQYiReRPfwpvLx91VJhkbtky6oiklCgzECkTPXuGNptf+lKYT5itN3ekQJQZiBSphQtDJdRe\nveC222CvvaKOSIqdMgORMtS7NyxfDnvsEbKEuXOjjkjKmTIDkRLw+OOhq1rfvuFltT33jDoiKUbK\nDETK3AknwAsvhAJ4nTvD/PlRRyTlJq+ZgZkNAL4K7AFMAj4EvkYond3R3b/cyHHKDEQaMW8ejBoV\n+ibcfDPsvnvUEUmxKPpyFGbWAviZu49MbQ8A9nX3uxr5vAYDkW1491349rfhiSdCJdTjj486IikG\neX9MZGaTzWy9mdU32F9rZqvMbLWZjdnGKX4AjE/bPg+4b/vDFRGAFi3CIDB+PAwdCpdeCh98EHVU\nUsoynTOYAtSm7zCzKsIP+FqgIzDYzA4zs6FmNs7MWltwIzDX3ZeljjsA+F931z9dkSx99athLuGd\nd6BrV3jqqagjklKV0WDg7guBdxrs7gG87O5Jd98IzAAGuPs0d/+Wu78BXAKcCJxlZhemjjsfmJyb\n8EVk771h2rQwf3DOOaHD2ocfRh2VlJpseiC3Adakba8FatI/4O6/AH7RYN81mZw8Ho8Ti8WIxWKf\n6fMpIp91+ulw7LEwenSoijp1KtTUfP5xUroSiQSJRIJkMkkymczqXBlPIJtZDJjj7p1S2wOBWncf\nldoeAtS4+yVZRYQmkEWy9cADcMkloc7Rj34EzZtHHZEUQlTvGawDqtO2qwnZgYhE7Oyzw1zC6tWh\n9/LixVFHJMUum8FgMdDezGJmthMwCFBZLZEise++8JvfwPe/Hyaaf/hD2LAh6qikWGW6tHQ6UAd0\nMLM1ZjbC3TcBo4F5wApgpruvzF+oIrK9zGDw4FAJdfny0Dth2bKoo5JipNpEIhXCPaw6+u53wyTz\n974HO+4YdVSSS0X/BvL20mAgkj/r1sHIkfDmm2HF0RFHRB2R5IoK1YlIxtq0gT/8AS6+OJSxuP56\n2LQp6qgkasoMRCrY66/DBRfAv/4F99wDhx0WdUSSDWUGItIkBxwAjzwSOqoddxz87GeweXPUUUkU\nlBmICAB//WtooLNxY8gS2rePOiLZXsoMRCRrBx0UOqqdey707Bn6Lm/ZEnVUUijKDETkM1avDo+O\nqqpgypQwUEjxU2YgIjnVvj0sWACnnQY9esCECcoSyp0yAxHZplWrYNiw0F5z0iQ48MCoI5LGKDMQ\nkbw59FB4+mk48UTo1g3uvju8zSzlJa+ZgZkdClwGtCTUMHoc+D6wp7ufvY3jlBmIFKE//zlkCfvs\nEwaF/fePOiJJV7SZgbuvcveLgHOBk939VXcfmc9rikj+HHEEPPMM9OoFRx0Vylno97bykGnV0slm\ntt7M6hvsrzWzVWa22szGNHJsf+BhQltMESlxO+4IV10VXlb7+c/DJPPf/hZ1VJKtTDODKYTG9x8z\nsypgfGp/R2CwmR1mZkPNbJyZtQZw9znu3g8YlsO4RSRiXbvCokWhxWbXrvA//6MsoZRl0/ayJ3C1\nu9emtscCuPsNacf0Ac4EdgZWAvcC1wEnAne7+42NXEtzBiIlZPHiMJdwyCFwxx2hsY4UXjZzBjtk\ncd02wJq07bXAp9pvu/sCYEGD476Rycnj8TixWIxYLEY8Hicej2cRqojkU7dusGQJXHMNdO4M48fD\nWWdFHVX5SyQSJBIJkskkyWQyq3NlkxkMBGrdfVRqewhQ4+6XZBURygxEStkzz8Dw4eHR0S9/CS1b\nRh1R5YhqNdE6oDptu5qQHYhIBTvmGFi6NPRN6NQJfve7qCOSTGQzGCwG2ptZzMx2AgYBs3MTloiU\nsl12gVtugfvvh+98B4YOhXfeiToq2ZZMl5ZOB+qADma2xsxGuPsmYDThZbIVwEx3X5m/UEWk1Bx7\nLCxfDnvtFbKEhx+OOiJpjGoTiUhBJBKhX0I8DuPGwZ57Rh1R+SnaN5BFRD4Sj4csoXnzkCU88kjU\nEUk6ZQYiUnDz58PIkXDyyWFuYffdo46oPCgzEJGSctJJUF8f3lju1Cl0WJNoKTMQkUjNnQtf/zoM\nGAA33AC77RZ1RKVLmYGIlKx+/eCFF+C998KLagsXRh1RZVJmICJFY/ZsuOgiOOcc+OlPYdddo46o\ntCgzEJGycNppIUtYvz5kCXV1UUdUOZQZiEhRmjULRo8Oby//+Mew885RR1T8lBmISNkZODBkCa++\nGrqqLVoUdUTlrdA9kF8GrgX+DMxIlbje2nHKDEQECMtP778fLr00vJtw1VXhxTX5rKLNDBr2QAYc\neA9ojiqcikgGzGDQoPD28p//HHonPP981FGVn0L3QF7o7qcAY4EfZRm7iFSQ/faD3/4WrrgCamtD\nI50NG6KOqnwUtAdy2rOfdwnZgYhIxszChPLSpfDcc1BTE+YVJHuF7oH8GuFxUQtggrs/2ci1NGcg\nItvkDpMnw9ixcPnlMGYM7JBNI98yUGo9kB/K5OTqgSwi22IGF1wQ6hxdcEF4hDR1KnTsGHVkhZPL\nHsjZDAZ5/dU9kUjk8/QiUiYOOCCUw544EY47LswpfOc7UFUVdWT51/AXZbMmJQWAeiCLSBkwg298\nI7yLMHcu9O4Nf/lL1FGVFvVAFpGy0bYtPPYYnHce9OoFt94KW7ZEHVVpUA9kESkrzZqFMhbPPAO/\n+U3osPbKK1FHVfxUm0hEytbmzXDbbXDddfCjH4WKqM3KuAhPNquJNBiISNl76SUYPjyUxJ40CWKx\nqCPKj6ItRyEiUgwOOQSeeir0XO7eHe68M7ynIJ9QZiAiFeXFF0OWsPfecPfdUF39uYeUDGUGIiIZ\nOvzw0DSnd+9QGnvKFGUJoMxARCrY8uUhS2jTJjw6at066oiyo8xARKQJunSBZ58NZbG7doVf/7py\nswRlBiIihB4Jw4ZBu3ahtEWrVlFHtP2UGYiIZOmoo2Dx4jCn0KULzJwZdUSFpcxARKSB554LcwmH\nHw4TJsA++0QdUWaUGYiI5FCPHuGxUdu20LkzPPhg1BHlX14zAzM7FLgMaEmoYfQn4GrgbeAxd5/V\nyHHKDESkKNTVhSyhWze4/XZo2TLqiBpXtJmBu69y94uAcwkdzmqB2939YuD/5fPaIiK50KsXLFsW\nJpQ7d4Y5c6KOKD8yrVo62czWm1l9g/21ZrbKzFab2ZhGju0PPAxMB6YB55rZTYRsQUSk6O26K4wb\nB9Onhxabw4fDu+9GHVVuZZoZTCH8Vv8xM6sCxqf2dwQGm9lhZjbUzMaZWWsAd5/j7v2A4e7+lruP\nBr4H/CNnX4WISAEcd1x4UW3XXaFTp9BIp1xkPGdgZjFgjrt3Sm33BK5299rU9lgAd78h7Zg+wJnA\nzsBKQv/jK4EvABPcva6Ra2nOQESK2qOPftKD+ec/hz32iDqi7OYMsumB3AZYk7a9FqhJ/4C7LwAW\nNDjuwkxOHo/HicVixGKxz/T5FBGJWt++UF8f+i136hRKY/ftW9gYEokEiUSCZDJJMpnM6lzZZAYD\ngVp3H5XaHgLUuPslWUWEMgMRKS3z5sHIkdC/P9x0E+y2WzRxRLWaaB2QXvy1mpAdiIhUlJNPDlnC\nhx+GFUeJRNQRbb9sBoPFQHszi5nZTsAgYHZuwhIRKS0tWoRy2LfdBl/7Glx2GXzwQdRRZS7TpaXT\ngTqgg5mtMbMR7r4JGE14mWwFMNPdV+YvVBGR4te/f8gS3n47VEJ9+umoI8qMahOJiOTJQw/BxReH\nTOEnP4Fddsnv9Yr2DWQRkUp2xhnwwgvw+uuhKuqzz0YdUeOUGYiIFMDMmWEe4fzz4eqroXnz3F9D\nmYGISJEbNCi8vbxyZSh69/zzUUf0aRoMREQKpFWrUA577Fjo1y9kCBs2RB1VoMFARKSAzMKE8tKl\nsGQJ1NSEjCFqGgxERCLQunUoh33ZZaGMxbXXwsaN0cWjCWQRkYitWRPKWbz9NkydGtptNoUmkEVE\nSlh1Nfzxj3DhhRCPw403wubNhY1BmYGISBF57bWw/PTf/4Z77oFDDsn8WGUGIiJl4sADYf58GDoU\nvvzl0CuhEFlC3jMDM/sCkACuAV4BLiO0vJzn7pMaOUaZgYhUvFdegREjwD0UwTv44G1/vtgzgyuA\nmQDuvsrdLwLOBU4uwLVLXqIUa+Hmie7FJ3QvPlHO96Jdu1AO+6yz4Jhj4PbbYcuW/Fwr06qlk81s\nvZnVN9hfa2arzGy1mY3ZynEnESqavpW2rz/wMDAju9ArQzn/Q99euhef0L34RLnfi2bNwvLTujqY\nMQNOOAFefTUP18nwc1MIje8/ZmZVwPjU/o7AYDM7zMyGmtk4M2sN9AGOAc4DRll4/jPH3fsBw3L2\nVYiIlLkOHeDJJ+HUU6FHD7jjjvD4KFcy6oHs7gtTbS/T9QBedvckgJnNAAa4+w3AtNRnfpD6u2GE\n7OA4MzsT2Bl4ItvgRUQqSVUVfPe78NWvwrBh8PvfhxfXrEmzBJ+WTQ/ks4CT89UDOdtziIhUoqZO\nIGeUGTR2zSyO3faJm/jFiIhI02SzmmgdUJ22XQ2szS4cERGJQjaDwWKgvZnFzGwnYBAwOzdhiYhI\nIWW6tHQ6UAd0MLM1ZjbC3TcBo4F5hOWjM919ZTbBfN5S1XJmZtVm9oSZvWhmfzazS1P79zaz+Wb2\nFzN7xMzjJ7PbAAADRklEQVRaRB1roZhZlZktNbM5qe2KvBdm1sLMfmNmK81shZnVVPC9+F7qe6Te\nzO4zs+aVci+2tsR/W1976l6tTv1M/crnnr9Y3vRNLVV9CehLeAS1CBic7QBTKsxsP2A/d19mZrsB\nS4DTgRHAP9z9ptQAuZe7j40y1kIxs28DRwO7u/tpZnYTFXgvzGwqsMDdJ5vZDsAXgO9TYfcitYjl\nceAwd/+vmc0E/gAcTgXcCzPrDbwP3Ju2kGer3xNm1hG4D+gOtAEeBTq4e6OvrBVTbaKPl6q6+0bC\nS2kDIo6pYNz97+6+LPXf7wMrCf8TTwOmpj42lTBAlD0z2x84Bbgb+GhBQcXdCzPbE+jt7pMB3H2T\nu/8vFXgvgH8BG4FdU4PirsAbVMi9cPeFwDsNdjf2tQ8Aprv7xtTy/5cJP2MbVUyDQRtgTdr22tS+\nipP6DehI4FmglbuvT/3VeqBVRGEV2jjg/wPpv8lU4r1oC7xlZlPM7HkzuytV76vi7oW7/xO4BXid\nMAi86+7zqcB7kaaxr701n17Q87k/T4tpMCiO51URSz0imgVc5u7vpf9dqnpf2d8nMzsVeNPdl/JJ\nVvAplXIvCMu/jwImuPtRwAfApx6BVMq9MLN2wOVAjPDDbrfU+00fq5R7sTUZfO3bvC/FNBhU/FJV\nM9uRMBBMc/ffpnavT80nYGZfAt6MKr4C6gWcZmavAtOBE8xsGpV5L9YCa919UWr7N4TB4e8VeC+6\nAXXu/nZqAcuDQE8q8158pLHviYY/T/dP7WtUMQ0GFb1U1cwMmASscPdb0/5qNp/UcRoG/LbhseXG\n3a9092p3b0uocPu4uw+lMu/F34E1ZtYhtasv8CIwhwq7F8Aq4Bgz2yX1/dKXsJKxEu/FRxr7npgN\nnGtmO5lZW6A98Nw2z+TuRfMH6EdYUfQy8L2o4ynw134s4fn4MmBp6k8tsDdhJcBfgEeAFlHHWuD7\n0geYnfrvirwXQBfC6rrlhN+G96zge3EFYTCsJ0yY7lgp94KQJb8BbCDMr47Y1tcOXJn6WbqKUDpo\nm+cvmqWlIiISnWJ6TCQiIhHRYCAiIhoMREREg4GIiKDBQERE0GAgIiJoMBARETQYiIgI8H/0Vphc\nVKbqVAAAAABJRU5ErkJggg==\n",
              "text": [
                "<matplotlib.figure.Figure at 0x7fd353203c18>"
              ]
            }
          ],
          "prompt_number": 4
        },
        {
          "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",
          "collapsed": false,
          "input": [
            "k = 10\n",
            "p = 5"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 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",
          "collapsed": false,
          "input": [
            "Omega = np.random.randn(n, k+p)\n",
            "Y = A.dot(Omega)\n",
            "Q, _ = la.qr(Y)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 6
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "As an alternative to the above, use a few iterations of the power method on `Omega`:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "Omega = np.random.randn(n, k+p)\n",
            "Y = A.dot(Omega)\n",
            "Y = A.dot(A.T.dot(Y))\n",
            "Q, _ = la.qr(Y)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 7
        },
        {
          "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",
          "collapsed": false,
          "input": [
            "C = Q.T.dot(A)\n",
            "print(C.shape)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "output_type": "stream",
              "stream": "stdout",
              "text": [
                "(15, 100)\n"
              ]
            }
          ],
          "prompt_number": 8
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Sanity-check that form:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "Aapprox1 = Q.dot(C)\n",
            "la.norm(A - Aapprox1, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 9,
              "text": [
                "3.9724863710473977e-06"
              ]
            }
          ],
          "prompt_number": 9
        },
        {
          "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",
          "collapsed": false,
          "input": [
            "UC, sigmaC, VTC = la.svd(C, full_matrices=False)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 10
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Reconstruct the SVD of $A$: $A\\approx QU_C \\Sigma_C V_C^T$"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "UAapprox = Q.dot(UC)\n",
            "sigmaAapprox = sigmaC\n",
            "VTAapprox = VTC"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [],
          "prompt_number": 11
        },
        {
          "cell_type": "markdown",
          "metadata": {},
          "source": [
            "Compare the 2-norm of the reconstructed $A$ with the original $A$:"
          ]
        },
        {
          "cell_type": "code",
          "collapsed": false,
          "input": [
            "Aapprox = (UAapprox*sigmaAapprox).dot(VTAapprox)\n",
            "la.norm(A - Aapprox, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 12,
              "text": [
                "3.972486371045497e-06"
              ]
            }
          ],
          "prompt_number": 12
        },
        {
          "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",
          "collapsed": false,
          "input": [
            "omega = np.random.randn(n)\n",
            "Aomega = A.dot(omega)\n",
            "err = Aomega - Q.dot(Q.T.dot(Aomega))\n",
            "la.norm(err, 2) / la.norm(omega, 2)"
          ],
          "language": "python",
          "metadata": {},
          "outputs": [
            {
              "metadata": {},
              "output_type": "pyout",
              "prompt_number": 18,
              "text": [
                "3.3931213029085818e-07"
              ]
            }
          ],
          "prompt_number": 18
        },
        {
          "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?"
          ]
        }
      ],
      "metadata": {}
    }
  ]
}