{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Pseudoinverse and Least Squares"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "\n",
        "np.set_printoptions(precision=4, linewidth=100)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A = np.random.randn(5, 3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compute the SVD of $A$. Note that `numpy.linalg.svd` returns $V^T$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "U, singval, VT = la.svd(A)\n",
        "V = VT.T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's first understand the shapes of these arrays:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "(5, 5)\n",
            "(3,)\n",
            "(3, 3)\n"
          ]
        }
      ],
      "source": [
        "print(U.shape)\n",
        "print(singval.shape)\n",
        "print(V.shape) "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check the orthogonality of $U$ and $V$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  1.0000e+00,   0.0000e+00,  -6.9389e-17,   5.5511e-17,   5.5511e-17],\n",
              "       [  0.0000e+00,   1.0000e+00,  -1.9429e-16,  -6.9389e-17,  -5.5511e-17],\n",
              "       [ -6.9389e-17,  -1.9429e-16,   1.0000e+00,  -1.6653e-16,   5.5511e-17],\n",
              "       [  5.5511e-17,  -6.9389e-17,  -1.6653e-16,   1.0000e+00,   0.0000e+00],\n",
              "       [  5.5511e-17,  -5.5511e-17,   5.5511e-17,   0.0000e+00,   1.0000e+00]])"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U.T.dot(U)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  1.0000e+00,   3.8858e-16,   5.5511e-17],\n",
              "       [  3.8858e-16,   1.0000e+00,   1.9429e-16],\n",
              "       [  5.5511e-17,   1.9429e-16,   1.0000e+00]])"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "V.T.dot(V)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now build the matrix $\\Sigma$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 2.4741,  0.    ,  0.    ],\n",
              "       [ 0.    ,  1.3509,  0.    ],\n",
              "       [ 0.    ,  0.    ,  0.4768],\n",
              "       [ 0.    ,  0.    ,  0.    ],\n",
              "       [ 0.    ,  0.    ,  0.    ]])"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Sigma = np.zeros(A.shape)\n",
        "Sigma[:3, :3] = np.diag(singval)\n",
        "Sigma"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now piece $A$ back together from the factorization:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -1.6653e-16,  -2.0817e-16,   0.0000e+00],\n",
              "       [ -2.2204e-16,   1.2837e-16,   3.1919e-16],\n",
              "       [ -8.3267e-17,  -2.2204e-16,  -2.4980e-16],\n",
              "       [  6.6613e-16,  -3.0531e-16,  -4.4409e-16],\n",
              "       [ -2.2204e-16,   1.3878e-16,   1.6653e-16]])"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U.dot(Sigma).dot(V.T) - A"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "----------------\n",
        "Next, compute the pseudoinverse:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.4042,  0.    ,  0.    ,  0.    ,  0.    ],\n",
              "       [ 0.    ,  0.7402,  0.    ,  0.    ,  0.    ],\n",
              "       [ 0.    ,  0.    ,  2.0971,  0.    ,  0.    ]])"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "SigmaInv = np.zeros((3,5))\n",
        "SigmaInv[:3, :3] = np.diag(1/singval)\n",
        "SigmaInv"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A_pinv = V.dot(SigmaInv).dot(U.T)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-------------\n",
        "Now use the pseudoinverse to \"solve\" $Ax=b$ for our tall-and-skinny $A$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "b = np.random.randn(5)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0.2425, -3.5347,  0.9537])"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "A_pinv.dot(b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---------------\n",
        "Compare with the solution from QR-based Least Squares:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0.2425, -3.5347,  0.9537])"
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Q, R = la.qr(A, \"complete\")\n",
        "la.solve(R[:3], Q.T.dot(b)[:3])"
      ]
    }
  ],
  "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 0
}