15_guest/15_01_svd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SVD Math\n",
    "\n",
    "Let's explore the usage of SVD decomposition in numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sb\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Singular Value Decomposition\n",
    "\n",
    "#### SVD\n",
    "\n",
    "If A is a $N \\times n$ matrix, then the _Singular Value Decomposition_ (SVD) of $A$ is given by \n",
    "\n",
    "$$A = U \\Sigma V^{T}$$\n",
    "\n",
    "where\n",
    "\n",
    "- $U$ and $V$ are orthogonal matrices, meaning that $UU^T = \\textbf{1}$ and $VV^T = \\textbf{1}$ are identity matrices, and\n",
    "- $\\Sigma$ is a diagonal matrix, meaning that the values lie *only* on the diagonal (This does not mean that $\\Sigma$ is also a square matrix!)\n",
    "\n",
    "\n",
    "- U shape: $N \\times N$\n",
    "- $\\Sigma$ shape: $N \\times n$\n",
    "- VT shape: $n \\times n$\n",
    "\n",
    "In particular, $U$ and $V$ are square matrices consisting of the eigenvectors of $AA^T$ and $A^TA$, respectively. Note that $\\Sigma$ is _not_ necessarily a square matrix (since generally $N \\neq n$). Its values are non-negative real numbers, and are the _singular values_ of $A$, which are defined as the square root of the eigenvalues of $AA^T$ (or $A^TA$, which are the same).\n",
    "\n",
    "Let's see that in action. Let's start with some arbitraty matrix $M$ and compute its SVD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.3588728 ,  0.39529603,  0.69319986],\n",
       "       [ 0.98972943,  0.38322409,  0.84533543],\n",
       "       [ 0.33124072,  0.63393481,  0.22036526],\n",
       "       [ 0.7278726 ,  0.12967207,  0.09276765],\n",
       "       [ 0.21235531,  0.65255971,  0.52620262]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = np.random.random((5,3))\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.42959002 -0.23510539  0.44556072  0.34822001 -0.66362364]\n",
      " [-0.67927809  0.3796194   0.31560578 -0.47119554  0.26988542]\n",
      " [-0.33254454 -0.33050489 -0.6767384  -0.4432501  -0.3545915 ]\n",
      " [-0.29777215  0.60120123 -0.48517171  0.55828758 -0.05303043]\n",
      " [-0.39342033 -0.57440244 -0.09220735  0.3854373   0.5985133 ]]\n",
      "[ 1.94905419  0.69247374  0.43755181]\n",
      "[[-0.63461913 -0.48037942 -0.60538762]\n",
      " [ 0.71842635 -0.65540155 -0.23305018]\n",
      " [-0.28481948 -0.58282452  0.7610476 ]]\n"
     ]
    }
   ],
   "source": [
    "U, s, VT = np.linalg.svd(M)  # returns matrix, array (only diagonal)\n",
    "print U\n",
    "print s\n",
    "print VT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `numpy`'s `linalg.svd` returns `U, s, VT`, where\n",
    "- `s` is the actual diagonal of $\\Sigma$ (not the entire matrix) and\n",
    "- `VT` is the transpose of $V$.\n",
    "\n",
    "Note that $\\Sigma$ is a diagonal matrix, but not a square matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.94905419,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.69247374,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.43755181],\n",
       "       [ 0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = np.zeros(M.shape)\n",
    "S[:min(M.shape), :min(M.shape)] = np.diag(s)\n",
    "S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check if the SVD equation holds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.isclose(M, U.dot(S).dot(VT)).all() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good!\n",
    "\n",
    "Let's check if the values in $\\Sigma$ are indeed the square roots of the eigenvalues of $AA^T$ resp. $A^TA$.  \n",
    "\n",
    "We sort the eigenvalues, so we can compare them easily. We also round down to some decimal to clearly show the values that are nearly zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.7988122  0.4795199  0.1914516 -0.        -0.       ]\n",
      "[ 3.7988122  0.4795199  0.1914516]\n",
      "[ 3.7988122  0.4795199  0.1914516]\n"
     ]
    }
   ],
   "source": [
    "values1 = np.sort(np.linalg.eig(M.dot(M.T))[0])[::-1].round(7)  \n",
    "values2 = np.sort(np.linalg.eig(M.T.dot(M))[0])[::-1].round(7)\n",
    "squares = np.sort(np.square(s))[::-1].round(7)\n",
    "print values1\n",
    "print values2\n",
    "print squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's inspect $U$ check if its columns are indeed equal to the eigenvectors of $AA^T$.  \n",
    "\n",
    "We sort the eigenvectors according to their eigenvalue, so we can compare them easily again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.42959002 -0.23510539  0.44556072  0.34822001 -0.66362364]\n",
      " [-0.67927809  0.3796194   0.31560578 -0.47119554  0.26988542]\n",
      " [-0.33254454 -0.33050489 -0.6767384  -0.4432501  -0.3545915 ]\n",
      " [-0.29777215  0.60120123 -0.48517171  0.55828758 -0.05303043]\n",
      " [-0.39342033 -0.57440244 -0.09220735  0.3854373   0.5985133 ]]\n",
      "[[-0.42959002 -0.23510539 -0.44556072  0.36595791 -0.31059127]\n",
      " [-0.67927809  0.3796194  -0.31560578 -0.47828983 -0.07582197]\n",
      " [-0.33254454 -0.33050489  0.6767384  -0.43354423 -0.55192107]\n",
      " [-0.29777215  0.60120123  0.48517171  0.55951279  0.30051763]\n",
      " [-0.39342033 -0.57440244  0.09220735  0.36918625  0.70912367]]\n",
      "[[ 1.          1.         -1.          0.95153024  2.13664616]\n",
      " [ 1.          1.         -1.          0.98516738 -3.55946202]\n",
      " [ 1.          1.         -1.          1.02238727  0.64246777]\n",
      " [ 1.          1.         -1.          0.99781022 -0.17646363]\n",
      " [ 1.          1.         -1.          1.04401857  0.84401822]]\n"
     ]
    }
   ],
   "source": [
    "values, vectors = np.linalg.eig(M.dot(M.T))\n",
    "vectors = vectors[:, values.argsort()[::-1]]\n",
    "print U\n",
    "print vectors\n",
    "print U / vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that they're equal up to a possible sign (+/-), since the orientation of eigenvectors and eigenvalues is arbitrary. Also note that the last eigenvectors are just noise as their corresponding eigenvalues are near-zero.\n",
    "\n",
    "Let's inspect $V$ and $A^TA$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.63461913  0.71842635 -0.28481948]\n",
      " [-0.48037942 -0.65540155 -0.58282452]\n",
      " [-0.60538762 -0.23305018  0.7610476 ]]\n",
      "[[-0.63461913 -0.71842635 -0.28481948]\n",
      " [-0.48037942  0.65540155 -0.58282452]\n",
      " [-0.60538762  0.23305018  0.7610476 ]]\n",
      "[[ 1. -1.  1.]\n",
      " [ 1. -1.  1.]\n",
      " [ 1. -1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "V = VT.T\n",
    "values, vectors = np.linalg.eig(M.T.dot(M))\n",
    "vectors = vectors[:, values.argsort()[::-1]]\n",
    "print V\n",
    "print vectors\n",
    "print V / vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, that all seems to work!\n",
    "\n",
    "#### Truncated SVD\n",
    "\n",
    "How does this decomposition help us with data science?  As we saw in the example above already, if $\\Sigma$ is not a square matrix, a lot of the column vectors in $U$ are irrelevant since they correspond to zero-value eigenvalues. Moreover, in practice, you'd have many eigenvalues that are very small. Since the values in $\\Sigma$ are ordered from big to small, we could significantly reduce the dimensions of $\\Sigma$ by limiting this matrix to an actual square matrix of dimension $k$, only containing the first $k$ rows and columns of $\\Sigma$, and throwing away the smallest singular values.\n",
    "\n",
    "\\begin{align*}\n",
    "A = U \\Sigma V^{T} \\approx U' \\Sigma' V'^{T} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "We often ommit these accents and simply write $A = U \\Sigma V^{T}$ when we in fact mean the _truncated SVD_ decomposition. Since we truncate only the smallest singular values, this loss in accuracy is only small if we choose $k$ large enough. This is a small price for a huge gain in insight in our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.3588728   0.39529603  0.69319986]\n",
      " [ 0.98972943  0.38322409  0.84533543]\n",
      " [ 0.33124072  0.63393481  0.22036526]\n",
      " [ 0.7278726   0.12967207  0.09276765]\n",
      " [ 0.21235531  0.65255971  0.52620262]]\n",
      "[[ 0.41440004  0.50892111  0.54482914]\n",
      " [ 1.02906126  0.46370859  0.74023941]\n",
      " [ 0.24690337  0.46135574  0.44571762]\n",
      " [ 0.66740892  0.00594556  0.25432874]\n",
      " [ 0.20086413  0.62904537  0.55690746]]\n"
     ]
    }
   ],
   "source": [
    "n_samples, n_features = M.shape\n",
    "k = 2 \n",
    "M_svd = U[:n_samples, :k].dot(S[:k, :k]).dot(V[:n_features, :k].T)\n",
    "print M\n",
    "print M_svd "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to measure the quality of this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M is a 5x3-matrix, but we only choose the largest 2 singular values\n",
      "\n",
      "We kept 85% of our singular values, weighted by their size:\n",
      "- before: [ 1.94905419  0.69247374  0.43755181]\n",
      "- after:  [ 1.94905419  0.69247374]\n",
      "\n",
      "Relative differences with original matrix (0.10 on average):\n",
      "[[ 0.15  0.29 -0.21]\n",
      " [ 0.04  0.21 -0.12]\n",
      " [-0.25 -0.27  1.02]\n",
      " [-0.08 -0.95  1.74]\n",
      " [-0.05 -0.04  0.06]]\n"
     ]
    }
   ],
   "source": [
    "s_kept = s[:k]\n",
    "print \"M is a %dx%d-matrix, but we only choose the largest %d singular values\" % (n_samples, n_features, k)\n",
    "print \"\\nWe kept %d%% of our singular values, weighted by their size:\" % (100. * s_kept.sum() / s.sum())\n",
    "print \"- before:\", s\n",
    "print \"- after: \", s_kept\n",
    "diff = ((M_svd - M) / M).round(2)\n",
    "print \"\\nRelative differences with original matrix (%.2f on average):\" % diff.mean()\n",
    "print diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an individual data point, the difference can be substantial, but for the population as a whole, we are pretty close to our initial matrix. \n",
    "\n",
    "\n",
    "#### Interpretation\n",
    "\n",
    "Let's $A$ be our $N \\times n$ feature matrix with $N$ samples and $n$ features, and let's write the SVD as follows, trunacted to some dimension $k < n$:\n",
    "\n",
    "\\begin{align}\n",
    "A = U \\Sigma V^{T}\n",
    "\\end{align}\n",
    "\n",
    "Each of these matrices represents some grouping the data:\n",
    "- $U$ ($N \\times k$) represents how each datapoint relates to some latent concept,\n",
    "- $V$ ($n \\times k$) represents how each concept relates to each original feature, and\n",
    "- $\\Sigma$ ($k \\times k$) represents the strength of each of those concepts in the dataset.\n",
    "\n",
    "This reduces the number of dimensions from $n$ to $k$, since instead of modelling our data for each feature, we now model our data for each concept. Our intuition is that datapoints with similar relevant features will get similar concepts.\n",
    "\n",
    "A very concrete example of this is text classification, in which case this method is often called _Latent semantic Analysis_ (LSA).  As always with Natural Language Processing (NLP), we transform our $N$ documents into a $N \\times n$ feature matrix, where each column represents a certain word or n-gram (e.g., using `CountVectorizer`).  Please refer to the **notebook about LSA** in this folder for a clear demonstration about this.\n",
    "\n",
    "### Principal Component Analysis\n",
    "\n",
    "Principal Component Analysis (PCA) is directlty related to SVD by the following expression:\n",
    "\n",
    "\\begin{align*}\n",
    "A &= U \\Sigma V^{T} \\\\\n",
    "AA^T &= (U \\Sigma V^{T})(V \\Sigma U^{T}) \\\\\n",
    "&= U \\Sigma^2 U^{T}\n",
    "\\end{align*}\n",
    "\n",
    "Note that we used that $V$ is an orthogonal matrix, and hence $V^TV = \\textbf{1}$.  Then we have:\n",
    "- $AA^T$ is called the _covariance matrix_ (assuimg that all features $A$ were scaled with mean zero),\n",
    "- $U$ is a square matrix consisting of the eigenvectors of $AA^T$, and\n",
    "- $\\Sigma^2$ is a square diagonal matrix with the eigenvalues of $AA^T$.\n",
    "\n",
    "This means they are unique set of points that can be scaled by the covariance matrix but not rotated.  Therefore, they represent a new set of axes for our space.  This means our data can be shifted form our current axes or perspective to a new coordinate system - we can move our data points from samples in our feature space to samples in some underlying concept space. Pease refer to the **notebook about PCA** for an in-depth look into the math behind this.\n",
    "\n",
    "### Further reading\n",
    "\n",
    "- [Linear algebra explained in four pages](http://cnd.mcgill.ca/~ivan/miniref/linear_algebra_in_4_pages.pdf)\n",
    "- [Latent Semantic Analysis](https://en.wikipedia.org/wiki/Latent_semantic_indexing) on Wikipedia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}