{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Least Squares using the SVD"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "import scipy.linalg as spla"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# tall and skinny w/nullspace\n",
        "np.random.seed(12)\n",
        "A = np.random.randn(6, 4)\n",
        "b = np.random.randn(6)\n",
        "A[3] = A[4] + A[5]\n",
        "A[1] = A[5] + A[1]\n",
        "A[2] = A[3] + A[1]\n",
        "A[0] = A[3] + A[1]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Part I: Singular least squares using QR\n",
        "\n",
        "Let's see how successfully we can solve the least squares problem **when the matrix has a nullspace** using QR:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Q, R = la.qr(A)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[-4.526,  3.492, -0.204, -3.647],\n",
              "       [ 0.   ,  0.796,  0.034,  0.603],\n",
              "       [ 0.   ,  0.   , -1.459,  0.674],\n",
              "       [ 0.   , -0.   , -0.   ,  0.   ]])"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "R.round(3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can choose `x_qr[3]` as we please:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x_qr = np.zeros(A.shape[1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x_qr[3] = 0"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "QTbnew = Q.T.dot(b)[:3,] - R[:3, 3] * x_qr[3]\n",
        "x_qr[:3] = spla.solve_triangular(R[:3,:3], QTbnew, lower=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's take a look at the residual norm and the norm of `x_qr`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ -4.44089210e-16,   0.00000000e+00,  -1.11022302e-16,\n",
              "        -7.35042876e-01])"
            ]
          },
          "execution_count": 22,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "R.dot(x_qr)-Q.T.dot(b)[:4]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "2.1267152888030982"
            ]
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(A.dot(x_qr)-b, 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.82393512974131577"
            ]
          },
          "execution_count": 24,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(x_qr, 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Choose a different `x_qr[3]` and compare residual and norm of `x_qr`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "--------------\n",
        "### Part II: Solving least squares using the SVD\n",
        "Now compute the SVD of $A$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "U, sigma, VT = la.svd(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a matrix `Sigma` of the correct size:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Sigma = np.zeros(A.shape)\n",
        "Sigma[:4,:4] = np.diag(sigma)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And check that we've actually factorized `A`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0., -0.,  0.,  0.],\n",
              "       [ 0., -0.,  0.,  0.],\n",
              "       [ 0., -0.,  0.,  0.],\n",
              "       [ 0., -0., -0.,  0.],\n",
              "       [ 0., -0., -0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.]])"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "(U.dot(Sigma).dot(VT) - A).round(4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now define `Sigma_pinv` as the \"pseudo-\"inverse of `Sigma`, where \"pseudo\" means \"don't divide by zero\":"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.147,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  0.624,  0.   ,  0.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  0.   ,  1.055,  0.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ]])"
            ]
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Sigma_pinv = np.zeros(A.shape).T\n",
        "Sigma_pinv[:3,:3] = np.diag(1/sigma[:3])\n",
        "Sigma_pinv.round(3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compute the SVD-based solution for the least-squares problem:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x_svd = VT.T.dot(Sigma_pinv).dot(U.T).dot(b)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "2.1267152888030982"
            ]
          },
          "execution_count": 30,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(A.dot(x_svd)-b, 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.77354943014895838"
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(x_svd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What do you observe about $\\|\\text{x_svd}\\|_2$ compared to $\\|\\text{x_qr}\\|_2$?\n",
        "* Is $\\|\\text{x_svd}\\|_2$ compared to $\\|\\text{x_qr}\\|_2$?"
      ]
    }
  ],
  "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 0
}