notes/SVD-eigenproblem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The SVD as an eigenproblem\n",
    "\n",
    "Notice that if $A = U\\Sigma V^H$, then\n",
    "\n",
    "$$\n",
    "A^H A = V \\Sigma U^H U \\Sigma V^H = V \\Sigma^2 V^H\n",
    "$$\n",
    "\n",
    "That is, to multiply $A^H A x$, you (1) compute $V^H x$ (the $V$ components of $x$), then (2) multiply each component by $\\sigma^2$, and finally (3) multiply the coefficients by $V$ and add up.  It follows that:\n",
    "\n",
    "* The singular values $\\sigma^2$ are the **nonzero eigenvalues** of $A^H A$ and the corresponding **eigenvectors are the right singular vectors** $V$.\n",
    "\n",
    "Similarly,\n",
    "\n",
    "$$\n",
    "A A^H = U \\Sigma V^H V \\Sigma U^H = U \\Sigma^2 U^H\n",
    "$$\n",
    "\n",
    "so\n",
    "\n",
    "* The singular values $\\sigma^2$ are the **nonzero eigenvalues** of $A A^H$ and the corresponding **eigenvectors are the left singular vectors** $U$.\n",
    "\n",
    "We can easily check this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 Array{Float64,2}:\n",
       "  0.202935  -0.810741  0.379812 \n",
       "  0.317852   1.17222   0.0789665\n",
       " -1.58283   -0.524304  0.949145 \n",
       "  0.122448   1.57466   0.693527 \n",
       " -0.496476   1.13621   0.511883 "
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = randn(5,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in this case, $A$ is a $5 \\times 3$ matrix of rank 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 6.31957 \n",
       " 4.01707 \n",
       " 0.443596"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U, σ, V = svd(A)\n",
    "\n",
    "σ.^2    # the σ² values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$A^H A$ is a $3 \\times 3$ matrix of rank 3 with three nonzero eigenvalues that equal the singular values squared:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.443596\n",
       " 4.01707 \n",
       " 6.31957 "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(A'*A) # AᴴA has the same eigenvals!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$AA^H$ is a $5 \\times 5$ matrix of rank 3 (recall that the ranks of $A$, $AA^H$, and $A^H A$ are all equal!).   It has 3 nonzero eigenvalues that equal the $\\sigma^2$ values, and 2 zero eigenvalues corresponding to the **two-dimensional** nullspace\n",
    "$$\n",
    "N(AA^H) = N(A^H) = C(A)^\\perp = C(U)^\\perp\n",
    "$$\n",
    "That is, the zero eigenvectors are those perpendicuar to the left singular vectors $U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " -1.71137e-16\n",
       "  8.87578e-16\n",
       "  0.443596   \n",
       "  4.01707    \n",
       "  6.31957    "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(A*A') # the same *nonzero* eigenvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the eigenvectors, e.g. the eigenvectors of $A^H A$ should match V:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  0.564298  -0.817673    0.113922\n",
       " -0.203239  -0.00384465  0.979122\n",
       "  0.800163   0.57567     0.168353"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvecs(A'A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -0.113922  -0.817673    -0.564298\n",
       " -0.979122  -0.00384465   0.203239\n",
       " -0.168353   0.57567     -0.800163"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes, they match, up to an overall sign flip.\n",
    "\n",
    "(Note that the columns are in reverse order, because `svdvals` are by default sorted in *descending* order in Julia, whereas `eigvals` of Hermitian matrices are sorted in *ascending* order.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remarks\n",
    "\n",
    "* This, in principle, finally gives us a way to compute the SVD of a matrix: just find the eigenvectors and eigenvalues of $A^H A$.  (Note that $AV = U\\Sigma$, so that once you have $V$ and $\\Sigma$ you can get $U$.)\n",
    "\n",
    "* In practice, computers use a different way to compute the SVD.  (The most famous practical method is called \"Golub-Kahan bidiagonalization.\")  In 18.06, we are content to let the `svd` function be a \"black box\", much like `eig`.\n",
    "\n",
    "* The fact that the singular values/vectors are related to eigenvalues of $A^H A$ and $A A^H$ has lots of important applications.  Perhaps most famously, it means that the SVD diagonalizes the \"covariance matrix\" in statistics, which gives rise to the statistical method of [principal component analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.3",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}