{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Interpolative Decomposition"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import scipy.linalg as la\n",
        "import matplotlib.pyplot as pt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Obtain a low-rank matrix"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "2.3513490975788167e-13\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f863620a210>]"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "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": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "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) @ VT0 - A0))\n",
        "\n",
        "sigma = np.exp(-np.arange(n))\n",
        "\n",
        "A = (U0 * sigma).dot(VT0)\n",
        "pt.semilogy(sigma)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run the factorization"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import scipy.linalg.interpolative as sli"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute a fixed-rank factorization:\n",
        "\n",
        "(There's also an adaptive, fixed-precision mode.)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "k = 20\n",
        "idx, proj = sli.interp_decomp(A, k)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Examine `idx`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([49, 97, 46, 22, 18, 35, 13, 88, 99, 83, 59, 61, 71,  3, 55, 67, 73,\n",
              "        0, 40,  2, 20, 21,  6, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,\n",
              "       34,  5, 36, 37, 38, 39,  4, 41, 42, 43, 44, 45, 19, 47, 48, 17, 50,\n",
              "       51, 52, 53, 54, 14, 56, 57, 58, 10, 60, 11, 62, 63, 64, 65, 66, 15,\n",
              "       68, 69, 70, 12, 72, 16, 74, 75, 76, 77, 78, 79, 80, 81, 82,  9, 84,\n",
              "       85, 86, 87,  7, 89, 90, 91, 92, 93, 94, 95, 96,  1, 98,  8],\n",
              "      dtype=int32)"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "idx"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What does `numpy.argsort` do?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sort_idx = np.argsort(idx)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
              "       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,\n",
              "       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,\n",
              "       51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,\n",
              "       68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,\n",
              "       85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99],\n",
              "      dtype=int32)"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "idx[sort_idx]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reconstruct the matrix:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "B = A[:,idx[:k]]\n",
        "P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]\n",
        "Aapprox = B@P"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "7.211609036219504e-09"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(A - Aapprox, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What's the structure of $P$?\n",
        "\n",
        "(ignoring the column permuation)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7f8634713e10>"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABlCAYAAABdnhjZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXl8VeW57593z5lnEjKQACEgg8wgg6AijgjWczzF2pbTo/ba4aOtxx61p/fctre9p6ftdei9tqc9rfdqa6XWWZxARETmIcxDCCEJgYRMZM5OdrLf80fC+j3PPiIIGHTn+f7Du7Oevda73vXuxXp/6xmMtZYURVGUzz+uS90BRVEU5eKgN3RFUZQoQW/oiqIoUYLe0BVFUaIEvaEriqJECXpDVxRFiRL0hq4oihIlXNAN3RhzgzHmkDGm1Bjz8MXqlKIoivLJMecbWGSMcRNRCREtJKIqItpKRHdYa/dfvO4piqIo54rnAr47g4hKrbVlRETGmOVEtISIznhDT0112by8vkMebB5yAYeOAnqN0zQR/6eKj74w2j0RCyoXLP3+kNMO9brl/jrw2frwHZcnLOyoHXbhGNkp44Ktm3XY6+4Vdt29mFLhMM7RsvMlIjI97DPrxvD0WmF3tAHzxNWNv/vTuoRdR9CHfXefeWzdCRin3lbvR/6diMjDzjfY5hfbYuJx7GAztoX98mDuDvTDk4LOd7f5hB2/4Pz6UFiOmcffA7tmjLMrWfbdxU66q8tLZ4JfA8unVuS84NfOyj65O5lZPI7r8cp5QY3ob372SaddVp8pzGwAxzZBdMqX0C3sQs0YQxuH79iIMeNzwcWGqTdifruYXTiG9aEr4jfH+mdDchvfP58Lxh0xL9jc6unEuPB5RSR/S52l1fXW2gw6CxdyQ88homPscxURzfy4L+TleeidN9OJiGjm2/ddwKE//3hOsR+k/D1SmN2Pw7lBfKiXN5ZwDH40RYXVTruqKUna7cbnYB5+GLEpncLObIFdx7ig2BaIxffiYzDxsuNbhF1FUwqO1Y2bCb/xERF567HN04of0zN3Py7svvTs/U47oRx/L1h2WNhtLylw2jFH8WN3y/sApSzAODWsHeq00+ZXC7vkAMbm0IfDxbYJ83DsA28VOe2OkfJgqVtwjhlfrHTaFR/ky06x3zu/Pq5W+fNMHdXotHv6f0dERAlLZN/9Htz4S4+yG2bEf6q+euw/xG7GJlXeWMLtOA8TkvtI3YmbWsMcTOSsoafkPv6M/5j/438+5rSX/uEBYdc1GuPuPxjjtPOurhR2J94ahu9Ma3PaPUH5H5j3GOZCbA363jShR9jFVmAsgmPRB09FQNjRqHanGaqNEZtiqvHD7SjAWHgT5XgmxeO31bQvzWlPmFUq7I6ewrbdi39SQefAhWjo5iP+9l/0G2PM140x24wx2xoawh/xFUVRFOVicCEa+iwi+qG19vr+z48QEVlr//VM3/EX5Nqs/973ZL75BvkkNvPt75xXPz6vGB+erlNS28S21Fg8IbR24cm2aauUqeJO4No1jUY7qVA+HbW0xjrtGQXlTnvHu5cJu69+YbXTXvG/rhbbahfjqcJVgSeT7Kny6bBhZbbTTi7DOZ68XT7xj8854bSPteCp3vvHVGFXtwTfm8n6frxdrkIqa/C9mH3oX1eKnN+x1XgOSb35OPrdHivseregT6OuPyK27TqEp8PYtA6nHaxIEHajJ+OpsvH3+E7LcPkcFVOLPjZfhWvfG5RP6H72tDn2aqwSdpYNE3YpG2HXMhJ/NxHPU3YYjlU879+d9pxtXxN2rZWJTjt3tJTE6lvjnHZXZbzTDtTLc3TPwJxsrYOdt0GeY08COhlThSdeX6vse9N4PGEH0nAeoco4YReoQz8yr61y2hU7s4VdbxqeqGNK8ZsLJUZIj8NwvXt75TlymYlLoNOzjgm7/b8a77RrFuA80tfL1UXqnfjeu9c8vt1aO43OwoU8oW8lolHGmOHGGB8RLSWi1y5gf4qiKMoFcN4aurW2xxjzbSJ6h4jcRPSUtXbfReuZoiiK8om4kJeiZK19k4jevEh9URRFUS6A89bQzweuoUfCNfXBoKdzd6h/ufZlse3Hmxc5bVct9LzkQ3If3g5cu5qroN/dO+t9YffGCWh29R/Cs8NMahZ2U7Oh2cW4pevNmncnOW3u8pVc2CjsEvzw0qjeimN5W+Q79I489NcVxLaM7cKMGsczzXsfjtswXu7POxreNl3l0LIjvOwoex320bIM3wn1SFfP8Vl4N1C8oUhsSz6AdtNotH3NEedYBO8Gtx/nG45wd4s5CE+KzizmFueXonfBy/jc9E28d2nfnyLsPG3MBY95SHZlntmzw8eclXqk8wZ15DA3viHSY4P7hbrK8cWeHGk3ZBU6UjcFf7feCJc+5o7ay7y6EoZJb6oYH+Zn50q8W2rPidjfMHil5D+BcamfKN+ZWDdzdb2uwWmPSZPvDNp60Kd9lUPFNlc1rqO/EP3lbopERFfn4f3Hvibso6Yl4Yx2T0597lPX0BVFUZTPEHpDVxRFiRIuSEO/mHCZZTDILzYOy99/23O92OY/iqVbcCiWlqfGS1kgdReWiTPHwbWu0H9S2DW/udBph+fA/ysckvtbVzzGaZs4uTx3sQhGyyLfuEskEVHHNgS89CRhqRnKlUvwkbl1Tjv4JFzIxjy0V9itf22i007ZBXnnyz/YIOx++f6NTtvLov6sPEVqXgaZqa0Cro+5Y+SYbT8KV8D4E1JK6bkV/bClyejffLmPjuNwpQzX4ZpGSim+K7C/7oOQT3oSZbRl1QLIFq7dsMuYLo9bfQQBhZYFnyUVywjVEFvht+WzwKK8DmHnLYErYE+vDBBLY3OwiSlTmUOknOcOYl48c+tvnPZdf/qWsLtyBmSGzdtx7R+5/m1h9/iRBeh7AcYzbVSDsGvcj+NWXc2uY4QU153E5ncz5vTG46OEXf5ISDCpKe1iW2ZejdPedzDPaa+44Qlht6L1cqe97vnpTvvKu4uF3Z5T0rXyXNAndEVRlChBb+iKoihRwmfGy+VMRGtEafpQLEnrj8uoR38NIsZ6R2H56z0g5Y3gEBZtWoBIvMYauT9vPDxPetjSf/hYGeVZxvJ+RCYlcrfhc2RiI056IZa83HPk5nwZorB81Vyn7crFOZoyeY5Dp6OPlSXoX/xRqaW0jmeSThDbYqqkqthZyOyYt4kJRCSTYiQUy3we/kacf91VGFtvQMpU3AMk1AqpwhVh99C0d5z2v66Fh5NIzEZEnjrMi55U7MNXK88xlMw9ZXBe3joZiegaySKUDyJ6sytP5qQZNxxRvSUbCsS23liWqI15MhXMlblX8uKanPbqTRPQh4gh4xIZ91CyEWMxejSifCvWIjdOsCDCC6eL7ZDtIzZJ5jGKC+CcW9pxvb0RScbam+HJM7tIRhBvPIKcP+HOM6vZ8RmQatpbcKxJw2VEafE+7K/y3n9SLxdFUZTBhN7QFUVRogS9oSuKokQJnxm3xTMRqZlHi0sj1+kiSZ8J96fm97KcdtvwCMGRaYIJ/w7dvKdIaqXDb4WeuasJemPZ4Sxh98hVK5z2h02FYtvGdePwgbktBk5KLbu5DW5i3Rno77ONM4Td9FklTntZ1nqn/eM3Zaa/qxfBrioJro479lwu7BJSoMMPS4ZeO3SKjDBcfQCumUWF0IYPlUdE/bFsmK2FUkf17GLaO8tZ3uOS7xb8sazoCCsm8dXLNwu7J/Yjs6WnBeNphskMlT1x2PbNWe857Wd/J91e3ewdQnAkrkFktG5iAnNPPAjXxJpsabd/N1w4p80rEdt+nY98fPN//T2nfaRG1mLIG4lrMur+TU675HfTpV1+vdOuqobb55D3pMvlIRdc+jJnwpXQrJcZSZff86jTfr4Jx3q9fLywS4+Frl1XBVfU78+XUdw/3LTEaW+pkHnt7Sn0MZFFtgb3JUu7dIzv5QV4F3DgPeki6Q188veb+oSuKIoSJegNXVEUJUr4zLstfhyfZ/klcByyiGuCjKrrbIeLG0+a33tSZk2yHiZ9sPJXwSwpEQRqsI277ZkOqbh5WpiUEOHFZwuxJPXthItb0lU1wu5EBcpm+eqw/9iIaMvWOXAbC7PSaPGJ0p2siy1Xe1mdxhmzZKayTYdHOO0Aq804J++osPugHBUf4mMhafg88oTrd2PpPnpmudhW2YQ+WeZb19YoXS5drO5nOIWVJIuRic+y/4jrXfcPkEFcG6X7addUXIOeRhZ5GlGzcnQRlvFlmyGXeEbJKhGB1QgVvfpuyEAvbZXecTyJF084RkSUsRYyw/z7IKW88u4Vws7mYKwzUiFHNK+XNUW70iEjhpMwThNHVAm7XaWIxDRefGdYlkwWV9OEc+yuwfXxN0qpMJjNrgnbX+CojIztLmLl6Y5K2TTQgLnQMgZS171z1gi7p1671mlz18xQVmQtSjQr73pY3RYVRVEGE3pDVxRFiRI+15IL5/MWUcoll8VLZKKpv+5Bwmg/W56HWTQfEdH0BUjMveOtsU67J05e016WWOtXi/+/036w+HZhV/BTyA4HvyVrM/KK8a44thTOl0vh4hK8+S8cjqRRtSvyhF17LtaTvanY3y0Tdgu7NcvhmcDzcicels8iLUVsfcpyaruOy2XxrPmIWN28Gp47w2fLyEbOuCQZUbuhFhF8TRshGQTz5JI5eRfzNloIKWBKphyzCQn4/OoJJKRq6pASm5fJQvVH4QHiz5TJtHqOYJ7kToUnT1WxTPZkcyEfJK3FsToXSmmmsxHbYiqlB1VoLIvyZfnQfWOkd1HcK6hLapbCW6mxRcpU84cj+nJqQrnT3tA8UtjteAVeKh3ZLF97mpSEYmIRARr7KvqQ9FV5DapXYn4Gh2B/7qCUCrvTcA1MTITXWRPkJ087vhdKk3Y3T8EcX/PSVOx7gryOsVswNnsfe0AlF0VRlMGE3tAVRVGiBL2hK4qiRAmf+UjRc+XzFlHazTLibWSaLBGR7cb/s8EgNGBXgtTGq9rgPhfMgLZnE6RmZ1qge/6i7AZ8p01G3x38Bo7rT5JRiuFebAvEQJfcvUVqm/4OaIeNmdBUWybIDH48c6CbRby+cUBG8FmWRdE04Tvtw+RYTJ5aij5V5TjtJddtEnYvr4Y7XZiNZ8mBXHlcVhii8oiMCOzMwfiaQub66Jfj3p2A/vp57U0j+76jBa6F5WVwl8zOl8UaRiYhivIg+3tTq9TaDXudUH4EGr93mCzIEGD9bRqNeRa3MVHYmWHYYbBIzosvjIYevGo7xvae0R8Ku1+Nvclp+zrhChg+LjX0Y0Mwp9/dg6Ilpl3equwYzIvMTLj91jXIupzhMHMlvBnZJesqZUSpeyL06/hN6FNLUYT/LhvbQJyc093M3XHEBLwnKKmUEdkudoG6kzEXelrl+4mY61g908fonNAndEVRlChBb+iKoihRQtS4LX4cn0X5xdOIJWRSqdzWMANL4cT9WIa1jpDLPx9zV3PvxFKzY5hc+qdkY0ka+gCRnD1ytUvZ8+DKVVYql4mJh9DfEPted4osPODNx7K+twfPCzyykYjI14BtXUPR3/jDctmZdi3c7o7Voo4m1coIvpxxcJFs6jxz4rPwJuwjmUW51rdIN03vVoxn5xD5GwnH4JxvmrnTab+xQyYMM+z8idVoTdom+97FAkKT5uI8ao6mCTv++JVdAPmlrkm6s3IZw5OH69F1KmJcwhGFNU8fJqK4iU2GO6btld8xHlbPMw2SRmJASjP8mrQcwHnFjz4l7Lq2wR0zaTbGonGblEhcrB+xUzEWrfvkmCWMg2w1KQNz6b29Y4SdOxbXp5cVxYiUeiiezdU98jq2jsY4FY7E3KpeJV12O5hc9JXJkARfena+sGsrwLEqv6EFLhRFUQYVZ72hG2OeMsbUGmP2sr+lGmNWGWMO9/+b8nH7UBRFUT59ziq5GGPmEVEbET1jrR3f/7efE1GjtfZnxpiHiSjFWvvQ2Q52qSQXzmcmopQnpCqTyzouhYy8Csml4r0yCm7LkQKn/ee5/+G0v7n3S8LO8yKWsadYWnN/oYzmS/grZAZ3SM6LxqVYunfWo4OFo2QUZWYMogzX72c51SMSSA17Ecvaqi9iqfrTGa8Iu188utRph2/C8jyyRulbv0aN0tYC/D1ri5Spjv8Nkw8asWSOqZHPNh35WO5mRXib/GgUcoB/Y9OXsb9m6TVkWZ1SfxW2deVK74hAJbZFJkXjjLgWc6HuqQKn3byoTdglxkHuSPBjzszNkDUw//IKlvhv/f3PnfbC9d8Wdi4Xq1FaIaWpghU41rjH9zjtd/8q8993XQ55cMuVv3baU1bJ+4H7FH4LfibL3XvnG8Lu0fXXOe3Yo8xjaqaUcLJ+jrFt+mf0ISUgk8BxWroxL0IvSamn14vfrZX5vagjG3M8thp24Wtkn3K/h2vSMhF5449fL+XLQBXOq+SHFylS1Fr7ARE1Rvx5CRE93d9+mohuPdt+FEVRlE+X89XQM6211URE/f8OOZOhMebrxphtxphtva3tZzJTFEVRLpBP/aWotfZ31tpp1tpp7oS4s39BURRFOS/OyW3RGFNARCuYhn6IiK6y1lYbY4YS0fvW2tFn289nQUOP5FK5NHoToKM9NHGl2BZi4tyvnkMNw2CmFFiHrkX75ExodtNmybqP29fh0vSkQ0OeXFQh7JJ80EPXlckIUJ6JP24rIhMn/N1+YbZ+V5HT9jbhPELJsu9FRXAhq22D213zMVnUwZ3KIkVZNr/hM48Ju5IyVhOUTWlPvMyA6N3PXlBMwjuEzroIH04WvWo6pFh6+1wUg1hbjfcEJ2tk7UhiEb+BGmjDGbPke4d4H86xbB2iUnk0MRGRO4j93bZwo9N+6cAkYZfyDsbJvRTRhidLZJ1Py8+R1Ty1cfJaxZZChw6myz4VXM6u40pE25rZUjcO7YTfBD//6mLpHvu9Ja867Z+tXOy0fUPl6r63HHMm83Lm3rhO7i9+NiI268vxLskVkUXx6rnQ/9ccQW3PmGI5L3zz4CLZWC3nKq+Hyt007xkls6k+ugUFLvjvKjOrSdidrMKYVd7z0KfqtvgaES3rby8jolc/xlZRFEUZAM7FbfE5ItpIRKONMVXGmLuI6GdEtNAYc5iIFvZ/VhRFUS4hZ03OZa294wybFlzkvlwSuMwykC6N3t1YMv6k+Ra5kRcaZImgYo7Jy9UdD23BVwB3wZIGubQetgpuckdvwz56Ivyu1m5HkQxXp1yS9qZBuuhmK82NZTKx2JIZO5z2ioMT0Pf9EfVQ/5SOYz0MdzLrkhKgPYGlaw8rFJDkj3A7Y26gEy+DlLR7b4EwS74Cy+L6EziR7DXy2cbXDGmh9uvSLXBnI6SFiemo37myMlXYGRZx2ROL8/K6paRxoBTJxOJYcrOwR/YpexqkihXPz8ZxkuSY1U/BZ/cu+CvEjJJuqp4Pcf7heVju9+yU0lFwLMb6trE7xbYPHkNCrsQ7ER0ZDMm5On4hirEcbcE4ubvkPFt+P5J4BWbi/G+Zt0fY7clAsY6KlQVOe9rivcJu0xr46QbY2PZOkNc0zgPZKzUJ8k6bX0ou8b/F2HT+vSxIUcPkozCr97siYYKwS9zJ3CJZLrGT4QjJzivlrXNBI0UVRVGiBL2hK4qiRAl6Q1cURYkSBkW2xfPl03Rp5NndYoojigG345qcmgrt2gQjYo2Zbjx8HNzHyqtlxjlXNXToWVciZL74ZVlMwt+A47ZeJ93EujtYFkSWRXD6ZWXCblsZ3O5siD0vRGTpG7IOGmsjkxhnzj0g7Ha+Cl2/lw/TOFnI2LsJYmQwA+cR9kZo8pnQSr1lGJeubOneGJMMF870P0od9fUnn3DaU57/Lo4VKzXP9M24XjFLoS+f3CZd63i4fyLLHHhqb7qwC2ezYhpHMBg82yARUcCLuXVFRrnTPtwm4//2lLJ3AYVwA921Txb0iBkCrbizQc7VrLU4x5prcdzkNKlR+16EC97mn/3GaRetXSbsYjciViXEkkjOWCw19PXvY+6GUnBcT7PU7tPGw22x9UN2/lObhV14F94n8AyigVr5zNvD3luFEuT1Hvoh5ngNXi2Q9cg56OpmWn4SLv6d02Uxlj+vneO0K+57ULMtKoqiDCb0hq4oihIlqORyjlxs+cVfjaVh2l55DWpugSzgP4gl7pzFu4TdphcmOu3OLCz/fPlyuevZgBqRY2475LS37pXRoFdORKXKHa9KOaYjh9UsZZkTc0bI5X7DRsgJnklwhevdLl2yAjOQwfAUcx8kj1zG+hLgctndAImERzYSEQWyWTbIRoxZIEUWWjCsnmfvQcg03K2QiMiVhe/Fx8l9XJkjsxaeZtXr08XnYC5kHB+LFA2PlC6XOc9CzmoYi3bbCFmohMMjbQ8fyxTbYhLQ364uJpVVywIX3hY8z337i6877bJO6fb6yjpkTgyclM+AnWNwrFgW/dx9UNYlzZsB987bsoud9v9ec5OwSzwECSeuBnPhxDVyXsQcx3jOuhl1Tbc/K4uMdM+Hq2ZhBubqnkOy6ERMKq6Jewv67pJKHGUtqnTapXtkHdowq0PrZoUwZhSUC7vtVTi2bxt0pbbLZBZOHvFc+Q8Pq+SiKIoymNAbuqIoSpSgkst5cDEiSk0n/i+9b8E7YtsLx6Y47ePl8HSIqZJv8IduwBI37n9gSfulrM3C7l9eRpGIUCKWhSm58k2/m0Vp1tcliG3+OCwHu5ik4U6Qa9JwExI5xRzH8nnoBilbHFnGJJNO2CWWyHNsK8BSO/sD9C/t/nJhV9cJ74jGVrR5XVMiot5q9N2m4ZwCh6Uc0ZmLJfNXZ60X27bcCbecpl/g/GvKpXdRXCZkoKTlGE8uqRERecpx7DArU2kjSn4mlOMPLYUYl6wN8jfcPJKN51HYeYJStvA34PxLvwJpxhUrpZ6cF7At8bsyKVrla4gU7pwGbxjfHukZ1JWGY4cD6G/hc3Isym/B9ekZim2Lxkkvl7vT1zntO377gNP+8df+JOwe2YFSDSPvQYEQV7qM6j3wIKTCseMgq1ScksXYug5AHpxzjYxKbe7GdSzeh3H55YLlwu57KxF872tk0cQj5W/kL7N/67RnFFSq5KIoijKY0Bu6oihKlKA3dEVRlChBNfSLwPm4NBpW/MCVLN2VxuayAgCtcKHq7pGRojfkI6py53+Drlt2W7ywGzcL0ZwuAy1z74eFwo7/986j74iIYiugo05ehKIWfpe02/AW3MZCicyV8lSEu1s+tOfsPLgwev6vjI48MQ/nPIQVMggtl656vX+DfTQ1QUM3EY8scdug0bZOhmbp9cvzSEmAHlyQJEvq7n5rjNPmbm024lg8SyHfv98nj9VeKgslnCaxSBaJ6NoIjZ5rxQ++t1TY+epZYREW9ehJk+6SV+SXO+0N6xGRG/bJewJ3bwwlyG2/uPHPTvsHu1CMJRSSc9XFipMk4LDUtkBGJAf8GNDOTryPiVsnq5352tCP2vn4TtJOWah7zFK44m4uYZlBw/IFRdpGzO+OodiWd1WlsCuvh/YeqpVRsykFuF6JAej/nU8PFXYjvgHX4R3H4fqY9Lr83TahVgyV/uAfVUNXFEUZTOgNXVEUJUpQyeUic67yS0omItgyfiKXifWTsfQ6NR5L5iEjG6TdfsgTeRMh01QckMmf3OlY/vV0wi3Q3SRdBHuTWZKjOq/YlszKlDZMZbUoE6VcZJnbYkIZlt3L7npb2D25c77THjEUEXzlJ6Xrn48twXnUYzjCHTFhB1zGWkfANTMyotSbiyV+dwf6+l9+BiwR2rA35CbPd5Boq6oBEbDhCikL9AYwTvwa2IiIzYRCRNTmJ2PZfvhtGcnbkY/rw8+raGyVsDt2ihVhaIMfpEiWRkSGFVCw7HyNPyLJWAbman2pvD7eFlaQg0XA9rTJ+ZOZi/NqaML8HpEpI41LKjB34w/g+rQPk0VB0kdABqurYZJVxDm64jB//nY8IlS9Ru7vhRLUZeX1Sm2ulKm8JXDHtGNlgjj+PX7tA7UR8hOT6bomQtpLTJAFMwJ/hLyz6S+anEtRFGVQoTd0RVGUKEEll0+Rj4sozc9HnubqLfIteOo+XJPb/3ml0/4/H1wr7BKzseTrZjUcQ+XybTlfXfZkMokkYnk6/HksE4/eHvF/PVMuXO1YQoaTI7IXsaV7IB3L1YRYGQXXUIKlO/eimH9jsbDbuHwyjsucQ1qmyf3ZNiYfsRzovlopKy26AVG0L22fiv4NkQnNvjAcCZ+e2TpLbLtyHPSnvc+gZmXzbNmneOZRM+fLqLX67urJwi5jErx3wiw8tL1LSnHjMyD1bN48Gt9JlF4zE0ZBgtm7h+Wnj8jLfcdM5N9+fjVqlHrb5LUP5mHO5ORIj5+67fA24l4+t1wno5U3/BIJvuLvRlRzpMTm8WKyBnnudZ+UgWJLMTa5C+CJUtcuZa/Qh9h/B0tgN2LicWF35ABqlPqZRBLMl5JiOpNKG0tltKmvCQPQVSAjYDmzi5DcbdsqVsd3nKz5GmSSYMVXvq+Si6IoymBCb+iKoihRgt7QFUVRogTV0AcQrqlfs/XrTtv7nowUdIXYNVkEV8WOoF/YBVvxOSUdenrrQantseBQij0BjfbyO2S2OB4t6O6U7n4iCyArcOEKSbvUy/FuYFgiXNUCbqm1b1uBAhpdl0FrH51zUtgdb8bYXJmDiNcN1QXCrrkUWfGSCnHctr1yLIZMwf5PNiIKt7dRjq2X6aExJ+U5tkxlLojc3S8gXeGIFdOwXbDznZS6/sIboK+vLkd4YLBeRiIWFsE1tWotiiRMvzHiOm7AdUwqguZ9Tc5hYffiRujaCTnQb9uPyPn46OJnnPYDW/9ObOtlbqrZrNhJfbN8j/OjSSig8cj7f+u0U7Nlxs/OzXDF5Rkv07dI17/6uazWLsvWaSMyRaZsQf9uufcDp72mpkjYJQcwBxs74ZpYXRfx2zzOMmPmRrzHOYVjedoxfwK1cv50sKyUvR2YC7FH5DsTD/Ni3Pv4A6qhK4qiDCb0hq4oihIlDKjkYoypI6J2Iqo/m+0gIZ10LE6jYwF0LIAJGKeXAAAD7UlEQVSORR/51tqMsxkN6A2diMgYs+1ctKDBgI4F0LEAOhZAx+KToZKLoihKlKA3dEVRlCjhUtzQf3cJjvlZRccC6FgAHQugY/EJGHANXVEURfl0UMlFURQlStAbuqIoSpQwoDd0Y8wNxphDxphSY8zDA3nsS40xJs8Ys8YYc8AYs88Yc3//31ONMauMMYf7/005276iBWOM2xhTbIxZ0f95uDFmc/9Y/MUY4zvbPqIBY0yyMeYFY8zB/vkxa7DOC2PMd/t/H3uNMc8ZYwKDdV6cDwN2QzfGuInoSSK6kYjGEtEdxpixH/+tqKKHiP7RWnsZEV1BRN/qP/+HiWi1tXYUEa3u/zxYuJ+IDrDP/0ZEj/WPxSkiuuuS9GrgeYKI3rbWjiGiidQ3JoNuXhhjcojoPiKaZq0dT0RuIlpKg3defGIG8gl9BhGVWmvLrLXdRLSciJYM4PEvKdbaamvtjv52K/X9aHOobwye7jd7mohuvTQ9HFiMMblEdDMR/b7/syGia4johX6TQTEWxphEIppHRH8gIrLWdltrm2iQzgsi8hBRjDHGQ0SxRFRNg3BenC8DeUPPIaJj7HNV/98GHcaYAiKaTESbiSjTWltN1HfTJ6Ihl65nA8rjRPRPRHQ6F2QaETVZa0+nyxss82MEEdUR0f/rl59+b4yJo0E4L6y1x4nol0RUSX038mYi2k6Dc16cFwN5Qzcf8bdB5zNpjIknoheJ6DvW2paz2UcjxphFRFRrrd3O//wRpoNhfniIaAoR/cZaO5n6ch1FvbzyUfS/J1hCRMOJKJuI4qhPoo1kMMyL82Igb+hVRJTHPucS0YkBPP4lxxjjpb6b+bPW2pf6/3zSGDO0f/tQIqq9VP0bQOYQ0WJjTDn1SW/XUN8Te3L/Upto8MyPKiKqstaeLsL5AvXd4AfjvLiWiI5aa+ustSEieomIZtPgnBfnxUDe0LcS0aj+N9Y+6nvZ8doAHv+S0q8R/4GIDlhrH2WbXiOiZf3tZUT06kD3baCx1j5irc211hZQ3zx4z1p7JxGtIaLT1Q8Gy1jUENExY8zpqs8LiGg/DcJ5QX1SyxXGmNj+38vpsRh08+J8Gej0uTdR35OYm4iestb+dMAOfokxxswlonVEtIegG3+f+nT054loGPVN6NuttY0fuZMoxBhzFRE9aK1dZIwZQX1P7KlEVExEX7bWnrl8epRgjJlEfS+HfURURkRfo76HrUE3L4wxPyKiL1KfV1gxEd1NfZr5oJsX54OG/iuKokQJGimqKIoSJegNXVEUJUrQG7qiKEqUoDd0RVGUKEFv6IqiKFGC3tAVRVGiBL2hK4qiRAn/CQeaKGlTo9DMAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "pt.imshow(np.hstack([np.eye(k), proj]))"
      ]
    },
    {
      "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
}