Linear_Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Algebra via Examples\n",
    "\n",
    "We demonstrate a few important LA concepts via examples. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenvalues and Eigenvectors\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "M = np.mat(\"3 -2; 1 0\")\n",
    "w, V = np.linalg.eig(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.,  1.])"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.89442719,  0.70710678],\n",
       "        [ 0.4472136 ,  0.70710678]])"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify the property that $$\\mathbf{M} \\mathbf{v} = \\lambda \\mathbf{v} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.89442719],\n",
       "        [ 0.4472136 ]])"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V[:, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 2.68328157],\n",
       "        [ 1.34164079]])"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3 * V[:, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i = 1\n",
      "Left  =  [[ 1.78885438  0.89442719]]\n",
      "Right =  [[ 1.78885438  0.89442719]]\n",
      "i = 2\n",
      "Left  =  [[ 0.70710678  0.70710678]]\n",
      "Right =  [[ 0.70710678  0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    print('i = {}'.format(i+1))\n",
    "    print('Left  = ', (M @ V[:, i]).reshape(1,2))\n",
    "    print('Right = ', (w[i] * V[:, i]).reshape(1,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Impact to basis vectors\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 3, -2],\n",
       "        [ 1,  0]])"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.mat(\"1 0; 0 1\")\n",
    "V = M @ U\n",
    "V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1, -1, -1,  0,  0,  0,  1,  1,  1],\n",
       "       [-1,  0,  1, -1,  0,  1, -1,  0,  1]])"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U2 = np.hstack([np.array([i-1, j-1]).reshape(2,1) for i in range(3) for j in range(3)])\n",
    "U2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-1, -3, -5,  2,  0, -2,  5,  3,  1],\n",
       "        [-1, -1, -1,  0,  0,  0,  1,  1,  1]])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V2 = M @ U2\n",
    "V2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAG3FJREFUeJzt3X+UXOV93/H3xxISTW2DhBYsJBmJ\noGJI3AprROy6jb2YH4KcIuxgLHIoAkOpleAc13WCCFHtCnMC9LT4uCGATPkRx+anj8O6MUcBdmlO\ncwBrVMsgcIWEsKu1FLT8tBOM+PXtH/fZcJ9lZmdHd2ZXy35e58yZuc997r3fufNoPnvvzOgqIjAz\nMxv2rokuwMzM9i8OBjMzyzgYzMws42AwM7OMg8HMzDIOBjMzyzgY7B1L0vWS1na6b4v1LJQUkqZX\nXZfZRJF/x2DWOZIWAk8DB0TE6xNbjdm+8RGDvSNJmjbRNZhNVg4GmzQkHSPpQUkvSnpc0umlebdI\nuk7S9yX9A9Cb2r5a6vOHknZL2iXpwnTK56jS8l9Njz8uaVDSf5S0Jy1zfmk9vyXph5J+LmmnpK+0\n8RwukfQzSb+QtFXSJ1L7VyTdLemONO//SPoXpeXWSHoqzXtC0idHrPffSfpxaf6HUvvhkr4jaUjS\n05J+v939blOPg8EmBUkHAN8D/ho4FPg88C1JR5e6/Q5wBfAe4H+PWH458EXgROAo4GMtNvk+4CBg\nHnABcK2kWWnePwDnAgcDvwWslnTGGJ7D0cDFwLKIeA9wCvCTUpcVwF3AbODbwF+m5w3wFPCvU03/\nGfgLSXPTej8NfCXV9F7gdOA5Se+i2Gc/Ss/jE8AXJJ3Sqlab2hwMNll8GHg3cGVEvBoR/cD/BM4u\n9bknIv42It6MiFdGLH8WcHNEPB4RL1O8uY7mNWBdRLwWEd8H/h44GiAiHoyIx9J2HgVuo3XQALwB\nzASOlXRARPwkIp4qzd8UEXdHxGvAfwMOTM+biLgrInalbd4BbAOOT8tdCFwdERujsD0ifgosA3oi\nYl3aZzuAbwArx1CrTWEOBpssDgd2RsSbpbafUvwlPGxnq+XH2BfguREfHr9MEUxI+g1JA+n0zEvA\n54A5rZ5ARGwHvkDx1/0eSbdLOrxRTel5Dqa6kXSupM3pNNqLwK+XtrmA4ohipCOAw4eXScv9EXBY\nq1ptanMw2GSxC1iQTo8Mez/ws9L0aF+x2w3ML00vqFDLt4E+YEFEHARcD2gsC0bEtyPiX1G8aQdw\nVaOa0vOcD+ySdATFX/oXA4dExMHAltI2dwK/2mBzO4GnI+Lg0u09EXFaG8/VpiAHg00Wj1Cc2/9D\nSQdI+jjwb4Dbx7j8ncD56QPsXwH+U4Va3gM8HxGvSDqe4rONliQdLekESTOBV4BfUpxeGrZU0qfS\nbyC+AOwFHgb+KUWIDKX1nE9xxDDsRuBLkpaqcFQKkx8AP08feP8TSdMk/bqkZRWeu00BDgabFCLi\nVYoPVU8FngX+DDg3Iv7vGJe/F/g6MABsBx5Ks/buQzm/C6yT9AuKgLlzjMvNBK6kqP/vKD5E/6PS\n/HuAzwAvAP8W+FT6jOMJ4L+mmp8BPgj8bem53UXxofu3gV8AfwnMjog3KMJzCcVvK56lCJGD2n/K\nNpX4B242JUk6huJ0zMz94Ydo6SuvR0XEORNdi5mPGGzKkPRJSTPS106vAr63P4SC2f7GwWBTyb+n\nOE//FMW5/dUTW47Z/smnkszMLOMjBjMzy0zK/xp4zpw5sXDhwokuw8xsUtm0adOzEdHTqt+kDIaF\nCxdSr9cnugwzs0lF0k/H0s+nkszMLONgMDOzjIPBzMwyDgYzM8s4GMzMLNORYJB0U7oE4pYm8yXp\n65K2S3p0+LKDad4qSdvSbVUn6hnp6qthYCBvGxgo2s2q8viybpqI8dWpI4ZbgOWjzD8VWJxuFwHX\nAUiaDXwZ+A2Kq1F9uXT5xI5ZtgzOOuutnTswUEwv838+bB3g8WXdNCHjKyI6cgMWAluazLsBOLs0\nvRWYS3FZxhua9Wt2W7p0abSrvz9izpyItWuL+/7+tldh1pTHl3VTp8YXUI8xvJ+P12cM88gvpTiY\n2pq1v42kiyTVJdWHhobaLqC3F1avhssvL+57e9tehVlTHl/WTeM9vsYrGBpd9jBGaX97Y8T6iKhF\nRK2np+Uvut9mYACuuw7Wri3uR56zM6vC48u6abzH13gFwyD5NXbnU1zDt1l7Rw2fk7vzTli3rrgv\nn7Mzq8Ljy7ppIsbXeAVDH3Bu+nbSh4GXImI3sAE4WdKs9KHzyamtozZuLHbm8OFXb28xvXFjp7dk\nU5HHl3XTRIyvjlyPQdJtwMeBORTXpP0ycABARFwvScCfUnxz6WXg/Iiop2U/y1vXvb0iIm5utb1a\nrRb+T/TMzNojaVNE1Fr168j/rhoRZ7eYH8DvNZl3E3BTJ+owM7Pq/MtnMzPLOBjMzCzjYDAzs4yD\nwczMMg4GMzPLOBjMzCzjYDAzs4yDwczMMg4GMzPLOBjMzCzjYDAzs4yDwczMMg4GMzPLOBjMzCzj\nYDAzs4yDwczMMh0JBknLJW2VtF3Smgbzr5G0Od2elPRiad4bpXl9najHzMz2XeUruEmaBlwLnAQM\nAhsl9UXEE8N9IuI/lPp/HjiutIpfRsSSqnWYmVlndOKI4Xhge0TsiIhXgduBFaP0Pxu4rQPbNTOz\nLuhEMMwDdpamB1Pb20g6AlgE9JeaD5RUl/SwpDOabUTSRalffWhoqANlm5lZI50IBjVoiyZ9VwJ3\nR8Qbpbb3R0QN+B3ga5J+tdGCEbE+ImoRUevp6alWsZmZNdWJYBgEFpSm5wO7mvRdyYjTSBGxK93v\nAB4k//zBzMzGWSeCYSOwWNIiSTMo3vzf9u0iSUcDs4CHSm2zJM1Mj+cAHwWeGLmsmZmNn8rfSoqI\n1yVdDGwApgE3RcTjktYB9YgYDomzgdsjonya6RjgBklvUoTUleVvM5mZ2fhT/j49OdRqtajX6xNd\nhpnZpCJpU/pMd1T+5bOZmWUcDGZmlnEwmJlZxsFgZmYZB4OZmWUcDGZmlnEwmJlZxsFgZmYZB4OZ\nmWUcDGZmlnEwmJlZxsFgZmYZB4OZmWUcDGZmlnEwmJlZxsFgZmaZjgSDpOWStkraLmlNg/nnSRqS\ntDndLizNWyVpW7qt6kQ9Zma27ypf2lPSNOBa4CRgENgoqa/BJTrviIiLRyw7G/gyUAMC2JSWfaFq\nXWZmtm86ccRwPLA9InZExKvA7cCKMS57CnBfRDyfwuA+YHkHajIzs33UiWCYB+wsTQ+mtpF+W9Kj\nku6WtKDNZZF0kaS6pPrQ0FAHyjYzs0Y6EQxq0BYjpr8HLIyIfw7cD9zaxrJFY8T6iKhFRK2np2ef\nizUzs9F1IhgGgQWl6fnArnKHiHguIvamyW8AS8e6rJmZja9OBMNGYLGkRZJmACuBvnIHSXNLk6cD\nP06PNwAnS5olaRZwcmozM7MJUvlbSRHxuqSLKd7QpwE3RcTjktYB9YjoA35f0unA68DzwHlp2ecl\nXU4RLgDrIuL5qjWZmdm+U0TDU/r7tVqtFvV6faLLMDObVCRtiohaq37+5bOZmWUcDGZmlnEwmJlZ\nxsFgZmYZB4OZmWUcDGZmlnEwmJlZxsFgZmYZB4OZmWUcDGZmlnEwmJlZxsFgZmYZB4OZmWUcDGZm\nlnEwmJlZpiPBIGm5pK2Stkta02D+FyU9IelRSQ9IOqI07w1Jm9Otb+SyZmY2vipfwU3SNOBa4CSK\nazhvlNQXEU+Uuv0QqEXEy5JWA1cDn0nzfhkRS6rWYWZmndGJI4bjge0RsSMiXgVuB1aUO0TEQES8\nnCYfBuZ3YLtmZtYFnQiGecDO0vRgamvmAuDe0vSBkuqSHpZ0RrOFJF2U+tWHhoaqVWxmZk1VPpUE\nqEFbwwtJSzoHqAEfKzW/PyJ2SToS6Jf0WEQ89bYVRqwH1kNxzefqZZuZWSOdOGIYBBaUpucDu0Z2\nknQicBlwekTsHW6PiF3pfgfwIHBcB2oyM7N91Ilg2AgslrRI0gxgJZB9u0jSccANFKGwp9Q+S9LM\n9HgO8FGg/KG1mZmNs8qnkiLidUkXAxuAacBNEfG4pHVAPSL6gP8CvBu4SxLA/4uI04FjgBskvUkR\nUleO+DaTmZmNM0VMvtP1tVot6vX6RJdhZjapSNoUEbVW/fzLZzMzyzgYzMws42AwM7OMg8HMzDIO\nBjMzyzgYzMws42AwM7OMg8HMzDIOBjMzyzgYzMws42AwM7OMg8HMzDIOBjMzyzgYzMws42AwM7OM\ng8HMzDIdCQZJyyVtlbRd0poG82dKuiPNf0TSwtK8S1P7VkmndKKeka6+GgYG8raBgaLdrCqPL+um\niRhflYNB0jTgWuBU4FjgbEnHjuh2AfBCRBwFXANclZY9luIa0b8GLAf+LK2vo5Ytg7POemvnDgwU\n08uWdXpLNhV5fFk3Tcj4iohKN+AjwIbS9KXApSP6bAA+kh5PB54FNLJvud9ot6VLl0a7+vsj5syJ\nWLu2uO/vb3sVZk15fFk3dWp8AfUYw/t6J04lzQN2lqYHU1vDPhHxOvAScMgYlwVA0kWS6pLqQ0ND\nbRfZ2wurV8Pllxf3vb1tr8KsKY8v66bxHl+dCAY1aIsx9hnLskVjxPqIqEVEraenp80Si8Ov666D\ntWuL+5Hn7Myq8Piybhrv8dWJYBgEFpSm5wO7mvWRNB04CHh+jMtWNnxO7s47Yd264r58zs6sCo8v\n66aJGF+dCIaNwGJJiyTNoPgwuW9Enz5gVXp8JtCfznf1ASvTt5YWAYuBH3SgprzAjcXOHD786u0t\npjdu7PSWbCry+LJumojxpeL9ueJKpNOArwHTgJsi4gpJ6yg+6OiTdCDwTeA4iiOFlRGxIy17GfBZ\n4HXgCxFxb6vt1Wq1qNfrles2M5tKJG2KiFrLfp0IhvHmYDAza99Yg8G/fDYzs4yDwczMMg4GMzPL\nOBjMzCzjYDAzs4yDwczMMg4GMzPLOBjMzCzjYDAzs4yDwczMMg4GMzPLOBjMzCzjYDAzs4yDwczM\nMg4GMzPLOBjMzCxTKRgkzZZ0n6Rt6X5Wgz5LJD0k6XFJj0r6TGneLZKelrQ53ZZUqcfMzKqresSw\nBnggIhYDD6TpkV4Gzo2IXwOWA1+TdHBp/h9ExJJ021yxHjMzq6hqMKwAbk2PbwXOGNkhIp6MiG3p\n8S5gD9BTcbtmZtYlVYPhsIjYDZDuDx2ts6TjgRnAU6XmK9IppmskzRxl2Ysk1SXVh4aGKpZtZmbN\ntAwGSfdL2tLgtqKdDUmaC3wTOD8i3kzNlwIfAJYBs4FLmi0fEesjohYRtZ4eH3CYmXXL9FYdIuLE\nZvMkPSNpbkTsTm/8e5r0ey/wV8AfR8TDpXXvTg/3SroZ+FJb1ZuZWcdVPZXUB6xKj1cB94zsIGkG\n8F3gzyPirhHz5qZ7UXw+saViPWZmVlHVYLgSOEnSNuCkNI2kmqQbU5+zgN8EzmvwtdRvSXoMeAyY\nA3y1Yj1mZlaRImKia2hbrVaLer0+0WWYmU0qkjZFRK1VP//y2czMMg4GMzPLOBjMzCzjYDAzs4yD\nwczMMg4GMzPLOBjMzCzjYDAzs4yDwczMMg4GMzPLOBjMzCzjYDAzs4yDwczMMg4GMzPLOBjMzCxT\nKRgkzZZ0n6Rt6X5Wk35vlC7S01dqXyTpkbT8Helqb2ZmNoGqHjGsAR6IiMXAA2m6kV9GxJJ0O73U\nfhVwTVr+BeCCivWYmVlFVYNhBXBrenwrxXWbxyRd5/kE4O59Wd7MzLqjajAcFhG7AdL9oU36HSip\nLulhScNv/ocAL0bE62l6EJjXbEOSLkrrqA8NDVUs28zMmpneqoOk+4H3NZh1WRvbeX9E7JJ0JNAv\n6THg5w36Nb0AdUSsB9ZDcc3nNrZtZmZtaBkMEXFis3mSnpE0NyJ2S5oL7Gmyjl3pfoekB4HjgO8A\nB0uano4a5gO79uE5mJlZB1U9ldQHrEqPVwH3jOwgaZakmenxHOCjwBMREcAAcOZoy5uZ2fiqGgxX\nAidJ2gaclKaRVJN0Y+pzDFCX9COKILgyIp5I8y4BvihpO8VnDv+jYj1mZlaRij/cJ5darRb1en2i\nyzAzm1QkbYqIWqt+/uWzmZllHAxmZpZxMJiZWcbBYGZmGQeDmZllHAxmZpZxMJiZWcbBYGZmGQeD\nmZllHAxmZpZxMJiZWcbBYGZmGQeDmZllHAxmZpZxMJiZWcbBYGZmmUrBIGm2pPskbUv3sxr06ZW0\nuXR7RdIZad4tkp4uzVtSpR4zM6uu6hHDGuCBiFgMPJCmMxExEBFLImIJcALwMvDXpS5/MDw/IjZX\nrMfMzCqqGgwrgFvT41uBM1r0PxO4NyJerrhdMzPrkqrBcFhE7AZI94e26L8SuG1E2xWSHpV0jaSZ\nzRaUdJGkuqT60NBQtarNzKyplsEg6X5JWxrcVrSzIUlzgQ8CG0rNlwIfAJYBs4FLmi0fEesjohYR\ntZ6ennY2bWZmbZjeqkNEnNhsnqRnJM2NiN3pjX/PKKs6C/huRLxWWvfu9HCvpJuBL42xbjMz65Kq\np5L6gFXp8SrgnlH6ns2I00gpTJAkis8ntlSsx8zMKqoaDFcCJ0naBpyUppFUk3TjcCdJC4EFwP8a\nsfy3JD0GPAbMAb5asR4zM6uo5amk0UTEc8AnGrTXgQtL0z8B5jXod0KV7ZuZWef5l89mZpZxMJiZ\nWcbBYGZmGQeDmZllHAxmZpZxMJiZWcbBYGZmGQeDmZllHAxmZpZxMJiZWcbBYGZmGQeDmZllHAxm\nZpZxMJiZWcbBYGZmGQeDmZllKgWDpE9LelzSm5Jqo/RbLmmrpO2S1pTaF0l6RNI2SXdImlGlnmau\nvhoGBvK2gYGi3awqjy/rpokYX1WPGLYAnwL+plkHSdOAa4FTgWOBsyUdm2ZfBVwTEYuBF4ALKtbT\n0LJlcNZZb+3cgYFietmybmzNphqPL+umCRlfEVH5BjwI1JrM+wiwoTR9aboJeBaY3qjfaLelS5dG\nu/r7I+bMiVi7trjv7297FWZNeXxZN3VqfAH1GMN77Hh8xjAP2FmaHkxthwAvRsTrI9obknSRpLqk\n+tDQUNtF9PbC6tVw+eXFfW9v26swa8rjy7ppvMdXy2CQdL+kLQ1uK8a4DTVoi1HaG4qI9RFRi4ha\nT0/PGDf9loEBuO46WLu2uB95zs6sCo8v66bxHl/TW3WIiBMrbmMQWFCang/sojiNdLCk6emoYbi9\n44bPyd15Z5G0vb35tFkVHl/WTRMxvsbjVNJGYHH6BtIMYCXQl853DQBnpn6rgHu6UsDGfCf29hbT\nGzd2Y2s21Xh8WTdNxPhS8f68jwtLnwT+O9ADvAhsjohTJB0O3BgRp6V+pwFfA6YBN0XEFan9SOB2\nYDbwQ+CciNjbaru1Wi3q9fo+121mNhVJ2hQRTX9a8I/9qgTDRHEwmJm1b6zB4F8+m5lZxsFgZmYZ\nB4OZmWUcDGZmlpmUHz5LGgJ+uo+Lz6H4DcX+xnW1x3W1x3W1551a1xER0fIXwpMyGKqQVB/Lp/Lj\nzXW1x3W1x3W1Z6rX5VNJZmaWcTCYmVlmKgbD+okuoAnX1R7X1R7X1Z4pXdeU+4zBzMxGNxWPGMzM\nbBQOBjMzy7wjg0HSpyU9LulNSU2/2iVpuaStkrZLWlNqXyTpEUnbJN2R/rvwTtQ1W9J9ab33SZrV\noE+vpM2l2yuSzkjzbpH0dGnekvGqK/V7o7TtvlL7RO6vJZIeSq/3o5I+U5rX0f3VbLyU5s9Mz397\n2h8LS/MuTe1bJZ1SpY59qOuLkp5I++cBSUeU5jV8TceprvMkDZW2f2Fp3qr0um+TtGqc67qmVNOT\nkl4szevK/pJ0k6Q9krY0mS9JX081PyrpQ6V5nd9XY7n+52S7AccARzP6tainAU8BRwIzgB8Bx6Z5\ndwIr0+PrgdUdqutqYE16vAa4qkX/2cDzwK+k6VuAM7uwv8ZUF/D3TdonbH8B/wxYnB4fDuwGDu70\n/hptvJT6/C5wfXq8ErgjPT429Z8JLErrmTaOdfWWxtDq4bpGe03Hqa7zgD9tsOxsYEe6n5Uezxqv\nukb0/zzFpQK6vb9+E/gQsKXJ/NOAeymufPlh4JFu7qt35BFDRPw4Ira26HY8sD0idkTEqxTXhVgh\nScAJwN2p363AGR0qbUVa31jXeyZwb0S83KHtN9NuXf9oovdXRDwZEdvS413AHorrg3Raw/EySr13\nA59I+2cFcHtE7I2Ip4HtaX3jUldEDJTG0MMUV0vstrHsr2ZOAe6LiOcj4gXgPmD5BNV1NnBbh7bd\nVET8DcUfgc2sAP48Cg9TXP1yLl3aV+/IYBijecDO0vRgajsEeDGKy42W2zvhsIjYDZDuD23RfyVv\nH5RXpEPJayTNHOe6DpRUl/Tw8Okt9qP9Jel4ir8Cnyo1d2p/NRsvDfuk/fESxf4Zy7LdrKvsAoq/\nPIc1ek3Hs67fTq/P3ZKGLwG8X+yvdMptEdBfau7W/mqlWd1d2Vctr/m8v5J0P/C+BrMui4ixXCJU\nDdpilPbKdY11HWk9c4EPAhtKzZcCf0fx5rceuARYN451vT8idqm48l6/pMeAnzfoN1H765vAqoh4\nMzXv8/5qtIkGbSOfZ1fGVAtjXrekc4Aa8LFS89te04h4qtHyXajre8BtEbFX0ucojrZOGOOy3axr\n2Erg7oh4o9TWrf3VyriOrUkbDBFxYsVVDAILStPzgV0U/0HVwZKmp7/6htsr1yXpGUlzI2J3eiPb\nM8qqzgK+GxGvlda9Oz3cK+lm4EvjWVc6VUNE7JD0IHAc8B0meH9Jei/wV8Afp8Ps4XXv8/5qoNl4\nadRnUNJ04CCK0wNjWbabdSHpRIqw/ViULp/b5DXtxBtdy7oi4rnS5DeAq0rLfnzEsg92oKYx1VWy\nEvi9ckMX91crzeruyr6ayqeSNgKLVXyjZgbFIOiL4hOdAYrz+wCrgLEcgYxFX1rfWNb7tnOb6c1x\n+Lz+GUDDbzB0oy5Js4ZPxUiaA3wUeGKi91d67b5Lcf71rhHzOrm/Go6XUeo9E+hP+6cPWKniW0uL\ngMXADyrU0lZdko4DbgBOj4g9pfaGr+k41jW3NHk68OP0eANwcqpvFnAy+ZFzV+tKtR1N8WHuQ6W2\nbu6vVvqAc9O3kz4MvJT+8OnOvurGJ+wTfQM+SZGke4FngA2p/XDg+6V+pwFPUiT+ZaX2Iyn+4W4H\n7gJmdqiuQ4AHgG3pfnZqrwE3lvotBH4GvGvE8v3AYxRvcH8BvHu86gL+Zdr2j9L9BfvD/gLOAV4D\nNpduS7qxvxqNF4pTU6enxwem57897Y8jS8telpbbCpza4fHeqq7707+D4f3T1+o1Hae6/gR4PG1/\nAPhAadnPpv24HTh/POtK018BrhyxXNf2F8UfgbvTWB6k+Czoc8Dn0nwB16aaH6P0bctu7Cv/lxhm\nZpaZyqeSzMysAQeDmZllHAxmZpZxMJiZWcbBYGZmGQeDmZllHAxmZpb5/2McAyjZ+RRrAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x151a55d6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(U2[0,:], U2[1,:], 'bx')\n",
    "plt.title('original space')\n",
    "plt.draw()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFphJREFUeJzt3X+0ZWV93/H3x0FAo5HBGRWZ0cGE\nVk20UE/QlsYfyAg1liGr/gCXOmostY1Gm6BiaBZTKg1JmmhtXBVCRkdwiYpxMUZdBFBM0oiZM4r8\nssA4ahgH5SqithgQ+PaPs2et81zvnTt3zrn3cOe+X2uddfZ+nmef/d3D4nzOfvY5d6eqkCRpj4dN\nugBJ0kOLwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqSGwaADTpJvJjkzyfVJfpjk\no0kOHep/SZLrktyd5O+SPLNrf12STw2N25HkY0Prtyc5Zob9HZrkkiTf715zW5LHd33XJPn9JH/f\n1XJ5ksOHtv14ku90fX+d5JeG+h6R5I+TfKvr/9skj+j6ntPVfneSryZ5/pj/GbWMGQw6UL0cOBk4\nCngm8FqAJP8c2Az8e+CxwAXA1iSHAF8AfjXJw5IcATwcOL7b7inAo4DrZ9jXRuAxwNruNd8I/GSo\n/zXA64EnAvcD7x3q+yxwNPA44MvAh4f6/jvwLOBfAocDbwceTHIk8GngXV37mcAnkqyex7+PNCuD\nQQeq91bV7qq6C/gUsOeT/r8DLqiqL1XVA1W1BbgXeE5V7QR+3I19HnAF8O0kT+3W/6aqHpxhXz9l\nEAi/2L3m9qr60VD/xVV1Y1X9P+D3gJcnWQFQVZur6sdVdS+wCfhnSR6T5GEMwuQtVfXt7nX/rhv3\nKuAzVfWZqnqwqq4E+sCLx/WPp+XNYNCB6jtDy/cw+LQP8GTgd7opmLuT3M3gk/4Tu/4vAM8Hntst\nX8MgFJ7Xrc/kYgYhcmmS3Un+MMnDh/pvH1r+FoMzkVVJViQ5P8nXk/wI+GY3ZlX3OBT4+gz7ezLw\nsmnH8K+AI2b915DmwWDQcnM7cF5VHTb0eGRVfaTr3xMMv9otf4E5gqGqflpV/6Wqns5g2uclDKaP\n9lg7tPwkBmcY3wNeCWwATmQwFbWuG5Ou/x+BX5jlGC6edgw/V1Xnz+PfQZqVwaDl5s+ANyZ5dgZ+\nLsmvJXl01/8F4AXAI6pqF/A3DK5VPBb4ykwvmOQFSZ7RTQ/9iMEb/wNDQ16V5OlJHgmcC1xWVQ8A\nj2YwjfV94JHAf9uzQTdltRn4kyRP7M4u/kV3LeQS4N8kOalrPzTJ85OsGdO/kZY5g0HLSlX1GVxn\n+FPgB8AOugvTXf+twP9lEAh01wp2Av+7ezOfyROAyxiEwtcYhMslQ/0XAx9kML11KPBbXfuHGEwt\nfRu4Gbh22uueCdwAbAPuAv4AeFhV3c7gTON3gSkGZxBvw/+fNSbxRj3SwklyDXBJVV006VqkfeUn\nDElSw2CQJDWcSpIkNTxjkCQ1Dpp0Aftj1apVtW7dukmXIUlLyvbt279XVXP+6ZQlGQzr1q2j3+9P\nugxJWlKSfGtfxjmVJElqGAySpIbBIElqGAySpIbBIElqjCUYkmxOcmeSG2fpT5L3drdKvL67i9ae\nvo1JbuseG8dRj7Qsbdo06QoW33I75kU63rH88jnJcxn8RcoPVdUvz9D/YuDNDO4w9Wzgf1TVs7t7\n3/aBHlDAduBZVfWDve2v1+uVX1eVpklguf0lg+V2zCMeb5LtVdWba9xYzhiq6q8Z/Fng2WxgEBpV\nVdcCh3X31D0JuLKq7urC4EoGf/tekjQhi3WN4Uja2xvu6tpma/8ZSc5I0k/Sn5qaWrBCpSVl06bB\np8hksL5n+UCeYlluxzyB4x3bH9FLsg74y1mmkj4N/H5V/W23fjXwduAE4JCqelfX/nvAPVX1x3vb\nl1NJ0gyW27QKLL9jXkpTSftgF+19b9cAu/fSLkmakMUKhq3Aa7pvJz0H+GFV3QFcAbwoycokK4EX\ndW2S5uuccyZdweJbbse8SMc7rm8lfQR4PrAK+C5wDvBwgKp6f5IwuMfuycA9wOu6e++S5PUM7l0L\ncF5VfWCu/TmVJEnzt69TSWP566pVdfoc/QX85ix9m4HN46hDkjQ6f/ksSWoYDJKkhsEgSWoYDJKk\nhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkxliC\nIcnJSW5JsiPJWTP0vzvJdd3j1iR3D/U9MNS3dRz1SJL238h3cEuyAngfsB7YBWxLsrWqbt4zpqr+\n09D4NwPHDr3ET6rqmFHrkCSNxzjOGI4DdlTVzqq6D7gU2LCX8acDHxnDfiVJC2AcwXAkcPvQ+q6u\n7WckeTJwFPC5oeZDk/STXJvk1Nl2kuSMblx/ampqDGVLkmYyjmDIDG01y9jTgMuq6oGhtidVVQ94\nJfCeJL8w04ZVdWFV9aqqt3r16tEqliTNahzBsAtYO7S+Btg9y9jTmDaNVFW7u+edwDW01x8kSYts\nHMGwDTg6yVFJDmbw5v8z3y5K8k+BlcAXh9pWJjmkW14FHA/cPH1bSdLiGflbSVV1f5I3AVcAK4DN\nVXVTknOBflXtCYnTgUurania6WnABUkeZBBS5w9/m0mStPjSvk8vDb1er/r9/qTLkKQlJcn27pru\nXvnLZ0lSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUM\nBklSw2CQJDUMBklSw2CQJDXGEgxJTk5yS5IdSc6aof+1SaaSXNc93jDUtzHJbd1j4zjqkSTtv5Fv\n7ZlkBfA+YD2wC9iWZOsMt+j8aFW9adq2hwPnAD2ggO3dtj8YtS5J0v4ZxxnDccCOqtpZVfcBlwIb\n9nHbk4Arq+quLgyuBE4eQ02SpP00jmA4Erh9aH1X1zbdv01yfZLLkqyd57YkOSNJP0l/ampqDGVL\nkmYyjmDIDG01bf1TwLqqeiZwFbBlHtsOGqsurKpeVfVWr16938VKkvZuHMGwC1g7tL4G2D08oKq+\nX1X3dqt/BjxrX7eVJC2ucQTDNuDoJEclORg4Ddg6PCDJEUOrpwBf65avAF6UZGWSlcCLujZJ0oSM\n/K2kqro/yZsYvKGvADZX1U1JzgX6VbUV+K0kpwD3A3cBr+22vSvJf2UQLgDnVtVdo9YkSdp/qZpx\nSv8hrdfrVb/fn3QZkrSkJNleVb25xvnLZ0lSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQ\nJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSYyzBkOTkJLck2ZHkrBn6fzvJzUmuT3J1\nkicP9T2Q5LrusXX6tpKkxTXyHdySrADeB6xncA/nbUm2VtXNQ8O+AvSq6p4k/wH4Q+AVXd9PquqY\nUeuQJI3HOM4YjgN2VNXOqroPuBTYMDygqj5fVfd0q9cCa8awX0nSAhhHMBwJ3D60vqtrm81vAJ8d\nWj80ST/JtUlOnW2jJGd04/pTU1OjVSxJmtXIU0lAZmib8UbSSV4F9IDnDTU/qap2J3kK8LkkN1TV\n13/mBasuBC6EwT2fRy9bkjSTcZwx7ALWDq2vAXZPH5TkROBs4JSqundPe1Xt7p53AtcAx46hJknS\nfhpHMGwDjk5yVJKDgdOA5ttFSY4FLmAQCncOta9Mcki3vAo4Hhi+aC1JWmQjTyVV1f1J3gRcAawA\nNlfVTUnOBfpVtRX4I+BRwMeTAPxDVZ0CPA24IMmDDELq/GnfZpIkLbJULb3p+l6vV/1+f9JlSNKS\nkmR7VfXmGucvnyVJDYNBktQwGCRJDYNBktQwGCRJDYNBktQwGCRJDYNBktQwGCRJDYNBktQwGCRJ\nDYNBktQwGCRJDYNBktQwGCRJDYNBktQYSzAkOTnJLUl2JDlrhv5Dkny06/9SknVDfe/s2m9JctI4\n6tFDxKZNk65g8XnMOgCMfAe3JCuAW4H1wC4G94A+ffgWnUn+I/DMqnpjktOAX6+qVyR5OvAR4Djg\nicBVwD+pqgf2tk/v4LZEJLAE7xA4Eo9ZD2GLeQe344AdVbWzqu4DLgU2TBuzAdjSLV8GvDCDmz9v\nAC6tqnur6hvAju71JEkTMo5gOBK4fWh9V9c245iquh/4IfDYfdwWgCRnJOkn6U9NTY2hbC2ITZsG\nnyCTwfqe5QN5usFjXh7HvIyMYyrpZcBJVfWGbv3VwHFV9eahMTd1Y3Z1619ncGZwLvDFqrqka/9z\n4DNV9Ym97dOppCViOU4xeMx6CFvMqaRdwNqh9TXA7tnGJDkIeAxw1z5uK0laROMIhm3A0UmOSnIw\ncBqwddqYrcDGbvmlwOdqcKqyFTit+9bSUcDRwN+PoSY9FJxzzqQrWHwesw4AI08lASR5MfAeYAWw\nuarOS3Iu0K+qrUkOBS4GjmVwpnBaVe3stj0beD1wP/DWqvrsXPtzKkmS5m9fp5LGEgyLzWCQpPlb\nzGsMkqQDiMEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoY\nDJKkhsEgSWoYDJKkhsEgSWqMFAxJDk9yZZLbuueVM4w5JskXk9yU5Pokrxjq+2CSbyS5rnscM0o9\nkqTRjXrGcBZwdVUdDVzdrU93D/Caqvol4GTgPUkOG+p/W1Ud0z2uG7EeSdKIRg2GDcCWbnkLcOr0\nAVV1a1Xd1i3vBu4EVo+4X0nSAhk1GB5fVXcAdM+P29vgJMcBBwNfH2o+r5tieneSQ/ay7RlJ+kn6\nU1NTI5YtSZrNnMGQ5KokN87w2DCfHSU5ArgYeF1VPdg1vxN4KvArwOHAO2bbvqourKpeVfVWr/aE\nQ5IWykFzDaiqE2frS/LdJEdU1R3dG/+ds4z7eeDTwH+uqmuHXvuObvHeJB8AzpxX9ZKksRt1Kmkr\nsLFb3ghcPn1AkoOBTwIfqqqPT+s7onsOg+sTN45YjyRpRKMGw/nA+iS3Aeu7dZL0klzUjXk58Fzg\ntTN8LfXDSW4AbgBWAe8asR5J0ohSVZOuYd56vV71+/1JlyFJS0qS7VXVm2ucv3yWJDUMBklSw2CQ\nJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDVG\nCoYkhye5Mslt3fPKWcY9MHSTnq1D7Ucl+VK3/Ue7u71JkiZo1DOGs4Crq+po4OpufSY/qapjuscp\nQ+1/ALy72/4HwG+MWI8kaUSjBsMGYEu3vIXBfZv3SXef5xOAy/Zne0nSwhg1GB5fVXcAdM+Pm2Xc\noUn6Sa5NsufN/7HA3VV1f7e+Czhyth0lOaN7jf7U1NSIZUuSZnPQXAOSXAU8YYaus+exnydV1e4k\nTwE+l+QG4EczjJv1BtRVdSFwIQzu+TyPfUuS5mHOYKiqE2frS/LdJEdU1R1JjgDunOU1dnfPO5Nc\nAxwLfAI4LMlB3VnDGmD3fhyDJGmMRp1K2gps7JY3ApdPH5BkZZJDuuVVwPHAzVVVwOeBl+5te0nS\n4ho1GM4H1ie5DVjfrZOkl+SibszTgH6SrzIIgvOr6uau7x3AbyfZweCaw5+PWI8kaUQZfHBfWnq9\nXvX7/UmXIUlLSpLtVdWba5y/fJYkNQwGSVLDYJAkNQwGSVLDYJAkNQwGSVLDYJAkNQwGSVLDYJAk\nNQwGSVLDYJAkNQwGSVLDYJAkNQwGSVLDYJAkNQwGSVJjpGBIcniSK5Pc1j2vnGHMC5JcN/T4xySn\ndn0fTPKNob5jRqlHkjS6Uc8YzgKurqqjgau79UZVfb6qjqmqY4ATgHuAvxoa8rY9/VV13Yj1SJJG\nNGowbAC2dMtbgFPnGP9S4LNVdc+I+5UkLZBRg+HxVXUHQPf8uDnGnwZ8ZFrbeUmuT/LuJIfMtmGS\nM5L0k/SnpqZGq1qSNKs5gyHJVUlunOGxYT47SnIE8AzgiqHmdwJPBX4FOBx4x2zbV9WFVdWrqt7q\n1avns2tJ0jwcNNeAqjpxtr4k301yRFXd0b3x37mXl3o58Mmq+unQa9/RLd6b5APAmftYtyRpgYw6\nlbQV2NgtbwQu38vY05k2jdSFCUnC4PrEjSPWI0ka0ajBcD6wPsltwPpunSS9JBftGZRkHbAW+MK0\n7T+c5AbgBmAV8K4R65EkjWjOqaS9qarvAy+cob0PvGFo/ZvAkTOMO2GU/UuSxs9fPkuSGgaDJKlh\nMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiS\nGgaDJKkxUjAkeVmSm5I8mKS3l3EnJ7klyY4kZw21H5XkS0luS/LRJAePUs8+2bRpwXfxkOMxS5qH\nVNX+b5w8DXgQuAA4s7tz2/QxK4BbGdz6cxewDTi9qm5O8jHgL6rq0iTvB75aVf9rrv32er3q939m\nV/taNIxwzEuSxywJSLK9qmb9EL/HSGcMVfW1qrpljmHHATuqamdV3QdcCmxIEuAE4LJu3Bbg1FHq\nkSSNbjGuMRwJ3D60vqtreyxwd1XdP619RknOSNJP0p+amppfBZs2DT5BJntebPA4kKcbPOblcczS\nAphzKinJVcATZug6u6ou78Zcw+xTSS8DTqqqN3Trr2ZwFnEu8MWq+sWufS3wmap6xlxFO5U0Tx6z\nJPZ9KumguQZU1Ykj1rILWDu0vgbYDXwPOCzJQd1Zw552SdIELcZU0jbg6O4bSAcDpwFba3Cq8nng\npd24jcDlC17NOecs+C4ecjxmSfMw6reSfh34n8Bq4G7guqo6KckTgYuq6sXduBcD7wFWAJur6ryu\n/SkMLkYfDnwFeFVV3TvXfkeaSpKkZWpfp5JGCoZJMRgkaf4W5euqkqQDj8EgSWoYDJKkhsEgSWos\nyYvPSaaAb026jnlaxeC3G8uJx7w8eMxLx5OravVcg5ZkMCxFSfr78m2AA4nHvDx4zAcep5IkSQ2D\nQZLUMBgWz4WTLmACPOblwWM+wHiNQZLU8IxBktQwGCRJDYNhApKcmaSSrJp0LQstyR8l+T9Jrk/y\nySSHTbqmhZLk5CS3JNmR5KxJ17PQkqxN8vkkX0tyU5K3TLqmxZBkRZKvJPnLSdeyUAyGRdbdqW49\n8A+TrmWRXAn8clU9E7gVeOeE61kQSVYA7wP+NfB04PQkT59sVQvufuB3quppwHOA31wGxwzwFuBr\nky5iIRkMi+/dwNuBZXHVv6r+aui+3tcyuFPfgeg4YEdV7ayq+xjcZ2TDhGtaUFV1R1V9uVv+MYM3\ny1nv234gSLIG+DXgoknXspAMhkWU5BTg21X11UnXMiGvBz476SIWyJHA7UPruzjA3ySHJVkHHAt8\nabKVLLj3MPhg9+CkC1lIc97zWfOT5CrgCTN0nQ38LvCixa1o4e3tmKvq8m7M2QymHj68mLUtoszQ\ntizOCpM8CvgE8Naq+tGk61koSV4C3FlV25M8f9L1LCSDYcyq6sSZ2pM8AzgK+GoSGEypfDnJcVX1\nnUUscexmO+Y9kmwEXgK8sA7cH87sAtYOra8Bdk+olkWT5OEMQuHDVfUXk65ngR0PnNLdqvhQ4OeT\nXFJVr5pwXWPnD9wmJMk3gV5VLcW/0LjPkpwM/AnwvKqamnQ9CyXJQQwurr8Q+DawDXhlVd000cIW\nUAafcLYAd1XVWyddz2LqzhjOrKqXTLqWheA1Bi20PwUeDVyZ5Lok7590QQuhu8D+JuAKBhdhP3Yg\nh0LneODVwAndf9vruk/TWuI8Y5AkNTxjkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1/j/Z\n60gNTpfBxAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x151ad11860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(V2[0,:], V2[1,:], 'r+')\n",
    "plt.title('new space')\n",
    "plt.draw()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGL1JREFUeJzt3X20ZXV93/H3RxBoa6ODMypPAhpq\nxJpCOLC0diVe5EmrQFMdIRrHp9BFNY2JGkGlM4I0SNpiU0mUGJSmRCAa6xjrQmAuSf8IdM4o8mSQ\nkaiMg3IRH2IxKPLtH2ePnn3nnrn3zjn3HO6d92utvc7ev/3b5/fd9661P3fvfc7dqSokSdrhcZMu\nQJL02GIwSJJaDAZJUovBIElqMRgkSS0GgySpxWDQipPkhUm2TbqOUVhJ+6Llw2DQHi3JR5O8d54+\nleTnx1TP2MaSBjEYpD1Ekr0nXYOWB4NBy1KSryY5N8mdSb6T5CNJ9hvQ99lJbkzy3SR3JDm1aT8L\neBXwu0l+kOTTc2z7183sF5s+r2zafyPJ1iQPJtmY5MC+bU5KcleS7yX5wyR/leSNfetfn+RLTd3X\nJjl0V2M1696a5P4k9yV5XV/7vkn+c5KvJ/lWkg8m+UfNuhcm2ZbkHUm+CXxkt37Y2uMYDFrOXgWc\nDDwT+GfAu2d3SPJ44NPA54CnAL8JXJnkWVV1GXAlcHFVPaGqXjZ7+6r65Wb2XzR9rk5yPPB7wFrg\nAOBrwFXNeKuBjwPnAk8G7gL+ZV89pwPvBH4VWAP8H+Bjg8Zqlp8GPBE4CHgDcGmSVc269zX7fhTw\n802f/9i3C08D9gcOBc4a+JOU+hgMWs4+UFX3VtWDwIXAmXP0eR7wBOCiqvpRVW0C/nJA34V6FXB5\nVX2+qh6mFwLPT3IY8BLgjqr6i6p6BPgD4Jt92/474Peq6kvN+v8EHLXjrGGAHwPnV9WPq+p/Az8A\nnpUkwG8Av11VD1bV3zfvd0bfto8C66vq4ar64RD7rD2I1xy1nN3bN/814MA5+hwI3FtVj87qe9AQ\n4x4IfH7HQlX9IMm3m/c8sL+uqqpZnyo6FPhvSf5LX1uabb82YLxvNyGyw0P0wm4N8I+BLb2M+Ol7\n7dXXd6aq/mER+yYZDFrWDumbfzqwfY4+24FDkjyuLxyeDny5md+dfy+8nd4BHoAk/4TeZaNvAPcB\nB/etS/8yvdC4sKqu3I1xZ3sA+CHwnKr6xoA+/vtkLZqXkrScvSnJwUn2p3fd/uo5+twM/D96N5gf\nn+SFwMto7gkA3wKeMc84s/v8GfC6JEcl2Zfe5Zubq+qrwGeA5yY5vfkU0JvoXeff4YPAuUmeA5Dk\niUlesYuxBmqC7o+BS5I8pXm/g5KcvJDtpUEMBi1nf0bvpvI9zbTT9xGq6kfAqcCL6f2F/YfAa6rq\nb5sufwIc2Xxi6X8NGGcDcEXTZ21V3QCcB3yC3hnCM2mu61fVA8ArgIuBbwNHAl3g4Wb9J+ndML4q\nyfeB25va5hxrAT+DdwBbgZua97seeNYCtpMGig/q0XKU5KvAG6vq+knXsitJHgdsA15VVdOTrkda\nCM8YpBFLcnKSJzWXmd5J74bwTRMuS1owg0EavecDX6F36eplwOl+VFTLiZeSJEktnjFIklqW5fcY\nVq9eXYcddtiky5CkZWXLli0PVNWa+foty2A47LDD6Ha7ky5DkpaVJIO+Xd/ipSRJUovBIElqMRgk\nSS0GgySpxWCQJLWMJBiSXN48dvD2AeuT5A+aRyHemuSX+tatS3J3M60bRT3aM118MUzP+m9E09O9\n9pU47k5jb9gw1rEfE5p93mOMa3+raugJ+GXgl4DbB6x/CfBZev8z5nn0/kUx9B45eE/zuqqZXzXf\neMccc0xJs23aVLV6de91ruWVNu5OY8FYx35M6D0Lac8x5P4C3VrIMX0hnRb0RnDYLoLhQ8CZfct3\n0XtW7pnAhwb1GzQZDBpkx4HxvPPGe4Cc1Lj9YxfsWaFQZTAsevOFBcO47jEcRPsxjNuatkHtO0ly\nVpJuku7MzMySFarlbWoKzj4bLrig9zo1tbLHZcMGpo4PMw/0Hu0580CYOj4r+xLLhg2Q9Cb42fxK\n3edJ7O9C0mMhE7s+Y/gM8K/6lm8AjgHeDry7r/084K3zjeUZgwbxjMEzhhVthZ0xbKP9fN6D6T03\nd1C7tGjT07B2LVxzDZx/fu917dqdbwyvlHFnjw3jHVsr17iCYSPwmubTSc8DvldV9wHXAiclWZVk\nFXBS0yYt2ubNvQPjjss4U1O95c2bV+a4O429fv1Yx35MWL9+0hWM15j2dyTPY0jyMeCFwGp6DzNf\nDzweoKo+mCTAB4BTgIeA11VVt9n29fSecgVwYVV9ZL7xOp1O+U/0JGlxkmypqs58/Uby31Wr6sx5\n1hfwpgHrLgcuH0UdkqTh+c1nSVKLwSBJajEYJEktBoMkqcVgkCS1GAySpBaDQZLUYjBIkloMBklS\ni8EgSWoxGCRJLQaDJKnFYJAktRgMkqQWg0GS1GIwSJJaRhIMSU5JcleSrUnOmWP9JUluaaYvJ/lu\n37qf9K3bOIp6JEm7b+gnuCXZC7gUOBHYBmxOsrGq7tzRp6p+u6//bwJH973FD6vqqGHrkCSNxijO\nGI4DtlbVPVX1I+Aq4LRd9D8T+NgIxpUkLYFRBMNBwL19y9uatp0kORQ4HNjU17xfkm6Sm5KcPmiQ\nJGc1/bozMzMjKFuSNJdRBEPmaKsBfc8APl5VP+lre3pVdYBfA96f5JlzbVhVl1VVp6o6a9asGa5i\nSdJAowiGbcAhfcsHA9sH9D2DWZeRqmp783oPcCPt+w+SpDEbRTBsBo5IcniSfegd/Hf6dFGSZwGr\ngL/pa1uVZN9mfjXwAuDO2dtKksZn6E8lVdUjSd4MXAvsBVxeVXckOR/oVtWOkDgTuKqq+i8zPRv4\nUJJH6YXURf2fZpIkjV/ax+nlodPpVLfbnXQZkrSsJNnS3NPdJb/5LElqMRgkSS0GgySpxWCQJLUY\nDJKkFoNBktRiMEiSWgwGSVKLwSBJajEYJEktBoMkqcVgkCS1GAySpBaDQZLUYjBIkloMBklSy0iC\nIckpSe5KsjXJOXOsf22SmSS3NNMb+9atS3J3M60bRT2SpN039KM9k+wFXAqcCGwDNifZOMcjOq+u\nqjfP2nZ/YD3QAQrY0mz7nWHrkiTtnlGcMRwHbK2qe6rqR8BVwGkL3PZk4LqqerAJg+uAU0ZQkyRp\nN40iGA4C7u1b3ta0zfZvk9ya5ONJDlnktiQ5K0k3SXdmZmYEZUuS5jKKYMgcbTVr+dPAYVX1i8D1\nwBWL2LbXWHVZVXWqqrNmzZrdLlaStGujCIZtwCF9ywcD2/s7VNW3q+rhZvGPgWMWuq0kabxGEQyb\ngSOSHJ5kH+AMYGN/hyQH9C2eCnypmb8WOCnJqiSrgJOaNknShAz9qaSqeiTJm+kd0PcCLq+qO5Kc\nD3SraiPwH5KcCjwCPAi8ttn2wSQX0AsXgPOr6sFha5Ik7b5UzXlJ/zGt0+lUt9uddBmStKwk2VJV\nnfn6+c1nSVKLwSBJajEYJEktBoMkqcVgkCS1GAySpBaDQZLUYjBIkloMBklSi8EgSWoxGCRJLQaD\nJKnFYJAktRgMkqQWg0GS1DKSYEhySpK7kmxNcs4c638nyZ1Jbk1yQ5JD+9b9JMktzbRx9raSpPEa\n+gluSfYCLgVOpPcM581JNlbVnX3dvgB0quqhJGcDFwOvbNb9sKqOGrYOSdJojOKM4Thga1XdU1U/\nAq4CTuvvUFXTVfVQs3gTcPAIxpUkLYFRBMNBwL19y9uatkHeAHy2b3m/JN0kNyU5fdBGSc5q+nVn\nZmaGq1iSNNDQl5KAzNE254Okk7wa6AC/0tf89KranuQZwKYkt1XVV3Z6w6rLgMug98zn4cuWJM1l\nFGcM24BD+pYPBrbP7pTkBOBdwKlV9fCO9qra3rzeA9wIHD2CmiRJu2kUwbAZOCLJ4Un2Ac4AWp8u\nSnI08CF6oXB/X/uqJPs286uBFwD9N60lSWM29KWkqnokyZuBa4G9gMur6o4k5wPdqtoI/D7wBODP\nkwB8vapOBZ4NfCjJo/RC6qJZn2aSJI1Zqpbf5fpOp1PdbnfSZUjSspJkS1V15uvnN58lSS0GgySp\nxWCQJLUYDJKkFoNBktRiMEiSWgwGSVKLwSBJajEYJEktBoMkqcVgkCS1GAySpBaDQZLUYjBIkloM\nBklSi8EgSWoZSTAkOSXJXUm2JjlnjvX7Jrm6WX9zksP61p3btN+V5ORR1KPHiA0bxjrcxRfD9HS7\nbXq61z62cZt9Hse4O43dGNfYPzXm37PGoKqGmug9zvMrwDOAfYAvAkfO6vPvgQ8282cAVzfzRzb9\n9wUOb95nr/nGPOaYY0rLAIx1uE2bqlav7r3OtTyWcWFs4+409hzLYzHm37N2H73HLc9/XF9Ip12+\nATwfuLZv+Vzg3Fl9rgWe38zvDTwAZHbf/n67mgyGZWICB4wdB8bzzhvvAXLHuAVjPzBPap9/ymBY\nNhYaDKO4lHQQcG/f8rambc4+VfUI8D3gyQvcFoAkZyXpJunOzMyMoGwtiQ0bIOlN8LP5MV1umJqC\ns8+GCy7ovU5NjWHQDRuYOj7MPNDb55kHwtTxK3+fJ/l71tIaRTBkjrZaYJ+FbNtrrLqsqjpV1Vmz\nZs0iS9TYbNgAvb8he8s75sd0wJiehj/6IzjvvN7r7OvvS2LDBqY3FWtW9/Z5zepietPK3+dJ/p61\ntEYRDNuAQ/qWDwa2D+qTZG/gicCDC9xWWpDpaVi7Fq65Bs4/v/e6du3SHyj7x4XxjTt77HHus1a2\nUQTDZuCIJIcn2YfezeWNs/psBNY18y8HNjXXuzYCZzSfWjocOAL4vyOoSY8F69ePdbjNm3sHxh2X\nUqamesubN49x3PXrxzbuTmMzvn1uGfPvWUsvVXNeuVncmyQvAd5P7xNKl1fVhUnOp3ejY2OS/YA/\nBY6md6ZwRlXd02z7LuD1wCPAW6rqs/ON1+l0qtvtDl23JO1Jkmypqs68/UYRDONmMEjS4i00GPzm\nsySpxWCQJLUYDJKkFoNBktRiMEiSWgwGSVKLwSBJajEYJEktBoMkqcVgkCS1GAySpBaDQZLUYjBI\nkloMBklSi8EgSWoxGCRJLUMFQ5L9k1yX5O7mddUcfY5K8jdJ7khya5JX9q37aJK/S3JLMx01TD2S\npOENe8ZwDnBDVR0B3NAsz/YQ8Jqqeg5wCvD+JE/qW//2qjqqmW4Zsh5J0pCGDYbTgCua+SuA02d3\nqKovV9Xdzfx24H5gzZDjSpKWyLDB8NSqug+geX3KrjonOQ7YB/hKX/OFzSWmS5Lsu4ttz0rSTdKd\nmZkZsmxJ0iDzBkOS65PcPsd02mIGSnIA8KfA66rq0ab5XOAXgGOB/YF3DNq+qi6rqk5Vddas8YRD\nkpbK3vN1qKoTBq1L8q0kB1TVfc2B//4B/X4O+Azw7qq6qe+972tmH07yEeBti6pekjRyw15K2gis\na+bXAZ+a3SHJPsAngf9RVX8+a90BzWvo3Z+4fch6JElDGjYYLgJOTHI3cGKzTJJOkg83fdYCvwy8\ndo6PpV6Z5DbgNmA18N4h65EkDSlVNekaFq3T6VS32510GZK0rCTZUlWd+fr5zWdJUovBIElqMRgk\nSS0GgySpxWCQJLUYDJKkFoNBktRiMEiSWgwGSVKLwSBJajEYJEktBoMkqcVgkCS1GAySpBaDQZLU\nMlQwJNk/yXVJ7m5eVw3o95O+h/Rs7Gs/PMnNzfZXN097kyRN0LBnDOcAN1TVEcANzfJcflhVRzXT\nqX3t7wMuabb/DvCGIeuRJA1p2GA4Dbiimb+C3nObF6R5zvPxwMd3Z3tJ0tIYNhieWlX3ATSvTxnQ\nb78k3SQ3Jdlx8H8y8N2qeqRZ3gYcNGigJGc179GdmZkZsmxJ0iB7z9chyfXA0+ZY9a5FjPP0qtqe\n5BnApiS3Ad+fo9/AB1BX1WXAZdB75vMixpYkLcK8wVBVJwxal+RbSQ6oqvuSHADcP+A9tjev9yS5\nETga+ATwpCR7N2cNBwPbd2MfJEkjNOylpI3AumZ+HfCp2R2SrEqybzO/GngBcGdVFTANvHxX20uS\nxmvYYLgIODHJ3cCJzTJJOkk+3PR5NtBN8kV6QXBRVd3ZrHsH8DtJttK75/AnQ9YjSRpSen+4Ly+d\nTqe63e6ky5CkZSXJlqrqzNfPbz5LkloMBklSi8EgSWoxGCRJLQaDJKnFYJAktRgMkqQWg0GS1GIw\nSJJaDAZJUovBIElqMRgkSS0GgySpxWCQJLUYDJKkFoNBktQyVDAk2T/JdUnubl5XzdFnKsktfdM/\nJDm9WffRJH/Xt+6oYeqRJA1v2DOGc4AbquoI4IZmuaWqpqvqqKo6CjgeeAj4XF+Xt+9YX1W3DFmP\nJGlIwwbDacAVzfwVwOnz9H858NmqemjIcSVJS2TYYHhqVd0H0Lw+ZZ7+ZwAfm9V2YZJbk1ySZN9B\nGyY5K0k3SXdmZma4qiVJA80bDEmuT3L7HNNpixkoyQHAc4Fr+5rPBX4BOBbYH3jHoO2r6rKq6lRV\nZ82aNYsZWpK0CHvP16GqThi0Lsm3khxQVfc1B/77d/FWa4FPVtWP+977vmb24SQfAd62wLolSUtk\n2EtJG4F1zfw64FO76Hsmsy4jNWFCktC7P3H7kPVIkoY0bDBcBJyY5G7gxGaZJJ0kH97RKclhwCHA\nX83a/soktwG3AauB9w5ZjyRpSPNeStqVqvo28KI52rvAG/uWvwocNEe/44cZX5I0en7zWZLUYjBI\nkloMBklSi8EgSWoxGCRJLQaDJKnFYJAktRgMkqQWg0GS1GIwSJJaDAZJUovBIElqMRgkSS0GgySp\nxWCQJLUYDJKklqGCIckrktyR5NEknV30OyXJXUm2Jjmnr/3wJDcnuTvJ1Un2GaaeBdmwYcmHeMzZ\nQ/b54otherpZaPZ5errXPrZxG+MYd9Jja+Ua9ozhduBXgb8e1CHJXsClwIuBI4EzkxzZrH4fcElV\nHQF8B3jDkPXM7z3vWfIhHnP2kH0+9lhYu7Y5UL7nPUxP95aPPXaM48LYxp302FrBqmroCbgR6AxY\n93zg2r7lc5spwAPA3nP129V0zDHH1G6D3d92udqD9nnTpqrVq6sKavXq3vI4xz3vvBrruJMeW8sL\n0K0FHGPHcY/hIODevuVtTduTge9W1SOz2ueU5Kwk3STdmZmZxVWwYQMkvan3Zr1pJV9i2UP3eer4\nMPNAb59nHghTx49nn6em4Oyz4YILeq9TU0s+5GNibK1Q8yUHcD29S0azp9P6+tzI4DOGVwAf7lv+\ndeC/A2uArX3thwC3LSTNPGNYpD1onz1j8IxBg7HAM4a9FxAcJwyZPduag/4OBwPb6V1GelKSvat3\n1rCjXdotO66vX3MNcHzvdcfyUv4V3T/u1FRvGse4kx5bK9c4LiVtBo5oPoG0D3AGsLFJr2ng5U2/\ndcCnlrya9euXfIjHnD1knzdv7jsgrl/P1FRvefPmMY4LYxt30mNr5Urv+LybGyf/hp9dFvoucEtV\nnZzkQHqXj17S9HsJ8H5gL+DyqrqwaX8GcBWwP/AF4NVV9fB843Y6nep2u7tdtyTtiZJsqaqBXy34\nab9hgmFSDAZJWryFBoPffJYktRgMkqQWg0GS1GIwSJJaluXN5yQzwNcmXccirab33Y09ifu8Z3Cf\nl49Dq2rNfJ2WZTAsR0m6C/k0wEriPu8Z3OeVx0tJkqQWg0GS1GIwjM9lky5gAtznPYP7vMJ4j0GS\n1OIZgySpxWCQJLUYDBOQ5G1JKsnqSdey1JL8fpK/TXJrkk8medKka1oqSU5JcleSrUnOmXQ9Sy3J\nIUmmk3wpyR1JfmvSNY1Dkr2SfCHJX066lqViMIxZkkOAE4GvT7qWMbkO+OdV9YvAl+k973vFSbIX\ncCnwYuBI4MwkR062qiX3CPDWqno28DzgTXvAPgP8FvClSRexlAyG8bsE+F1gj7jrX1Wfq5891/sm\nek/qW4mOo/eo2nuq6kf0njNy2oRrWlJVdV9Vfb6Z/3t6B8uBz21fCZIcDPxr4MOTrmUpGQxjlORU\n4BtV9cVJ1zIhrwc+O+kilshBwL19y9tY4QfJfkkOA44Gbp5sJUvu/fT+sHt00oUspXmf+azFSXI9\n8LQ5Vr0LeCdw0ngrWnq72ueq+lTT5130Lj1cOc7axihztO0RZ4VJngB8AnhLVX1/0vUslSQvBe6v\nqi1JXjjpepaSwTBiVXXCXO1JngscDnwxCfQuqXw+yXFV9c0xljhyg/Z5hyTrgJcCL6qV+8WZbcAh\nfcsHA9snVMvYJHk8vVC4sqr+YtL1LLEXAKc2jyreD/i5JP+zql494bpGzi+4TUiSrwKdqlqO/6Fx\nwZKcAvxX4FeqambS9SyVJHvTu7n+IuAbwGbg16rqjokWtoTS+wvnCuDBqnrLpOsZp+aM4W1V9dJJ\n17IUvMegpfYB4J8C1yW5JckHJ13QUmhusL8ZuJbeTdhrVnIoNF4A/DpwfPO7vaX5a1rLnGcMkqQW\nzxgkSS0GgySpxWCQJLUYDJKkFoNBktRiMEiSWgwGSVLL/wd3YbAYJYAIjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x151ad342e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(U2[0,:], U2[1,:], 'bx')\n",
    "plt.plot(V2[0,:], V2[1,:], 'r+')\n",
    "plt.title('plot together')\n",
    "plt.draw()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With some more efforts, one can draw the original coordinate systems and the transformed one in a more visual way. This will hopefully strengthen the idea that a matrix multiplication represents a **linear transformation**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5,  3],\n",
       "       [ 0, -1],\n",
       "       [12,  1],\n",
       "       [ 4,  2]])"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([ [5, 3], [0, -1], [12, 1], [4, 2] ])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We only do a reduced SVD: $$\\mathbf{A} = \\mathbf{U} \\mathbf{S} \\mathbf{V}^{\\top}, $$ \n",
    "so that $\\mathbf{U}$ is of the **same** shape as $\\mathbf{A}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, s, Vt = linalg.svd(A, full_matrices=False)\n",
    "S = np.diag(s)\n",
    "V = Vt.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.39607852,  0.69414082],\n",
       "       [ 0.01400816, -0.34519358],\n",
       "       [-0.86373785, -0.47426224],\n",
       "       [-0.31125955,  0.41723523]])"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 13.8536595 ,   2.84185124])"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 13.8536595 ,   0.        ],\n",
       "       [  0.        ,   2.84185124]])"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.98098881, -0.1940643 ],\n",
       "       [-0.1940643 ,  0.98098881]])"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Vt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.98098881, -0.1940643 ],\n",
       "       [-0.1940643 ,  0.98098881]])"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We demonstrate a few useful properties of SVD: \n",
    "\n",
    "* (1) Exact Recovery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 2) (2, 2) (2, 2)\n"
     ]
    }
   ],
   "source": [
    "print(U.shape, S.shape, Vt.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  5.,   3.],\n",
       "       [  0.,  -1.],\n",
       "       [ 12.,   1.],\n",
       "       [  4.,   2.]])"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reconA = np.round(U @ S @ Vt, 5)\n",
    "reconA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* (2) Approximation\n",
    "\n",
    "Another view of SVD is to decompose the original matrix into a sum of rank-1 matrices of decreasing \"importance\", i.e., \n",
    "$$\n",
    "    \\mathbf{A} = \\sum_{i} {\\sigma_i} \\left( \\mathbf{U}_i \\otimes \\mathbf{V}_i \\right), \n",
    "$$\n",
    "where $\\otimes$ is the vector outer product, i.e., \n",
    "$$\n",
    "\\mathbf{u} \\otimes \\mathbf{v} = \\mathbf{u} \\otimes \\mathbf{v}^{\\top}\n",
    "$$\n",
    "\n",
    "In the following, we will approximate $\\mathbf{A}$ by only one summand. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.39607852,  0.01400816, -0.86373785, -0.31125955])"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U1 = U[:, 0]\n",
    "U1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13.853659500444321"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s1 = s[0]\n",
    "s1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.98098881, -0.1940643 ])"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Vt1 = Vt[0, :]\n",
    "Vt1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  5.38282,   1.06486],\n",
       "       [ -0.19037,  -0.03766],\n",
       "       [ 11.73844,   2.32216],\n",
       "       [  4.23011,   0.83682]])"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reconA1 = np.round(s1 * np.outer(U1, Vt1), 5)\n",
    "reconA1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5,  3],\n",
       "       [ 0, -1],\n",
       "       [12,  1],\n",
       "       [ 4,  2]])"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* (3) Feature reduction\n",
    "\n",
    "We manually creat a matrix $\\mathbf{A}$ such that all but its last column is heavily correlated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.45990338  0.96384231  0.48995298  0.05508264  0.06889611  0.49147112\n",
      "  0.04126475  0.85105499]\n",
      "[ 0.30916353  0.86528242  0.46032828  0.48086873  0.83216986  0.59160002\n",
      "  0.28097732  0.89340609]\n",
      "\n",
      "[[ 0.44777444  0.46136587  0.45736377  0.43501112  0.30916353]\n",
      " [ 0.95876347  0.96037656  0.96619143  0.95710081  0.86528242]\n",
      " [ 0.49118062  0.48100437  0.51165972  0.49153326  0.46032828]\n",
      " [ 0.05763428  0.04672335  0.04842411  0.06212327  0.48086873]\n",
      " [ 0.06610218  0.07657449  0.0596183   0.05856341  0.83216986]\n",
      " [ 0.47279545  0.48355236  0.50616279  0.48448376  0.59160002]\n",
      " [ 0.03870542  0.0394704   0.04529884  0.04100914  0.28097732]\n",
      " [ 0.84777085  0.84486164  0.8449768   0.86221843  0.89340609]]\n"
     ]
    }
   ],
   "source": [
    "a = np.random.rand(8)\n",
    "b = np.random.rand(8)\n",
    "sigma = 0.01\n",
    "At = np.array( [a + sigma*np.random.randn(8) for _ in range(4)])\n",
    "At = np.vstack((At, b))\n",
    "A = At.T\n",
    "print(a)\n",
    "print(b)\n",
    "print()\n",
    "print(A)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approximate_one(A):\n",
    "    U, s, Vt = linalg.svd(A, full_matrices=False)\n",
    "    print('singular values are: \\n', s) \n",
    "    S = np.diag(s)\n",
    "    V = Vt.T\n",
    "    # only 1 summand\n",
    "    U1 = U[:, 0]\n",
    "    s1 = s[0]\n",
    "    Vt1 = Vt[0, :]\n",
    "    A1 = s1 * np.outer(U1, Vt1)\n",
    "    return A1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "singular values are: \n",
      " [ 3.44554418  0.82844891  0.03062468  0.02714454  0.01036905]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.41124606,  0.41233158,  0.41714834,  0.41314695,  0.44603497],\n",
       "       [ 0.91990736,  0.92233552,  0.93311004,  0.92415941,  0.99772591],\n",
       "       [ 0.47617328,  0.47743017,  0.48300741,  0.47837428,  0.51645464],\n",
       "       [ 0.14129983,  0.14167281,  0.1433278 ,  0.14195296,  0.15325294],\n",
       "       [ 0.22314108,  0.22373007,  0.22634364,  0.22417249,  0.24201745],\n",
       "       [ 0.49799175,  0.49930623,  0.50513902,  0.5002936 ,  0.54011882],\n",
       "       [ 0.09007498,  0.09031274,  0.09136776,  0.09049133,  0.09769478],\n",
       "       [ 0.840482  ,  0.84270051,  0.85254476,  0.84436693,  0.91158166]])"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A1 = approximate_one(A)\n",
    "A1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.03652838 -0.0490343  -0.04021542 -0.02186417  0.13687144]\n",
      " [-0.03885611 -0.03804104 -0.03308139 -0.03294139  0.13244349]\n",
      " [-0.01500734 -0.00357419 -0.02865231 -0.01315898  0.05612636]\n",
      " [ 0.08366556  0.09494946  0.09490369  0.07982969 -0.32761579]\n",
      " [ 0.15703889  0.14715559  0.16672534  0.16560909 -0.5901524 ]\n",
      " [ 0.0251963   0.01575387 -0.00102377  0.01580984 -0.05148121]\n",
      " [ 0.05136956  0.05084235  0.04606891  0.0494822  -0.18328254]\n",
      " [-0.00728885 -0.00216112  0.00756796 -0.01785149  0.01817556]]\n"
     ]
    }
   ],
   "source": [
    "print(A1 - A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.21958962603481635,\n",
       " 0.21589060023747919,\n",
       " 0.21229301357777805,\n",
       " 0.2249147927967553,\n",
       " -0.80891509158094188]"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[ np.sum( (A1-A)[:, i] ) for i in range(A.shape[1]) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see and expect, the error is pretty small except in the last column. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that scipy has two version of SVD: `scipy.linalg.svd` and `scipy.sparse.linalg.svds`. See http://fa.bianp.net/blog/2012/singular-value-decomposition-in-scipy/ for a detailed discussion of the differences. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}