16_LA_SVD_and_apps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli and Rajath Kumar Mysore Pradeep Kumar</td>\n",
    "</table>\n",
    "\n",
    "Note:  This material largely follows the text \"Numerical Linear Algebra\" by Trefethen and Bau (SIAM, 1997) and is meant as a guide and supplement to the material presented there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%precision 6 \n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Understanding the Singular Value Decomposition\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Factorizations:  A review\n",
    "\n",
    "**Matrix Factorization** is a fundamental concept in Linear Algebra and refers to the ability to write general matrices as the products of simpler matrices with special properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In particular, all of the great algorithms encountered in numerical linear Algebra can be succinctly described by their factorizations.  Examples include"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### **The LU decomposition**  \n",
    "$$PA=LU$$\n",
    "    \n",
    "* **Algorithms** *Gaussian Elimination with partial pivoting* \n",
    "* **Factorization** $P$ Permutation matrix (from pivoting), $L,U$: Lower and upper triangular matrices\n",
    "* **Application** direct solution of $A\\mathbf{x}=\\mathbf{b}$ for $A\\in\\mathbb{R}^{n\\times n}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### **The QR decomposition**  \n",
    "$$A=QR$$\n",
    "    \n",
    "* **Algorithms** Gram-Schmidt Orthogonalization, Householder Triangularization \n",
    "* **Factorization**  $Q$:  Orthonormal Basis for $C(A)$, $R$: upper triangular matrix\n",
    "* **Applications** Least-Squares solutions, Projection problems,  Eigenvalues $QR/RQ$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### **EigenProblems**  \n",
    "$$ AX = X\\Lambda, \\quad (A\\in\\mathbb{R}^{n\\times n})$$\n",
    "\n",
    "* Diagonalizable matrices $A = X\\Lambda X^{-1}$\n",
    "* Hermitian/Symmetric Matrices $A = Q\\Lambda Q^T$\n",
    "* Schur Factorization $A = QTQ^T$\n",
    "    \n",
    "    \n",
    "* **Algorithms**: Power Method, Inverse Power with shifts,  $QR/RQ$\n",
    "* **Factorization**-- $X$, $Q$: matrix of Eigenvectors,  $\\Lambda$ and diagonal matrix of eigenvalues \n",
    "* **Application** Solution of Dynamical systems,  Iterative Maps, Markov Chains,  Vibrational analysis$\\ldots$, Quantum Mechanics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "*But perhaps the most beautiful factorization, which contains aspects of all of these problems$\\ldots$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Singular Value Decomposition (SVD)\n",
    "$$\n",
    "    A = U\\Sigma V^{T}\n",
    "$$\n",
    "for *any* $A\\in\\mathbb{R}^{m\\times n}$, where\n",
    "* $U \\in \\mathbb R^{m \\times m}$ and is the orthogonal matrix whose columns are the eigenvectors of $AA^{T}$\n",
    "* $V \\in \\mathbb R^{n \\times n}$ and is the orthogonal matrix whose columns are the eigenvectors of $A^{T}A$\n",
    "* $\\Sigma \\in \\mathbb R^{m \\times n}$ and is a diagonal matrix with elements $\\sigma_{1}, \\sigma_{2}, \\sigma_{3}, ... \\sigma_{r}$ where $r = rank(A)$ corresponding to the square roots of the eigenvalues of $A^{T}A$. They are called the singular values of $A$ and are positive arranged in descending order. ($\\sigma_{1} \\geq \\sigma_{2} \\geq \\sigma_{3} \\geq ... \\sigma_{r} > 0$).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Or a picture is worth a thousand words...\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Existence and Uniqueness\n",
    "\n",
    "Every matrix $A \\in \\mathbb{C}^{m \\times n}$ has a singular value decomposition. Furthermore, the singular values $\\{\\sigma_{j}\\}$ are uniquely determined, and if $A$ is square and the $\\sigma_{j}$ are distinct, the left and right singular vectors $\\{u_{j}\\}$ and $\\{v_{j}\\}$ are uniquely determined up to complex signs (i.e., complex scalar factors of absolute value 1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Eigenvalue Decomposition vs. SVD Decomposition\n",
    "\n",
    "Let the matrix $X$ contain the eigenvectors of $A$ which are linearly independent, then we can write a decomposition of the matrix $A$ as\n",
    "$$\n",
    "    A = X \\Lambda X^{-1}.\n",
    "$$\n",
    "\n",
    "How does this differ from the SVD?\n",
    " - The basis of the SVD representation differs from the eigenvalue decomposition\n",
    "   - The basis vectors are not in general orthogonal for the eigenvalue decomposition where it is for the SVD\n",
    "   - The SVD effectively contains two basis sets.\n",
    " - All matrices have an SVD decomposition whereas not all have eigenvalue decompositions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Applications\n",
    "* Matrix Pseudo-Inverse for ill-conditions problems\n",
    "* Principal Component analysis\n",
    "* Total Least Squares\n",
    "* Image Compression\n",
    "* Data Analysis for interpretation of high-dimensional data\n",
    "* and more$\\ldots$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A conceptual computation of the SVD (not how it's really done)\n",
    "\n",
    "Given $A = U\\Sigma V^T$, we can form the matrix \n",
    "\n",
    "$$\n",
    "    A^TA = V\\Sigma^T\\Sigma V^T\n",
    "$$\n",
    "\n",
    "where $V\\in\\mathbb{R}^{n\\times n}$ is unitary and\n",
    "$$\n",
    "    \\Sigma^T\\Sigma = \\begin{bmatrix} \n",
    "                \\sigma_1^2 & & &    \\\\\n",
    "                    & \\sigma_2^2 & & \\\\\n",
    "                   & & \\ddots &   \\\\\n",
    "                   & & &   \\sigma_n^2  \\\\ \n",
    "                                   \\end{bmatrix}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "However, because $A^TA$ is clearly **Symmetric** it can also has the eigen factorization\n",
    "\n",
    "$$\n",
    "    A^TA = Q\\Lambda Q^T\n",
    "$$\n",
    "\n",
    "where $Q$ is unitary and $\\Lambda$ is real diagonal.  \n",
    "\n",
    "Moreover, we can show that $A^T A$ is at least Positive semi-definite so all the eigenvalues are $\\geq 0$. \n",
    "\n",
    "In particular,  if $\\mathrm{rank}(A)=r$ we know that $r$ eigenvalues are $>0$ and $n-r$ eigenvalues are equal to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given that\n",
    "\n",
    "$$\\begin{align}\n",
    "    A^TA &= V\\Sigma^T\\Sigma V^T \\\\\n",
    "        &= Q\\Lambda Q^T \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "If we simply order the eigenvalues from largest to smallest (and rearrange the columns of $Q$ to match), then we can make the association\n",
    "\n",
    "$$V = Q\\quad, \\Sigma^T\\Sigma = \\Lambda$$ or\n",
    "\n",
    "$$\\sigma_i^2 = \\lambda_i\\quad\\textrm{or}\\quad \\sigma_i = \\sqrt{\\lambda_i}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A similar argument can be made about $U$ and the eigenvectors of $AA^T$,  however, a better way to think about the relationship between $U$, $V$ and $\\Sigma$ is to rearrange the SVD as\n",
    "\n",
    "$$ \n",
    "    AV = U\\Sigma\n",
    "$$\n",
    "\n",
    "or we can look at this column-by-column \n",
    "\n",
    "$$\n",
    "    A\\mathbf{v}_i = \\sigma_i\\mathbf{u}_i\n",
    "$$\n",
    "\n",
    "(in the same way the Eigen decomposition $AX = X\\Lambda$ is equivalent to $A\\mathbf{x}_i = \\lambda_i\\mathbf{x}_i$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or rearranging \n",
    "$$\n",
    "    \\mathbf{u}_i = \\frac{A\\mathbf{v}_i}{\\sigma_i},\\quad i=1,2,\\ldots,r\n",
    "$$\n",
    "\n",
    "lets us solve for the first $r$ columns of $U$ given $V$ and $\\Sigma$.  \n",
    "\n",
    "As we'll show, we usually don't need to solve for the remaing $m-r$ columns of $U$ for reasons that relate to the \"4 fundamental subspaces of $A$\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Review:  the Singular Value Decomposition (SVD)\n",
    "$$\n",
    "    A = U\\Sigma V^{T}\n",
    "$$\n",
    "for *any* $A\\in\\mathbb{R}^{m\\times n}$, where\n",
    "* $U \\in \\mathbb R^{m \\times m}$ and is the orthogonal matrix whose columns are the eigenvectors of $AA^{T}$\n",
    "* $V \\in \\mathbb R^{n \\times n}$ and is the orthogonal matrix whose columns are the eigenvectors of $A^{T}A$\n",
    "* $\\Sigma \\in \\mathbb R^{m \\times n}$ and is a diagonal matrix with elements $\\sigma_{1}, \\sigma_{2}, \\sigma_{3}, ... \\sigma_{r}$ where $r = rank(A)$ corresponding to the square roots of the eigenvalues of $A^{T}A$. They are called the singular values of $A$ and are positive arranged in descending order. ($\\sigma_{1} \\geq \\sigma_{2} \\geq \\sigma_{3} \\geq ... \\sigma_{r} > 0$).\n",
    "\n",
    "But a better way to write and understand the SVD is as\n",
    "$$ \n",
    "    AV = U\\Sigma\n",
    "$$\n",
    "where the columns of $U$ and $V$ contain orthonormal bases for the 4 fundamental subspaces of $A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The Strangian view of the 4 Subspaces: a quick refresher\n",
    "\n",
    "<img align=center src=\"./images/Strang_4_subspaces.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Big idea:  The columns of $U$ and $V$ contain orthonormal bases for the 4 fundamental subspaces\n",
    "\n",
    "Returning to $AV = U\\Sigma$ it follows that \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "        A\\mathbf{v}_i &= \\sigma_i\\mathbf{u_i},\\quad\\mathrm{for}\\quad i=1,2\\ldots r \\\\\n",
    "        A\\mathbf{v}_i &= \\mathbf{0},\\quad\\mathrm{for}\\quad i=r+1,\\ldots n \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Therefore the last $n-r$ columns of $V$ must be in the Null space of $A$, and in fact must form an orthonormal basis for the Null space $N(A)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Since the first $r$ columns of $V$ are all orthogonal to the last $n-r$ columns,  they must form an orthonormal basis for the row space $C(A^T)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Likewise,  since the first $r$ columns of $U$ satisfy\n",
    "\n",
    "$$\n",
    "    \\mathbf{u}_i = \\frac{A\\mathbf{v}_i}{\\sigma_i},\\quad i=1,2,\\ldots,r\n",
    "$$\n",
    "\n",
    "they must form an orthonormal basis for $C(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Leaving only the last $m-r$ columns of $U$ as an orthonormal basis for the left null space $N(A^T)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again a picture is worth a lot here\n",
    "\n",
    "$$\n",
    "A\\begin{bmatrix} &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " \\mathbf{v}_1 & \\cdots & \\mathbf{v}_r & | &\\mathbf{v}_{r+1} & \\cdots & \\mathbf{v}_n \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    "  &  &  & |&  &  &  \\\\\n",
    "   \\end{bmatrix} = \\begin{bmatrix} &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " \\mathbf{u}_1 & \\cdots & \\mathbf{u}_r & | &\\mathbf{u}_{r+1} & \\cdots & \\mathbf{u}_m \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    "  &  &  & |&  &  &  \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix} \\sigma_1 &  &  & |  &  &    \\\\\n",
    " & \\ddots &  & |   &   \\mathbf{0}  \\\\\n",
    " &  & \\sigma_r &|    &     \\\\\n",
    " -&- &-  & | & - &  -\\\\\n",
    " & \\mathbf{0} &  & |   & \\mathbf{0}   \\\\\n",
    " &  &  & |  &    \\\\\n",
    "  &  &  & |  &     \\\\\n",
    "    \\end{bmatrix} \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## The Economy (or Skinny) SVD\n",
    "\n",
    "As it turns out, because of all the 0's  all we actually need to reconstruct $A$ is the the first $r$ columns of $U$, and $V$,  and the square sub-block of $\\Sigma$ with just the singular values.  i.e.\n",
    "\n",
    "$$\n",
    "A\\begin{bmatrix} &  &  & \\\\\n",
    " &  &  &    \\\\\n",
    " \\mathbf{v}_1 & \\cdots & \\mathbf{v}_r   \\\\\n",
    " &  &  &   \\\\\n",
    "  &  &  &    \\\\\n",
    "   \\end{bmatrix} = \\begin{bmatrix} &  &  &   \\\\\n",
    " &  &  &    \\\\\n",
    "  &  &  &    \\\\ \n",
    " \\mathbf{u}_1 & \\cdots & \\mathbf{u}_r     \\\\\n",
    " &  &  &    \\\\\n",
    " &  &  &    \\\\\n",
    "  &  &  &    \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix} \\sigma_1 &  &        \\\\\n",
    " & \\ddots &    \\\\\n",
    " &  & \\sigma_r       \\\\\n",
    "    \\end{bmatrix} \n",
    "$$\n",
    "or"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    A = U_r\\Sigma_r V_r^T\n",
    "$$\n",
    "\n",
    "And it is this object (The Economy (or Skinny) SVD), that helps us understand what the SVD does."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Spectral Theorem\n",
    "\n",
    "The Economy SVD can also be written as \n",
    "$$\n",
    "    A = U_r\\Sigma_r V_r^T = \\sum_{i=1}^r \\sigma_i \\mathbf{u}_i\\mathbf{v}^T_i\n",
    "$$\n",
    "\n",
    "Which says we can expand $A$ as a series of rank-1 matrices $\\mathbf{u}_i\\mathbf{v}^T_i$ weighted by the singular values...\n",
    "\n",
    "and it is this picture that leads to many of the important approximating properties of the SVD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Matrix Multiplication \n",
    "\n",
    "Consider the Matrix Vector product $A\\mathbf{x}$ which maps a vector $\\mathbf{x}\\in\\mathbb{R}^n$ to its image in column space $C(A)$.  In the context of the Economy SVD\n",
    "$$\n",
    "    A\\mathbf{x} = U_r\\Sigma_r V_r^T\\mathbf{x}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or (since $V_r^TV_r = I$)\n",
    "$$ \n",
    "    A\\mathbf{x} = U_r\\Sigma_r (V_r^T V_r) V_r^T\\mathbf{x} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "$$\n",
    "\\begin{align}\n",
    "    A\\mathbf{x} &=  (U_r\\Sigma_r V_r^T) (V_r V_r^T)\\mathbf{x} \\\\\n",
    "    & = A \\mathbf{x}^+ \n",
    "\\end{align}\n",
    "$$\n",
    "where $$\\mathbf{x}^+ = (V_r V_r^T)\\mathbf{x}$$\n",
    "is the projection of $\\mathbf{x}$ onto the row space $C(A^T)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Question**: If we wanted it, how would we  find the projection of $\\mathbf{x}$ onto the Null space of $A$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### SVD and the Matrix Inverse \n",
    "\n",
    "If a matrix is square and invertible, it implies that $m=n=r$ and $U$, $\\Sigma$ and $V$ are all square invertible matrices,  therefore if \n",
    "\n",
    "$$\n",
    "    A = U\\Sigma V^T\n",
    "$$\n",
    "then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "    A^{-1} = V\\Sigma^{-1} U^T\n",
    "$$\n",
    "where \n",
    "$$\n",
    "    \\Sigma^{-1} = \\begin{bmatrix} 1/\\sigma_1 &  &    &    \\\\\n",
    " & 1/\\sigma_2 &   & \\\\\n",
    "  &   &  \\ddots &  \\\\\n",
    " &  &  & 1/\\sigma_n       \\\\\n",
    "    \\end{bmatrix} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and clearly \n",
    "$$\n",
    "    A^{-1}\\mathbf{b} = \\mathbf{x}\n",
    "$$ \n",
    "maps any vector $\\mathbf{b}\\in C(A)$ back to a unique vector $\\mathbf{x}\\in C(A^T)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### SVD and the Matrix Pseudo-Inverse \n",
    "\n",
    "But suppose $A$ is not square.  Then it cannot have an inverse,  however, the SVD allows us to define a \"Pseudo-Inverse\"\n",
    "\n",
    "Given, the skinny SVD\n",
    "\n",
    "$$\n",
    "    A = U_r\\Sigma_r V_r^T\n",
    "$$\n",
    "\n",
    "we can define\n",
    "\n",
    "$$\n",
    "    A^{+} = V_r \\Sigma_r^{-1} U_r^T\n",
    "$$\n",
    "\n",
    "because $\\Sigma_r$ is square and invertible (but its size is set by the rank $r$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Action of the  Pseudo-Inverse \n",
    "\n",
    "In general, if $A$ is rank deficient ($r<\\min(m,n)$),  the problem $A\\mathbf{x} = \\mathbf{b}$ has either no solution or an infinite number of solutions.  However, it always has a minimal least squares solution given by \n",
    "\n",
    "$$\n",
    "    x^{+} = A^{+}\\mathbf{b}\n",
    "$$\n",
    "\n",
    "which maps any vector $\\mathbf{b}\\in \\mathbb{R}^m$ to a unique vector $\\mathbf{x}^+\\in C(A^T)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As we did with $A\\mathbf{x}$, we can consider the action of because $A^{+}$ on any vector $\\mathbf{b}$ as a projection problem.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    A^{+}\\mathbf{b} &= V_r\\Sigma^{-1}_r U_r^T\\mathbf{b} \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "  \\hspace{5em}  &= V_r\\Sigma^{-1}_r (U_r^T U_r) U_r^T\\mathbf{b} \\\\\n",
    "    & = (V_r\\Sigma^{-1}_r U_r^T) (U_r U_r^T)\\mathbf{b} \\\\\n",
    "    & = A^{+} \\mathbf{p}  \n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "where $$\\mathbf{p}  = (U_r U_r^T)\\mathbf{b}$$\n",
    "is projection of $\\mathbf{b}$ onto $C(A)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The minimal Least-Squares solution\n",
    "\n",
    "It should be clear that \n",
    "$$\n",
    "    \\mathbf{x}^+ = A^{+}\\mathbf{b}\n",
    "$$ \n",
    "lies entirely in the Row space of $A$ (Why?)\n",
    "\n",
    "Therefore it has no component in the $N(A)$ and is the shortest solution to $A\\hat{\\mathbf{x}}=\\mathbf{p}$ which is the least squares problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Pseudo-Inverse and ill-conditioned problems\n",
    "\n",
    "Technically, a square matrix is either invertible or singular,  however, in many real problems, a matrix can be close to singular or \"ill-conditioned\".\n",
    "\n",
    "Again, the condition number of a matrix defined by an induced $p$ norm is\n",
    "$$\n",
    "    \\kappa_p(A) = ||A||_p ||A^{-1}||_p\n",
    "$$\n",
    "\n",
    "In particular,  for $p=2$, it is easy to show (using the SVD) that for a square matrix\n",
    "$$\n",
    "    \\kappa_2(A) =\\frac{\\sigma_1}{\\sigma_n}\n",
    "$$\n",
    "\n",
    "and thus for a singular matrix with $r<n$, $\\kappa_2 =\\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "###  ill-conditioned matrices\n",
    "\n",
    "However, it is possible that while $\\sigma_n>0$,  it is significantly smaller than $\\sigma_1$, leading to a very large condition number.  In particular if the last $j$ singular values are very small, it suggest that the matrix has a \"near Null Space\",  with a basis defined by the last $j$ columns of $V$.  Thus, while strictly invertible,  if a direct solver is used to solve\n",
    "\n",
    "$$\n",
    "A\\mathbf{x} =\\mathbf{b}\n",
    "$$ \n",
    "\n",
    "any component in the near null space will be amplified by $1/\\sigma_j$ and completely pollute your answer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "However, a useful fix can be to approximate $A$ as a lower rank matrix \n",
    "\n",
    "$$A_k = U_k\\Sigma_kV_k^T$$\n",
    "\n",
    "which just keeps the first $k$ singular values (i.e. pretends $r=k$).  If chosen wisely,  The approximate solution\n",
    "\n",
    "$$\n",
    "    \\mathbf{x}^{+} = A^{+}_k\\mathbf{b}\n",
    "$$ \n",
    "\n",
    "can be considerably more accurate as it does not allow the near-null space to contaminate the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Full SVD example\n",
    "\n",
    "Consider the matrix\n",
    "$$\n",
    "    A = \\begin{bmatrix} \n",
    "        2 & 0 & 3 \\\\\n",
    "        5 & 7 & 1 \\\\\n",
    "        0 & 6 & 2 \n",
    "    \\end{bmatrix}.\n",
    "$$\n",
    "Confirm the SVD representation using `numpy` functions as appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.array([\n",
    "    [2.0, 0.0, 3.0],\n",
    "    [5.0, 7.0, 1.0],\n",
    "    [0.0, 6.0, 2.0]\n",
    "])\n",
    "\n",
    "U, sigma, V_T = numpy.linalg.svd(A, full_matrices=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('U:\\n{}\\n'.format(U))\n",
    "print('Sigma:\\n{}\\n'.format(numpy.diag(sigma)))\n",
    "print('V:\\n{}\\n'.format(V_T.T))\n",
    "print(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Aprime = numpy.dot(U, numpy.dot(numpy.diag(sigma), V_T))\n",
    "print('A - USV^T:\\n{}'.format(A - Aprime))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: a rank-1 Matrix \n",
    "\n",
    "$$A = \\mathbf{x}\\mathbf{y}^T$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.array([ 1., 2., 1. ])\n",
    "y = numpy.array([2., -1., 1. ])\n",
    "A = numpy.outer(x,y)\n",
    "print('A:\\n{}\\n'.format(A))\n",
    "U, sigma, V_T = numpy.linalg.svd(A, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print('U:\\n{}\\n'.format(U))\n",
    "print('sigma:\\n{}\\n'.format(sigma))\n",
    "print('V:\\n{}\\n'.format(V_T.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Aprime = numpy.dot(U, numpy.dot(numpy.diag(sigma), V_T))\n",
    "print('A - USV^T:\\n{}'.format(A - Aprime))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# rank-1 reconstruction\n",
    "k = 1\n",
    "A1 = numpy.dot(U[:,:k], numpy.dot(numpy.diag(sigma[:k]), V_T[:k,:]))\n",
    "print('A_k =\\n{}\\n'.format(A1))\n",
    "print('A - A_k:\\n{}'.format(A - A1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### An ill-conditioned example\n",
    "\n",
    "Consider the matrix\n",
    "$$\n",
    "    A = \\begin{bmatrix} \n",
    "        1 & 0 & 1 \\\\\n",
    "        1 & 1 & 2+\\epsilon \\\\\n",
    "        0 & 1 & 1\\\\\n",
    "    \\end{bmatrix}.\n",
    "$$\n",
    "where the last column is almost the sum of the first two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "epsilon = 1.e-5\n",
    "A = numpy.array([\n",
    "    [1.0, 0.0, 1.0],\n",
    "    [1.0, 1.0, 2+epsilon],\n",
    "    [0.0, 1.0, 1.0]\n",
    "])\n",
    "\n",
    "U, S, V_T = numpy.linalg.svd(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('U:\\n{}\\n'.format(U))\n",
    "print('S:\\n{}\\n'.format(S))\n",
    "print('V:\\n{}\\n'.format(V_T.T))\n",
    "print('sigma_1/sigma_n = {}, cond_2(A) = {}'.format(S[0]/S[2], numpy.linalg.cond(A, p=2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Try to solve $A\\mathbf{x}=\\mathbf{b}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "epsilon = 5*numpy.finfo(numpy.float64).eps\n",
    "A = numpy.array([\n",
    "    [1.0, 0.0, 1.0],\n",
    "    [1.0, 1.0, 2+epsilon],\n",
    "    [0.0, 1.0, 1.0]\n",
    "])\n",
    "x_true = numpy.array([1., 2., 3.])\n",
    "b = A.dot(x_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linalg.solve(A,b)\n",
    "print('x = {}: cond(A)={}, S ={}'.format(x, numpy.linalg.cond(A, p=2), S))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Rank 2 reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "U, S, VT = numpy.linalg.svd(A)\n",
    "print(S)\n",
    "k = 2\n",
    "A2 = numpy.dot(U[:,:k], numpy.dot(numpy.diag(S[:k]), VT[:k,:]))\n",
    "#Pseudo Inverse\n",
    "Ap2 = numpy.dot(VT.T[:,:k], numpy.dot(numpy.diag(1./S[:k]), U.T[:k,:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "xp = Ap2.dot(b)\n",
    "print('x_true\\t= {}'.format(x_true))\n",
    "print('x^+ \\t= {}'.format(xp))\n",
    "print('x \\t= {}: '.format(numpy.linalg.solve(A,b)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Review:  the Singular Value Decomposition (SVD)\n",
    "$$\n",
    "    A = U\\Sigma V^{T}\n",
    "$$\n",
    "for *any* $A\\in\\mathbb{R}^{m\\times n}$, where\n",
    "* $U \\in \\mathbb R^{m \\times m}$ and is the orthogonal matrix whose columns are the eigenvectors of $AA^{T}$\n",
    "* $V \\in \\mathbb R^{n \\times n}$ and is the orthogonal matrix whose columns are the eigenvectors of $A^{T}A$\n",
    "* $\\Sigma \\in \\mathbb R^{m \\times n}$ and is a diagonal matrix with elements $\\sigma_{1}, \\sigma_{2}, \\sigma_{3}, ... \\sigma_{r}$ where $r = rank(A)$ corresponding to the square roots of the eigenvalues of $A^{T}A$. They are called the singular values of $A$ and are positive arranged in descending order. ($\\sigma_{1} \\geq \\sigma_{2} \\geq \\sigma_{3} \\geq ... \\sigma_{r} > 0$).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Again a picture is worth a lot here\n",
    "\n",
    "$$\n",
    "A\\begin{bmatrix} &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " \\mathbf{v}_1 & \\cdots & \\mathbf{v}_r & | &\\mathbf{v}_{r+1} & \\cdots & \\mathbf{v}_n \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    "  &  &  & |&  &  &  \\\\\n",
    "   \\end{bmatrix} = \\begin{bmatrix} &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " \\mathbf{u}_1 & \\cdots & \\mathbf{u}_r & | &\\mathbf{u}_{r+1} & \\cdots & \\mathbf{u}_m \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    " &  &  & |&  &  &  \\\\\n",
    "  &  &  & |&  &  &  \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix} \\sigma_1 &  &  & |  &  &    \\\\\n",
    " & \\ddots &  & |   &   \\mathbf{0}  \\\\\n",
    " &  & \\sigma_r &|    &     \\\\\n",
    " -&- &-  & | & - &  -\\\\\n",
    " & \\mathbf{0} &  & |   & \\mathbf{0}   \\\\\n",
    " &  &  & |  &    \\\\\n",
    "  &  &  & |  &     \\\\\n",
    "    \\end{bmatrix} \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### $U$ and $V$ contain orthonormal bases for the 4 subspaces of $A$\n",
    "\n",
    "* The first $r$ columns of $U$ form an orthonormal bases for $C(A)\\in\\mathbb{R}^m$\n",
    "* The first $r$ columns of $V$ form an orthonormal bases for $C(A^T)\\in\\mathbb{R}^n$\n",
    "* The last $n-r$ columns of $V$ form an orthonormal basis for $N(A)\\in\\mathbb{R}^n$\n",
    "* The last $m-r$ columns of $U$ form an orthonormal basis for $N(A^T)\\in\\mathbb{R}^m$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## The Economy (or Skinny) SVD\n",
    "\n",
    "As it turns out, because of all the 0's  all we actually need to reconstruct $A$ is the the first $r$ columns of $U$, and $V$,  and the square sub-block of $\\Sigma$ with just the singular values.  i.e.\n",
    "\n",
    "$$\n",
    "A\\begin{bmatrix} &  &  & \\\\\n",
    " &  &  &    \\\\\n",
    " \\mathbf{v}_1 & \\cdots & \\mathbf{v}_r   \\\\\n",
    " &  &  &   \\\\\n",
    "  &  &  &    \\\\\n",
    "   \\end{bmatrix} = \\begin{bmatrix} &  &  &   \\\\\n",
    " &  &  &    \\\\\n",
    "  &  &  &    \\\\ \n",
    " \\mathbf{u}_1 & \\cdots & \\mathbf{u}_r     \\\\\n",
    " &  &  &    \\\\\n",
    " &  &  &    \\\\\n",
    "  &  &  &    \\\\\n",
    "    \\end{bmatrix} \\begin{bmatrix} \\sigma_1 &  &        \\\\\n",
    " & \\ddots &    \\\\\n",
    " &  & \\sigma_r       \\\\\n",
    "    \\end{bmatrix} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* $V_rV_r^T$ is an orthogonal projector onto $C(A^T)$ \n",
    "* $(I - V_rV_r^T)$ projects onto $N(A)$\n",
    "* $U_rU_r^T$ is an orthogonal projector onto $C(A)$\n",
    "* $(I - U_rU_r^T)$ projects onto $N(A^T)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## The Pseudo-Inverse\n",
    "\n",
    "Given the SVD,  every matrix $A\\in\\mathbb{R}^{m\\times n}$  has a 'pseudo-inverse'\n",
    "\n",
    "$$\n",
    "    A^{+} = V_r\\Sigma_r^{-1}U_r^T\n",
    "$$\n",
    "\n",
    "with the properties\n",
    "* if $A\\in\\mathbb{R}^{n\\times n}$ and $r=n$, then $A^{+}=A^{-1}$\n",
    "* $\\mathbf{x}^+ = A^{+}\\mathbf{b}$ is always the shortest least squares solution s.t. $(I - V_rV_r^T)\\mathbf{x}^+=\\mathbf{0}$\n",
    "* $AA^+$ projects $\\mathbf{b}$ onto ______?\n",
    "* $A^+A$ projects $\\mathbf{x}$ onto ______?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Matrix Properties via the SVD\n",
    "\n",
    " - The $\\text{rank}(A) = r$ where $r$ is the number of non-zero singular values.\n",
    " - The $\\text{range}(A) = \\mathrm{span}\\langle u_1, ... , u_r\\rangle$ and $\\text{null}(a) = \\mathrm{span}\\langle v_{r+1}, ... , v_n\\rangle$.\n",
    " - The $|| A ||_2 = \\sigma_1$ and $||A||_F = \\sqrt{\\sigma_{1}^{2}+\\sigma_{2}^{2}+...+\\sigma_{r}^{2}}$.\n",
    " - The nonzero singular values of A are the square roots of the nonzero eigenvalues of $A^{T}A$ or $AA^{T}$.\n",
    " - If $A = A^{T}$, then the singular values of $A$ are the absolute values of the eigenvalues of $A$.\n",
    " - For $A \\in \\mathbb{C}^{m \\times m}$ then $|det(A)| = \\Pi_{i=1}^{m} \\sigma_{i}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Low-Rank Approximations\n",
    "\n",
    " - $A$ is the sum of the $r$ rank-one matrices:\n",
    "$$\n",
    "    A = U \\Sigma V^T = \\sum_{j=1}^{r} \\sigma_{j}u_{j}v_{j}^{T}\n",
    "$$\n",
    " - For any $k$ with $0 \\leq k \\leq r$, define\n",
    "$$\n",
    "    A_{k} = \\sum_{j=1}^{k} \\sigma_{j}u_{j}v_{j}^{T}\n",
    "$$\n",
    "\n",
    "then \n",
    "\n",
    "$$\n",
    "     ||A - A_{k}||_{2} = ||\\sum_{j=k+1}^{r} \\sigma_{j}u_{j}v_{j}^{T}||_{2} = \\sigma_{k+1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It can then be shown that for all matrices $B\\in \\mathbb{C}^{m \\times n}$ with rank $k\\leq r$, then\n",
    "\n",
    "$$\n",
    "    ||A - A_{k}||_{2} \\leq || A-B||_{2}  \n",
    "$$\n",
    "\n",
    "i.e. that $A_k$ is the best rank-k approximation to $A$ in the 2-norm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A similar result can be shown for the Frobenious norm\n",
    "\n",
    "- For any $k$ with $0 \\leq k \\leq r$, the matrix $A_{k}$ also satisfies\n",
    "\n",
    "$$\n",
    "    ||A - A_{k}||_{F} = \\text{inf}_{B \\in \\mathbb{C}^{m \\times n}} \\text{rank}(B)\\leq v ||A-B||_{F} = \\sqrt{\\sigma_{v+1}^{2} + ... + \\sigma_{r}^{2}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Example:   Image compression\n",
    "\n",
    "How does this work in practice?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### the original Matrix (From Durer's Melancholia)\n",
    "\n",
    "<img align=center src=\"./images/Durer_Melancholia_I.jpg\" width=600>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib import image\n",
    "data = image.imread('images/melancholia-magic-square.png')\n",
    "m,n = data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('{}x{} pixel image'.format(m,n))\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.imshow(data,cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "u, s, vt = numpy.linalg.svd(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(s,'bo-')\n",
    "axes.set_ylabel('$\\sigma$', fontsize=16)\n",
    "axes.grid()\n",
    "axes.set_title('Spectrum of Singular Values')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### First 4 Modes\n",
    "\n",
    "where mode $i$ is the rank-1 matrix $\\sigma_i\\mathbf{u}_i\\mathbf{v}_i^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 4)\n",
    "fig.set_figheight(fig.get_figheight() * 1)\n",
    "for i in range(4):\n",
    "    mode = s[i]*numpy.outer(u[:,i], vt[i,:])\n",
    "\n",
    "    axes = fig.add_subplot(1, 4, i+1)\n",
    "    mappable = axes.imshow(mode, cmap='gray')\n",
    "    axes.set_title('Mode = {}, $\\sigma$={:3.3f}'.format(i+1,s[i]),fontsize=18)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Compressed Reconstructions\n",
    "\n",
    "We can now view compressed reconstructions that sum the first $k$ modes\n",
    "$$\n",
    "    A_k = \\sum_{i=1}^k \\sigma_i\\mathbf{u}_i\\mathbf{v}_i^T = U_k\\Sigma_kV_k^T\n",
    "$$\n",
    "\n",
    "the amount of memory required to store $A_k$ is $O(k(m +n + 1))$ whereas the storage of $A$ is $O(mn)$ so the relative compression is\n",
    "\n",
    "$$\n",
    "    \\frac{k(m+n + 1)}{mn}\n",
    "$$\n",
    "\n",
    "The question is how many modes are required to adequately reconstruct the original image?\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "fig.set_figheight(fig.get_figheight() *2)\n",
    "for i,k in enumerate([1, 10, 20]):\n",
    "    Ak = u[:,:k].dot(numpy.diag(s[:k]).dot(vt[:k,:]))\n",
    "    storage = 100.*k*(m+n + 1)/(m*n)\n",
    "    axes = fig.add_subplot(1, 3, i+1)\n",
    "    mappable = axes.imshow(Ak, vmin=0.0, vmax=1.0, cmap='gray')\n",
    "    axes.set_title('$A_{{{}}}$, $storage={:2.2f}\\%$'.format(k,storage),fontsize=16)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "fig.set_figheight(fig.get_figheight() *2)\n",
    "for i,k in enumerate([50, 100, 200]):\n",
    "    Ak = u[:,:k].dot(numpy.diag(s[:k]).dot(vt[:k,:]))\n",
    "    storage = 100.*k*(m+n + 1)/(m*n)\n",
    "    axes = fig.add_subplot(1, 3, i+1)\n",
    "    mappable = axes.imshow(Ak, vmin=0.0, vmax=1.0, cmap='gray')\n",
    "    axes.set_title('$A_{{{}}}$, $storage={:2.2f}\\%$'.format(k,storage),fontsize=16)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Other Applications\n",
    "\n",
    "* Total Least-squares\n",
    "* PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Total Least Squares\n",
    "\n",
    "**GOAL:** Demonstrate the use of the SVD to calculate total least squares regression and compare it to the classical least squares problem that assumes only errors in y."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Random data: \n",
    "\n",
    "We start by constructing a random data set that approximates a straight line but has random errors in both x and y coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# npoints uniformly randomly distributed points in the interval [0,3]\n",
    "npnts =100\n",
    "x = numpy.random.uniform(0.,3.,npnts)\n",
    "# set y = mx + b plus random noise of size err\n",
    "slope = 2.\n",
    "intercept = 1.\n",
    "err = .5\n",
    "\n",
    "y = slope*x + intercept \n",
    "y += numpy.random.normal(loc=y,scale=err)\n",
    "\n",
    "# add some random noise to x variable as well\n",
    "x += numpy.random.normal(loc=x,scale=err)\n",
    "\n",
    "# And plot out the data\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.scatter(x,y)\n",
    "axes.set_xlabel('x',fontsize=16)\n",
    "axes.set_ylabel('y', fontsize=16)\n",
    "axes.set_title('Data', fontsize=18)\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Classical Least Squares:  \n",
    "\n",
    "We first calculate the best fit straight line assuming all the error is in the $y$ variable using the a $QR$ decomposition of the Vandermonde matrix $\\begin{bmatrix} 1 & \\mathbf{x}\\\\\\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Vandermonde matrix\n",
    "A = numpy.array([ numpy.ones(x.shape), x]).T\n",
    "\n",
    "# solve  Ac = y using the QR decomposition via scipy\n",
    "c_ls,res,rank,s = numpy.linalg.lstsq(A,y, rcond=None)\n",
    "print('Best fit Linear Least Squares:')\n",
    "print('    slope={}'.format(c_ls[1]))\n",
    "print('    intercept={}'.format(c_ls[0]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# dummy variables \n",
    "t_ls = numpy.linspace(0,x.max())\n",
    "\n",
    "# And plot out the data\n",
    "fig  = plt.figure(figsize=(10,8))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.scatter(x,y)\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('y')\n",
    "# plot the least squares solution\n",
    "axes.plot(t_ls,c_ls[0]+t_ls*c_ls[1],'r-',label='Least Squares')\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Total Least Squares:  \n",
    "\n",
    "Suppose we wanted to find the line defined by its unit normal vector $\\mathbf{u}$ that instead minimized the sum of **orthogonal** distance between the data and the line, i.e."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# make up a line through the center of the data\n",
    "X = numpy.array([ x , y]).T\n",
    "X_mean = numpy.mean(X,0)\n",
    "v = numpy.array([1,2])\n",
    "v = v/numpy.linalg.norm(v, ord=2)\n",
    "u = numpy.array([-v[1],v[0]])\n",
    "#print(u.dot(v))\n",
    "t = numpy.linspace(0,x.max())\n",
    "t =  2*(t - numpy.mean(t_ls)) \n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.scatter(x,y)\n",
    "l = X_mean + numpy.outer(t,v)\n",
    "lu = X_mean + u\n",
    "axes.plot(l[:,0],l[:,1],'g')\n",
    "axes.plot([ X_mean[0], X_mean[0] + u[0]], [X_mean[1], X_mean[1] + u[1]], 'r')\n",
    "axes.text( X_mean[0] + 1.3*u[0],  X_mean[1] + 1.3*u[1], '$\\mathbf{u}$',fontsize=14)\n",
    "axes.set_aspect('equal')\n",
    "axes.set_xlabel('x',fontsize=16)\n",
    "axes.set_ylabel('y', fontsize=16)\n",
    "axes.set_title('Data', fontsize=18)\n",
    "\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This is equivalent to finding $\\mathbf{u}$ that minimizes\n",
    "\n",
    "$$\n",
    "    ||M\\mathbf{u}||_2\n",
    "$$ \n",
    "where \n",
    "\n",
    "$$\n",
    "    M = X - \\bar{X} =                                        \n",
    "\\begin{bmatrix} (x_1 - \\bar{x}) & (y_1 - \\bar{y}) \\\\\n",
    "                                        (x_2 - \\bar{x}) & (y_2 - \\bar{y}) \\\\\n",
    "                                        (x_3 - \\bar{x}) & (y_3 - \\bar{y}) \\\\\n",
    "                                        \\vdots & \\vdots \\\\\n",
    "                                        (x_n - \\bar{x}) & (y_n - \\bar{y}) \\\\\n",
    "                                       \\end{bmatrix}\n",
    "$$\n",
    "is a matrix of the de-meaned data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "but \n",
    "$$\n",
    "    ||M\\mathbf{u}||^2_2 = \\mathbf{u}^TM^TM\\mathbf{u}\n",
    "$$ \n",
    "\n",
    "or substituting in the SVD of $M=U\\Sigma V^T$ the problem becomes find $\\mathbf{u}$ that **minimizes**\n",
    "\n",
    "$$ \n",
    "    (\\mathbf{u}^TV) \\Sigma^T\\Sigma (V^T\\mathbf{u}) \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let $\\mathbf{y} = V^T\\mathbf{u} \\in \\mathbb{R}^2$ where $||\\mathbf{y}||_2 =1$\n",
    "\n",
    "then we seek $\\mathbf{y}$ that minimizes\n",
    "\n",
    "$$ \n",
    "    \\mathbf{y}^T\\Sigma^T\\Sigma \\mathbf{y} = \\sigma_1^2 y_1^2 + \\sigma_2^2 y_2^2 \\geq \\sigma_2^2(y_1^2 + y_2^2) \\geq \\sigma_2^2\n",
    "$$\n",
    "\n",
    "because $$||\\mathbf{y}||_2 =  y_1^2 + y_2^2 = 1$$\n",
    "\n",
    "hopefully it's clear that the vector \n",
    "\n",
    "$$\n",
    "    \\mathbf{y} = \\begin{bmatrix} 0 \\\\ 1 \\\\\\end{bmatrix}\n",
    "$$ \n",
    "provides the minimum value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "But  $\\mathbf{u} = V\\mathbf{y}$ \n",
    "\n",
    "So, the best fit line is the one that is orthogonal to $\\mathbf{v}_2$ which is the singular vector corresponding to the smallest singular value.\n",
    "\n",
    "The line we seek is orthogonal to that, so must lie in the direction of $\\mathbf{v}_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Total Least-Squares (a recipe)\n",
    "\n",
    "* form a matrix $M$ whose rows are the de-meaned data\n",
    "* find the SVD of $M = U\\Sigma V^T$, \n",
    "* set $\\mathbf{u} = \\mathbf{v}_2$ \n",
    "* and in 2-D,  the best fit line is parallel to $\\mathbf{v}_1$ and goes through the mean of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Prepare the data matrix\n",
    "X = numpy.array([ x , y]).T\n",
    "print('Shape of data Matrix: {}'.format(X.shape))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# and remove the mean\n",
    "X_mean = numpy.mean(X,0)\n",
    "print('Mean of data matrix={}'.format(X_mean))\n",
    "M = X - X_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# now calculate the SVD of the de-meaned data matrix\n",
    "U,S,VT = numpy.linalg.svd(M,full_matrices=False)\n",
    "V = VT.T\n",
    "print('Singular values = {}'.format(S))\n",
    "print('First Right singular vector V_1= {}'.format(V[:,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now plot and compare the two solutions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# dummy variables \n",
    "t_ls = numpy.linspace(0,x.max())\n",
    "t_svd =  2*(t_ls - numpy.mean(t_ls))\n",
    "\n",
    "# make figure\n",
    "plt.figure(figsize=(10,8))\n",
    "# plot data\n",
    "plt.scatter(x,y)\n",
    "# plot the least squares solution\n",
    "plt.plot(t_ls,c_ls[0]+t_ls*c_ls[1],'r-',label='Least Squares')\n",
    "\n",
    "# plot the total least Squares solution\n",
    "# plot the mean\n",
    "plt.plot(X_mean[0],X_mean[1],'go', markersize=12, label='X_mean')\n",
    "# calculate a line through the mean with the first principal component as a basis\n",
    "L_tls = X_mean + numpy.outer(t_svd,V[:,0])\n",
    "plt.plot(L_tls[:,0],L_tls[:,1],'c-',label='Total Least Squares')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title('Comparison Least Squares vs Total Least Squares')\n",
    "plt.legend(loc='best')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Notes\n",
    "\n",
    "This technique can be extended into higher dimensions.  For example, we could find a best fit plane through a 3-D cloud of data.\n",
    "\n",
    "For this problem again, we would find the SVD of the demeaned data matrix $M$, and the best fit plane would be normal to the last singular vector $\\mathbf{v}_3$ or  $\\mathrm{span}\\langle \\mathbf{v}_1, \\mathbf{v}_2\\rangle$ \n",
    "\n",
    "\n",
    "In general, for a demeaned data matrix $M$,  columns of the matrix $V$ form an orthonormal basis for the row space of $M$ and describe the axes of the best fit Ellipsoid describing the data. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Principal Component/EOF analysis\n",
    "\n",
    "**GOAL:** Demonstrate the use of the SVD to calculate principal components or \"Empirical Orthogonal Functions\" in a geophysical data set.  This example is modified from a paper by Chris Small (LDEO)\n",
    "\n",
    "Small, C., 1994. A global analysis of mid-ocean ridge axial topography. Geophys J Int 116, 64–84. [doi:10.1111/j.1365-246X.1994.tb02128.x](https://academic.oup.com/gji/article/116/1/64/638843/A-global-analysis-of-mid-ocean-ridge-axial)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## The Data\n",
    "\n",
    "Here we will consider a set of topography profiles taken across the global mid-ocean ridge system where the Earth's tectonic plates are spreading apart.\n",
    "<table>\n",
    "     <tr align=center><td><img align=center src=\"./images/World_OceanFloor_topo_green_brown_1440x720.jpg\"><td>\n",
    "</table>\n",
    "\n",
    "The data consists of 156 profiles from a range of spreading rates. Each profile contains 79 samples so is in effect a vector in $R^{79}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# read the data from the csv file\n",
    "data = pd.read_csv('data/m80_data.csv',header=None).values\n",
    "data_mean =  numpy.mean(data,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "# and plot out a few profiles and the mean depth.\n",
    "plt.figure(figsize=(10,8))\n",
    "rows = [ 9,59,99]\n",
    "labels = [ 'slow','medium','fast']\n",
    "for i,row in enumerate(rows):\n",
    "    plt.plot(data[row,:],label=labels[i])\n",
    "plt.plot(data_mean,'k--',label='mean')\n",
    "plt.xlabel('Distance across axis (km)')\n",
    "plt.ylabel('Relative Elevation (m)')\n",
    "plt.legend(loc='best')\n",
    "plt.title('Example cross-axis topography of mid-ocean ridges')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### EOF analysis\n",
    "\n",
    "While each profile lives in an 80 dimensional space, we would like to see if we can classify the variability in fewer components.  To begin we form a de-meaned data matrix $X$ where each row is a profile.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "X = data - data_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,8))\n",
    "plt.imshow(X, vmin=-1500, vmax=2000)\n",
    "plt.xlabel('Distance across axis (Km)', fontsize=16)\n",
    "plt.ylabel('Relative Spreading Rate', fontsize=16)\n",
    "plt.title('De-meaned Data', fontsize=18)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('Relative Depth (m)', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Applying the SVD\n",
    "\n",
    "We now use the SVD to factor the data matrix as $X = U\\Sigma V^T$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# now calculate the SVD of the de-meaned data matrix\n",
    "U,S,Vt = numpy.linalg.svd(X,full_matrices=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "And begin by looking at the spectrum of singular values $\\Sigma$.  Defining the variance as $\\Sigma^2$ then we can also calculate the cumulative contribution to the total variance as \n",
    "\n",
    "$$\n",
    "    g_k = \\frac{\\sum_{i=0}^k \\sigma_i^2}{\\sum_{i=0}^n \\sigma_i^2}\n",
    "$$\n",
    "\n",
    "Plotting both $\\Sigma$ and $g$ shows that $\\sim$ 80% of the total variance can be explained by the first 4-5 Components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# plot the singular values\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.semilogy(S,'bo')\n",
    "axes.grid()\n",
    "axes.set_title('Singular Values', fontsize=18)\n",
    "\n",
    "# and cumulative percent of variance\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "g = numpy.cumsum(S*S)/numpy.sum(S*S)\n",
    "axes.plot(g,'bx-')\n",
    "\n",
    "axes.set_title('% cumulative percent variance explained', fontsize=18)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Plotting the first 3 Singular Vectors in $V$,  shows them to reflect some commonly occuring patterns in the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "num_EOFs=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "for row in range(num_EOFs):\n",
    "    axes.plot(Vt[row,:],label='EOF{}'.format(row+1))\n",
    "\n",
    "axes.plot([0, X.shape[1]], [0., 0.], 'k--')\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('Distance (km)')\n",
    "axes.set_title('First {} EOFs '.format(num_EOFs))\n",
    "axes.legend(loc='best')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For example,  the first EOF pattern is primarily a symmetric pattern with an axial high surrounded by two off axis troughs (or an axial low with two flanking highs,  the EOF's are just unit vector bases for the row-space and can be added with any positive or negative coefficient).  The Second EOF is broader and all of one sign while the third EOF encodes assymetry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Reconstruction\n",
    "\n",
    "Using the SVD we can also decompose each profile into a weighted linear combination of EOF's i.e.\n",
    "\n",
    "$$\n",
    "    X = U\\Sigma V^T = C V^T\n",
    "$$\n",
    "\n",
    "where $C = U\\Sigma$ is a matrix of coefficients  that describes the how each data row is decomposed into the relevant basis vectors.  We can then produce a k-rank truncated representation of the data by\n",
    "\n",
    "$$ \n",
    "    X_k = C_k V_k^T\n",
    "$$\n",
    "\n",
    "where $C_k$ is the first $k$ columns of $C$ and $V_k$ is the first $k$ EOF's.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here we show the original data and the reconstructed data using the first 5 EOF's"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# recontruct the data using the first 5 EOF's\n",
    "k=5\n",
    "Ck = numpy.dot(U[:,:k],numpy.diag(S[:k]))\n",
    "Vtk = Vt[:k,:]\n",
    "data_k = data_mean + numpy.dot(Ck,Vtk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 8))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "im = axes.imshow(data_k, vmin=-1500, vmax=2000)\n",
    "fig.colorbar(im)\n",
    "axes.set_title('reconstructed data')\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "im = axes.imshow(data, vmin=-2000, vmax=2000)\n",
    "axes.set_title('Original data')\n",
    "fig.colorbar(im)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(Ck.T)\n",
    "plt.figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "And we can consider a few reconstructed profiles compared with the original data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# show the original 3 profiles and their recontructed values using the first k EOF's\n",
    "fig = plt.figure(figsize=(24,8))\n",
    "\n",
    "for i,row in enumerate(rows):\n",
    "    axes = fig.add_subplot(1, 3, i+1)\n",
    "    h = axes.plot(data_k[row,:])\n",
    "    h += axes.plot(data[row,:])\n",
    "    axes.grid()\n",
    "    Cstring = [ '{:3.0f},  '.format(Ck[row,i]) for i in range(k) ]\n",
    "    axes.set_title('Reconstruction profile {}:\\n C_{}='.format(row,k)+''.join(Cstring))\n",
    "    axes.legend(['k={}'.format(k), 'original data'], loc='best')\n",
    "\n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## projection of data onto a subspace\n",
    "\n",
    "We can also use the Principal Components to look at the projection of the data onto a lower dimensional space as the coefficients $C$, are simply the coordinates of our data along each principal component.  For example we can view the data in the 2-Dimensional space defined by the first 2 EOF's by simply plotting C_1 against C_2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# plot the data in the plane defined by the first two principal components\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.scatter(Ck[:,0],Ck[:,1])\n",
    "axes.set_xlabel('$V_1$')\n",
    "axes.set_ylabel('$V_2$')\n",
    "axes.grid()\n",
    "axes.set_title('Projection onto the first two principal components')\n",
    "\n",
    "# Or consider the degree of assymetry (EOF 3) as a function of spreading rate\n",
    "axes = fig.add_subplot(1,2,2)\n",
    "axes.plot(Ck[:,2],'o')\n",
    "axes.set_xlabel('Relative Spreading rate')\n",
    "axes.set_ylabel('$C_3$')\n",
    "axes.grid()\n",
    "axes.set_title('Degree of assymetry')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}