{"nbformat_minor": 0, "metadata": {"name": "", "signature": "sha256:eb6b2710f0d7a82b6076037280b773f7bc653538e39c69fd8d4b0b9d977ee1c6"}, "worksheets": [{"metadata": {}, "cells": [{"metadata": {}, "cell_type": "markdown", "source": ["# Pseudoinverse and Least Squares"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 1, "outputs": [], "collapsed": false, "input": ["import numpy as np\n", "import numpy.linalg as la\n", "\n", "np.set_printoptions(precision=4, linewidth=100)"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 2, "outputs": [], "collapsed": false, "input": ["A = np.random.randn(5, 3)"]}, {"metadata": {}, "cell_type": "markdown", "source": ["Now compute the SVD of $A$. Note that `numpy.linalg.svd` returns $V^T$:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 3, "outputs": [], "collapsed": false, "input": ["U, singval, VT = la.svd(A)\n", "V = VT.T"]}, {"metadata": {}, "cell_type": "markdown", "source": ["Let's first understand the shapes of these arrays:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 4, "outputs": [{"stream": "stdout", "output_type": "stream", "text": ["(5, 5)\n", "(3,)\n", "(3, 3)\n"]}], "collapsed": false, "input": ["print(U.shape)\n", "print(singval.shape)\n", "print(V.shape) "]}, {"metadata": {}, "cell_type": "markdown", "source": ["Check the orthogonality of $U$ and $V$:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 5, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 5, "text": ["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]])"]}], "collapsed": false, "input": ["U.T.dot(U)"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 6, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": ["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]])"]}], "collapsed": false, "input": ["V.T.dot(V)"]}, {"metadata": {}, "cell_type": "markdown", "source": ["Now build the matrix $\\Sigma$:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 7, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": ["array([[ 2.4741,  0.    ,  0.    ],\n", "       [ 0.    ,  1.3509,  0.    ],\n", "       [ 0.    ,  0.    ,  0.4768],\n", "       [ 0.    ,  0.    ,  0.    ],\n", "       [ 0.    ,  0.    ,  0.    ]])"]}], "collapsed": false, "input": ["Sigma = np.zeros(A.shape)\n", "Sigma[:3, :3] = np.diag(singval)\n", "Sigma"]}, {"metadata": {}, "cell_type": "markdown", "source": ["Now piece $A$ back together from the factorization:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 8, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": ["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]])"]}], "collapsed": false, "input": ["U.dot(Sigma).dot(V.T) - A"]}, {"metadata": {}, "cell_type": "markdown", "source": ["----------------\n", "Next, compute the pseudoinverse:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 9, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": ["array([[ 0.4042,  0.    ,  0.    ,  0.    ,  0.    ],\n", "       [ 0.    ,  0.7402,  0.    ,  0.    ,  0.    ],\n", "       [ 0.    ,  0.    ,  2.0971,  0.    ,  0.    ]])"]}], "collapsed": false, "input": ["SigmaInv = np.zeros((3,5))\n", "SigmaInv[:3, :3] = np.diag(1/singval)\n", "SigmaInv"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 10, "outputs": [], "collapsed": false, "input": ["A_pinv = V.dot(SigmaInv).dot(U.T)"]}, {"metadata": {}, "cell_type": "markdown", "source": ["-------------\n", "Now use the pseudoinverse to \"solve\" $Ax=b$ for our tall-and-skinny $A$:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 11, "outputs": [], "collapsed": false, "input": ["b = np.random.randn(5)"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 12, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": ["array([ 0.2425, -3.5347,  0.9537])"]}], "collapsed": false, "input": ["A_pinv.dot(b)"]}, {"metadata": {}, "cell_type": "markdown", "source": ["---------------\n", "Compare with the solution from QR-based Least Squares:"]}, {"language": "python", "metadata": {}, "cell_type": "code", "prompt_number": 13, "outputs": [{"metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": ["array([ 0.2425, -3.5347,  0.9537])"]}], "collapsed": false, "input": ["Q, R = la.qr(A, \"complete\")\n", "la.solve(R[:3], Q.T.dot(b)[:3])"]}]}], "nbformat": 3}