LinearAlgebra_SVD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "remove_cell"
    ]
   },
   "source": [
    "# Singular Value Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The singular value decomposition is closely associated with the eigenvalue-eigenvector factorization $\\mathbf{Q} \\Lambda \\mathbf{Q}^\\intercal$ of a positive definite matrix. However, we encounter many matrices that are not positive definite and the decomposition is not true. Consider a rectangular matrix (non-square) matrix, the eigenvalues are not defined in this scenario. If we allow the $\\mathbf{Q}$ on the left and the $\\mathbf{Q}$ on the right to be any two orthogonal matrices $\\mathbf{U}$ and $\\mathbf{V}^\\intercal$, then any matrix will have a decomposition $\\mathbf{A=U}\\Sigma\\mathbf{V}^\\intercal$. \n",
    "\n",
    "$\\Sigma$ has eigenvalues from $\\mathbf{A}^\\intercal\\mathbf{A}$, not from $\\mathbf{A}$. These positive entries will be $\\sigma_1,\\ldots,\\sigma_r$. They are the singular values of $\\mathbf{A}$. They fill the first $r$ places on the main diagonal of $\\Sigma$ when $\\mathbf{A}$ has rank $r$. The rest of $\\Sigma$ is zero. The columns of $\\mathbf{U}$ are the eigenvectors of $\\mathbf{AA}^\\intercal$, and the columns of $\\mathbf{V}$ are the eigenvectors of $\\mathbf{A}^\\intercal\\mathbf{A}$. The $r$ singular values on the diagonal of $\\Sigma$ are the square roots of the nonzero eigenvalues of both $\\mathbf{AA}^\\intercal$ and $\\mathbf{A}^\\intercal\\mathbf{A}$.\n",
    "\n",
    "The *SVD* is really good for numerically stable computations, because $\\mathbf{U}$ and $\\mathbf{V}$ are orthogonal matrices. They never change the length of a vector. $\\Sigma$ could multiply or divide a vector by a large $\\sigma$, but at least now we have an idea of exactly what is large and what is small. The ratio of $\\sigma_{\\mbox{max}} / \\sigma_{\\mbox{min}}$ is the *condition number* of a nonsingular matrix. For us as numerical analysts, this might be the most important use of the *SVD*. Other common uses of the *SVD* are listed below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " 1. Image processing\n",
    " 2. Effective rank\n",
    " 3. Polar decomposition\n",
    " 4. Least squares\n",
    " 5. Pseudoinverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and even... the algorithm that computes eHarmony.com's 29 Dimensions$\\circledR$ of compatability (U.S. Patent 6,735,568)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## SVD Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We already have elementary forms of all the tools required to calculate a symmetric matrix's *SVD*. We can use the $\\mathbf{QR}$ algorithm to find the eigenvalues of the $\\mathbf{AA}^\\intercal$ matrix. Then use the inverse Power Method to find the eigenvectors that correspond to the singular values for both the $\\mathbf{AA}^\\intercal$ and $\\mathbf{A}^\\intercal\\mathbf{A}$. Finally we can use the Gram-Schmidt method to find an orthogonal basis of these eigenvectors to define the $\\mathbf{U}$ and $\\mathbf{V}$ matrices.  This is the process that is performed in [`np.linalg.svd`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove_cell"
    ]
   },
   "source": [
    " 1. Use $\\mathbf{QR}$ method to find eigenvalues of $\\mathbf{A}^\\intercal\\mathbf{A}$\n",
    " 2. Use inverse Power Method to find the eigenvectors of $\\mathbf{A}^\\intercal\\mathbf{A}$ and $\\mathbf{A}\\mathbf{A}^\\intercal$\n",
    " 3. Use Gram-Schmidt to orthonormalize these eigenvectors to define $\\mathbf{U}$ and $\\mathbf{V}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More on Condition Number and Preconditioning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can use the condition number of a matrix to determine the accuracy of our solution.  For the matrix equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\mathbf{A}\\vec{x} = \\vec{b}$$\n",
    "\n",
    "The relative error is $\\Vert \\vec{x}\\Vert < \\sigma_{\\mbox{max}} / \\sigma_{\\mbox{min}} \\epsilon_m$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The goal of matrix preconditioner (we used them for iterative techniques in solving linear systems) is to reduce the condition number of the preconditioned matrix with respect to the original matrix.  These are most useful in iterative techniques because the rate of convergence of most iterative solvers increases as the condition number of the preconditioned matrix decreases.  Preconditioned iterative solvers usually outperform direct methods, especially for sparse matrices.  This is why we typically use iterative solvers for the Finite Element Methods, where the equations of motion generally result in very sparse matrices. \n",
    "\n",
    "Instead of solving the original equation above we would solve something like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\\mathbf{P}^{-1} \\left(\\mathbf{A}\\vec{x} - \\vec{b}\\right) = \\vec{0}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "where, $\\mathbf{P}$ is the preconditioner.  Usually there is some trade off for the computational expense of computing $\\mathbf{P}^{-1}$ that guides our selection of $\\mathbf{P}$.  The least computational expensive preconditioner would be $\\mathbf{P=I}$, however this does nothing to the original system.  The best preconditioner would be $\\mathbf{P}^{-1}=\\mathbf{A}^{-1}$, which has a condition number of 1, but this is as computationally expensive as solving the original system.  Therefore, our choice of a preconditioner is usually somewhere in between.  We strive for a minimal number of iterations for convergence and a minimal expense in computing $\\mathbf{P}^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Survey of Software for use in Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If you ever find the need to write software that uses a lot of matrix manipulations and/or linear systems solving techniques you can use the subroutine library LAPACK and it's accompanying lower-level operations package called BLAS. These package libraries are written in C and/or FORTRAN and can be called from C/C++ or even Python (although they are used internally by NumPy and SciPy) programs. They run well on shared memory and parallel processing machines. The following table shows equivalent or nearly equivalent commands between LAPACK, Python/NumPy/SciPy. This list is not complete, but shows a few of the operations we have covered while studying numerical techniques for linear systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| LAPACK   | Python/NumPy/SciPy |\n",
    "|:--------:|:--------:|\n",
    "|  `SGETRF`  | `scipy.lingalg.lu`\n",
    "|  `SGERTS`  | `np.linalg.solve`\n",
    "|  `SPOTRF`  | `scipy.linalg.cholesky`\n",
    "|  `SPOTRI`  | `np.linalg.inv`\n",
    "|  `SSTEQR`  | `np.linalg.eigvals`\n",
    "|  `STREVC`  |\t`np.linalg.eigvecs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide_input": true,
    "init_cell": true,
    "javascript_last_cell": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove_cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "function hideElements(elements, start) {\n",
       "for(var i = 0, length = elements.length; i < length;i++) {\n",
       "    if(i >= start) {\n",
       "        elements[i].style.display = \"none\";\n",
       "    }\n",
       "}\n",
       "}\n",
       "var prompt_elements = document.getElementsByClassName(\"prompt\");\n",
       "hideElements(prompt_elements, 0)\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "function hideElements(elements, start) {\n",
    "for(var i = 0, length = elements.length; i < length;i++) {\n",
    "    if(i >= start) {\n",
    "        elements[i].style.display = \"none\";\n",
    "    }\n",
    "}\n",
    "}\n",
    "var prompt_elements = document.getElementsByClassName(\"prompt\");\n",
    "hideElements(prompt_elements, 0)"
   ]
  }
 ],
 "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.8.1"
  },
  "rise": {
   "auto_select": null,
   "autolaunch": false,
   "controls": true,
   "enable_chalkboard": false,
   "footer": "",
   "progress": true,
   "scroll": true,
   "slideNumber": false,
   "theme": "simple"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}