{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear least squares "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(action=\"ignore\", module=\"scipy\", message=\"^internal gelsd\")\n",
    "import numpy as np\n",
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normal equations\n",
    "\n",
    "- Suppose $Ax = b$ is inconsistent. In this case, we can look instead for $\\widehat{x}$ which minimizes the distance between $Ax$ and $b$. In other words, we need to minimize $\\Vert Ax - b \\Vert^2$.\n",
    "- The minimum will occur when $\\langle Ax-b, Ax \\rangle = 0$ \n",
    "- Solving $(Ax)^T(Ax - b) = 0$ gives the normal equations $\\widehat{x} = (A^TA)^{-1}A^T b$\n",
    "- The corresponding vector in $C(A)$ is $A\\widehat{x} = A(A^TA)^{-1}A^T b$\n",
    "- Note that $P_A = A(A^TA)^{-1}A^T$ is the projection matrix for $C(A)$\n",
    "- This makes sense, since any solution to $Ax = b$ must be in $C(A)$, and the nearest such vector to $b$ must be the projection of $b$ onto $C(A)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 100\n",
    "n = 5\n",
    "A = np.random.rand(m, n)\n",
    "b = np.random.rand(m, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using least squares function (SVD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, res, rank, s, = la.lstsq(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.31192056],\n",
       "       [0.15414167],\n",
       "       [0.18783713],\n",
       "       [0.07997706],\n",
       "       [0.17878726]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A @ x).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using normal equations (projection) - for illustration only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.31192056],\n",
       "       [0.15414167],\n",
       "       [0.18783713],\n",
       "       [0.07997706],\n",
       "       [0.17878726]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "la.inv(A.T @ A) @ A.T @ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection matrices\n",
    "\n",
    "- Let $P$ be the projection matrix onto $C(A)$\n",
    "- Since it is a projection operator, $P^2$ = $P$ (check)\n",
    "- Check that $I-P$ is also idempotent (it turns out to also be a projection matrix)\n",
    "- Where does $I-P$ project to?\n",
    "- Show $C(I-P) = N(P)$\n",
    "- Show $\\rho(P) + \\nu(P) = n$\n",
    "- This implies $\\mathbb{R}^n = C(P) \\oplus C(I-P)$\n",
    "- Trivial example\n",
    "$$\n",
    "P = \\begin{bmatrix}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "0 & 0 & 0\n",
    "\\end{bmatrix}, (I-P) = \\begin{bmatrix}\n",
    "0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "- Geometry of $P$ and $I-P$ in linear least squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "P = A @ la.inv(A.T @ A) @ A.T "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 100)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x, res, rank, s, = la.lstsq(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q = np.eye(m) - P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.2763989 ],\n",
       "       [0.33523772],\n",
       "       [0.4619637 ],\n",
       "       [0.50292904],\n",
       "       [0.49148152]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(P @ b)[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.2763989 ],\n",
       "       [0.33523772],\n",
       "       [0.4619637 ],\n",
       "       [0.50292904],\n",
       "       [0.49148152]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A @ x)[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.834201864805834"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "la.norm(Q @ b)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9.83420186])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}