10.1_SVD.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",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular Value Decomposition\n",
    "\n",
    "Definition: Let $A \\in \\mathbb R^{m \\times n}$, then $A$ can be factored as\n",
    "$$\n",
    "    A = U\\Sigma V^{T}\n",
    "$$\n",
    "\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 non negative arranged in descending order. ($\\sigma_{1} \\geq \\sigma_{2} \\geq \\sigma_{3} \\geq ... \\sigma_{r} \\geq 0$).\n",
    "\n",
    "The SVD has a number of applications mostly related to reducing the dimensionality of a matrix."
   ]
  },
  {
   "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": "skip"
    }
   },
   "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)\n",
    "print(numpy.dot(U, numpy.dot(numpy.diag(sigma), V_T)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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": [
    "### 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": "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) = [u_1, ... , u_r]$ and $\\text{null}(a) = [v_{r+1}, ... , v_n]$.\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 n}$ 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 = \\sum_{j=1}^{k} \\sigma_{j}u_{j}v_{j}^{T}\n",
    "$$\n",
    "Let $k = min(m,n)$, then\n",
    "\n",
    "$$\n",
    "    ||A - A_{v}||_{2} = \\text{inf}_{B \\in \\mathbb{C}^{m \\times n}} \\text{rank}(B)\\leq k|| A-B||_{2} = \\sigma_{k+1}\n",
    "$$\n",
    "\n",
    "- For any $k$ with $0 \\leq k \\leq r$, the matrix $A_{k}$ also satisfies\n",
    "$$\n",
    "    ||A - A_{v}||_{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:   \"Hello World\"\n",
    "\n",
    "How does this work in practice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "data = numpy.zeros((15,40))\n",
    "\n",
    "#H\n",
    "data[2:10,2:4] = 1\n",
    "data[5:7,4:6] = 1\n",
    "data[2:10,6:8] = 1\n",
    "\n",
    "#E\n",
    "data[3:11,10:12] = 1\n",
    "data[3:5,12:16] = 1\n",
    "data[6:8, 12:16] = 1\n",
    "data[9:11, 12:16] = 1\n",
    "\n",
    "#L\n",
    "data[4:12,18:20] = 1\n",
    "data[10:12,20:24] = 1\n",
    "\n",
    "#L\n",
    "data[5:13,26:28] = 1\n",
    "data[11:13,28:32] = 1\n",
    "\n",
    "#0\n",
    "data[6:14,34:36] = 1\n",
    "data[6:8, 36:38] = 1\n",
    "data[12:14, 36:38] = 1\n",
    "data[6:14,38:40] = 1\n",
    "\n",
    "plt.imshow(data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "u, diag, vt = numpy.linalg.svd(data, full_matrices=True)\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "fig.set_figheight(fig.get_figheight() * 4)\n",
    "for i in range(1, 16):\n",
    "    diag_matrix = numpy.concatenate((numpy.zeros((len(diag[:i]) -1),), diag[i-1: i], numpy.zeros((40-i),)))\n",
    "    reconstruct = numpy.dot(numpy.dot(u, numpy.diag(diag_matrix)[:15,]), vt)\n",
    "    \n",
    "    axes = fig.add_subplot(5, 3, i)\n",
    "    mappable = axes.imshow(reconstruct, vmin=0.0, vmax=1.0)\n",
    "    axes.set_title('Component = %s' % i)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "u, diag, vt = numpy.linalg.svd(data, full_matrices=True)\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 3)\n",
    "fig.set_figheight(fig.get_figheight() * 4)\n",
    "for i in range(1, 16):\n",
    "    diag_matrix = numpy.concatenate((diag[:i], numpy.zeros((40-i),)))\n",
    "    reconstruct = numpy.dot(numpy.dot(u, numpy.diag(diag_matrix)[:15,]), vt)\n",
    "\n",
    "    axes = fig.add_subplot(5, 3, i)\n",
    "    mappable = axes.imshow(reconstruct, vmin=0.0, vmax=1.0)\n",
    "    axes.set_title('Component = %s' % i)\n",
    "    \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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}