{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Normal Equations vs Pseudoinverse"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here's a simple overdetermined linear system, which we'll solve using both the normal equations and the pseudoinverse:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false,
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "A = np.random.randn(5, 3)\n",
        "b = np.random.randn(5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Normal Equations\n",
        "\n",
        "Solve $Ax\\cong b$ using the normal equations:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0.15751782,  0.37883556,  1.02377709])"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "x1 = la.solve(A.T@A, A.T@b)\n",
        "x1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Pseudoinverse\n",
        "\n",
        "Solve $Ax\\cong b$ using the pseudoinverse:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[ 0.40189265  0.43776205 -0.32197721 -0.68693466 -0.26701711]\n",
            " [-0.60148796 -0.48057566 -0.36776893 -0.32823601 -0.40529791]\n",
            " [-0.50534567  0.25317047  0.34332919 -0.50930065  0.55069808]\n",
            " [-0.46931436  0.7148832  -0.11893232  0.37004742 -0.34293739]\n",
            " [-0.03262382  0.04751987 -0.79313105  0.15511173  0.5861408 ]]\n",
            "[ 2.48597339  2.0073724   1.09629545]\n",
            "[[-0.16519593 -0.85534195  0.49101981]\n",
            " [ 0.90601342 -0.3283188  -0.26710758]\n",
            " [ 0.38967935  0.40074545  0.82918821]]\n"
          ]
        }
      ],
      "source": [
        "U, sigma, VT = la.svd(A)\n",
        "print(U)\n",
        "print(sigma)\n",
        "print(VT)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.40225692,  0.        ,  0.        ,  0.        ,  0.        ],\n",
              "       [ 0.        ,  0.49816367,  0.        ,  0.        ,  0.        ],\n",
              "       [ 0.        ,  0.        ,  0.91216287,  0.        ,  0.        ]])"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Sigma_inv = np.zeros_like(A.T)\n",
        "Sigma_inv[:3,:3] = np.diag(1/sigma)\n",
        "Sigma_inv"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0.15751782,  0.37883556,  1.02377709])"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pinv = VT.T @ Sigma_inv @ U.T\n",
        "x2 = pinv @ b\n",
        "x2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "4.7428748402675471e-16"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(x1-x2)"
      ]
    },
    {
      "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.5.1+"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}