III_08_Left_and_right_inverses_Pseudoinverses.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ This notebook is part of lecture 32 *Left-, right-, and pseudoinverses* in the OCW MIT course 18.06 by Prof Gilbert Strang [1]\n",
    "+ Created by me, Dr Juan H Klopper\n",
    "    + Head of Acute Care Surgery\n",
    "    + Groote Schuur Hospital\n",
    "    + University Cape Town\n",
    "    + <a href=\"mailto:juan.klopper@uct.ac.za\">Email me with your thoughts, comments, suggestions and corrections</a> \n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" /></a><br /><span xmlns:dct=\"http://purl.org/dc/terms/\" href=\"http://purl.org/dc/dcmitype/InteractiveResource\" property=\"dct:title\" rel=\"dct:type\">Linear Algebra OCW MIT18.06</span> <span xmlns:cc=\"http://creativecommons.org/ns#\" property=\"cc:attributionName\">IPython notebook [2] study notes by Dr Juan H Klopper</span> is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\">Creative Commons Attribution-NonCommercial 4.0 International License</a>.\n",
    "\n",
    "+ [1] <a href=\"http://ocw.mit.edu/courses/mathematics/18-06sc-linear-algebra-fall-2011/index.htm\">OCW MIT 18.06</a>\n",
    "+ [2] Fernando Pérez, Brian E. Granger, IPython: A System for Interactive Scientific Computing, Computing in Science and Engineering, vol. 9, no. 3, pp. 21-29, May/June 2007, doi:10.1109/MCSE.2007.53. URL: http://ipython.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #ddd;\n",
       "    color: #000000;\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 2.2em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0, 80, 120);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 1.9em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(200,100,0);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(94,127,192);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "        font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML, Image\n",
    "css_file = 'style.css'\n",
    "HTML(open(css_file, 'r').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import init_printing, Matrix, symbols, sqrt, Rational\n",
    "from numpy import matrix, transpose, sqrt\n",
    "from numpy.linalg import pinv, inv, det, svd, norm\n",
    "from scipy.linalg import pinv2\n",
    "from warnings import filterwarnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "init_printing(use_latex = 'mathjax')\n",
    "filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Left- and right-sided inverses and pseudoinverses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The inverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Recall the four fundamental subspaces\n",
    "    + The rowspace (with **x**) and nullspace in &#8477;<sup>n</sup>\n",
    "    + The columnspace (with A**x**) and the nullspace of A<sup>T</sup> in &#8477;<sup>m</sup>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ The two-sided inverse gives us the following\n",
    "$$ {A}{A}^{-1}=I={A}^{-1}{A} $$\n",
    "    + For this we need *r* = *m* = *n* (i.e. full rank)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ For a left-inverse we have the following\n",
    "    + Full column rank, with *r* = *n* (but possibly more rows)\n",
    "    + The nullspace contains just the zero vector (columns are independent)\n",
    "    + The rows might not all be independent\n",
    "    + We thus have either no or only a single solution to A**x**=**b**\n",
    "    + A<sup>T</sup> will now also have full rank\n",
    "    + From (A<sup>T</sup>A)<sup>-1</sup>A<sup>T</sup>A = I follows the fact that (A<sup>T</sup>A)<sup>-1</sup>A<sup>T</sup> is a left-sided inverse (A<sup>-1</sup>)\n",
    "    + Note, though, that (A<sup>T</sup>A)<sup>-1</sup>A<sup>T</sup> is a *n* &times; *m* matrix and A is of size *m* &times; *n*, resulting in a *n* &times; *n* identity matrix\n",
    "    + We cannot do AA<sup>-1</sup> and have a *n* &times; *n* identity matrix, though, but instead will be a projection matrix (onto the columnspace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ For a right-inverse we have the following\n",
    "    + Full row rank, with *r* = *m* < *n*\n",
    "    + The nullspace of A<sup>T</sup> is the zero vector (rows are independent)\n",
    "    + Elimination will result in many solutions to A**x**=**b** (*n* - *m* free variables)\n",
    "    + Now there will be an A<sup>-1</sup> to the right of A to give I\n",
    "    + AA<sup>T</sup>(AA<sup>T</sup>)<sup>-1</sup> = I\n",
    "    + A<sup>-1</sup> is now A<sup>T</sup>(AA<sup>T</sup>)<sup>-1</sup>\n",
    "    + Putting the right-inverse on the left is also a projection (onto the rowspace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The pseudoinverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Consider a matrix where *r* is less than *m* and *n*\n",
    "+ Remember that the rowspace is in &#8477;<sup>r</sup> and the columnspace is also in &#8477;<sup>r</sup>\n",
    "+ The nullspace of the rowspace is in &#8477;<sup>n-r</sup> and the nullspace of A<sup>T</sup> is in &#8477;<sup>m-r</sup>\n",
    "+ The rowspace and columnspace are in the same dimension and every vector **x** in one translate to another vector in the other (one-to-one)\n",
    "    + If **y** in another vector in the rowspace (not same as **x**) then A**x** &ne; A**y**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ The pseudoinverse A<sup>+</sup>, then, maps **x** (or **y**) from the columnspace to the rowspace\n",
    "$$ y={A}^{+}{Ay} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Suppose A**x** = A**y** or A(**x**-**y**) = 0\n",
    "    + Now (**x**-**y**) is in the nullspace *and* in the rowspace, i.e. it has to be the zero vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the pseudoinverse A<sup>+</sup>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ One way is to start from the singular value decomposition\n",
    "$$  {A}={U}{\\Sigma}{V}^{T} $$\n",
    "+ &Sigma; has along the main diagonal all the square roots of the eigenvalues and *r* pivots, but *m* row and *n* columns which can be more than *r*\n",
    "+ &Sigma;<sup>+</sup> will have 1 over the square roots of the eigenvalues along the main diagonals and then (possibly) zero values further along, but be of size *n* &times; *m*\n",
    "+ &Sigma;&Sigma;<sup>+</sup> will have 1<sup>'s</sup> along the main diagonal, and then 0<sup>'s</sup> (if larger tha *r*)\n",
    "    + It will be of size *m* &times; *m*\n",
    "    + It is a projection onto the columnspace\n",
    "+ &Sigma;<sup>+</sup>&Sigma; will also have 1<sup>'s</sup> along the main diagonal as well, but be of size *n* &times; *n*\n",
    "    + It is a projection onto the rowspace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ We now have the following\n",
    "$$ {A}^{+}={V}{\\Sigma}^{+}{U}^{T} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's see how easy this is in python™"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(matrix([[3, 6],\n",
       "         [2, 4]]), 0.0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = matrix([[3, 6], [2, 4]]) # Not sympy\n",
    "A, det(A) # The det is zero, so no inverse exists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.04615385,  0.03076923],\n",
       "        [ 0.09230769,  0.06153846]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The numpy pinv() function use SVD\n",
    "Aplus = pinv(A)\n",
    "Aplus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.04615385,  0.03076923],\n",
       "       [ 0.09230769,  0.06153846]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The scipy pinv2() function also uses SVD\n",
    "# The scipy pinv() function uses least squares to approxiamte\n",
    "# the pseudoinverse and as matrices get BIG, this\n",
    "# becomes computationally expensive\n",
    "Aplus_sp = pinv2(A)\n",
    "Aplus_sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Example problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example problem 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Calculate the pseudoinverse of A=[1,2]\n",
    "+ Calculate AA<sup>+</sup>\n",
    "+ Calculate A<sup>+</sup>A\n",
    "+ If **x** is in the nullspace of A what is the effect of A<sup>+</sup>A on **x** (i.e. A<sup>+</sup>A**x**)\n",
    "+ If **x** is in the columnspace of A<sup>T</sup> what is A<sup>+</sup>A**x**?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1, 2]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = matrix([1, 2])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's use singular value decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "U, S, VT = svd(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-1.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.23606798])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-0.4472136 , -0.89442719],\n",
       "        [-0.89442719,  0.4472136 ]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "VT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Remember,\n",
    "$$ {A}^{+}={V}{\\Sigma}^{+}{U}^{T} $$\n",
    "+ &Sigma; must be of size 2 &times; 1, though"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "S = matrix([[sqrt(5)], [0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.],\n",
       "        [ 2.]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aplus = transpose(VT) * S * U\n",
    "Aplus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ This needs to be normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$2.2360679775$$"
      ],
      "text/plain": [
       "2.2360679775"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(Aplus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.4472136 ],\n",
       "        [ 0.89442719]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 / norm(Aplus) * Aplus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.2],\n",
       "        [ 0.4]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aplus = pinv(A)\n",
    "Aplus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A * Aplus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.2,  0.4],\n",
       "        [ 0.4,  0.8]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aplus * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's create a vector in the nullspace of A\n",
    "    + It will be any vector\n",
    "    $$ c\\begin{bmatrix}-2\\\\1\\end{bmatrix} $$\n",
    "+ Let's choose the constant *c* = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.],\n",
       "        [ 0.]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_vect_null_A = matrix([[-2], [1]])\n",
    "Aplus * A * x_vect_null_A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ This is now surprise as A<sup>+</sup>A reflects a vector onto the rowspace of A\n",
    "    + We chose **x** in the nullspace of A, so A**x** must be **0** and A<sup>+</sup>A**x** = **0**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ The columnsapce of A<sup>T</sup> is any vector\n",
    "$$ c\\begin{bmatrix}1\\\\2\\end{bmatrix} $$\n",
    "+ We'll choose *c* = 1 again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.],\n",
       "        [ 2.]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_vect_null_AT = matrix([[1], [2]])\n",
    "Aplus * A * x_vect_null_AT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ We recover **x** again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ For fun, let's just check what A<sup>+</sup> is when A is invertible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = matrix([[1, 2], [3, 4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-2. ,  1. ],\n",
       "        [ 1.5, -0.5]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pinv(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-2. ,  1. ],\n",
       "        [ 1.5, -0.5]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(A)"
   ]
  },
  {
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}