pca/week4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 4 Assessment: Principal Component Analysis (PCA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Objective\n",
    "\n",
    "In this notebook, we will implement PCA. We will implement the two versions of PCA as described in the lectures, which handles the when the dataset size exceeds the dataset dimensionality, as well as the case when we have the dimensionality greater than the size of the dataset.\n",
    "\n",
    "We will break down the task of implementing PCA into small components and combine them in the end to produce the final algorithm. We will apply PCA to the MNIST dataset and observe how the reconstruction changes as we change the number of principal components used.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import numpy as np\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import matplotlib as mpl\n",
    "mpl.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's plot a digit from the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "from sklearn.datasets import fetch_mldata\n",
    "MNIST = fetch_mldata('MNIST original', data_home='./MNIST')\n",
    "%matplotlib inline\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.imshow(MNIST.data[0].reshape(28,28), cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we implement PCA, we will need to do some data preprocessing. In this assessment, some of them \n",
    "will be implemented by you, others we will take care of. However, when you are working on real world problems, you will need to do all these steps by yourself!\n",
    "\n",
    "The preprocessing steps we will do are\n",
    "1. Convert unsigned interger 8 (uint8) encoding of pixels to a floating point number between 0-1.\n",
    "2. Subtract from each image the mean $\\mu$.\n",
    "3. Scale each dimension of each image by $\\frac{1}{\\sigma}$ where $\\sigma$ is the standard deviation of this dimension across the whole dataset.\n",
    "\n",
    "The steps above ensure that our images will have zero mean and one variance. These preprocessing\n",
    "steps are also known as [Data Normalization or Feature Scaling](https://en.wikipedia.org/wiki/Feature_scaling)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. PCA\n",
    "\n",
    "Now we will implement PCA. Before we do that, let's pause for a moment and\n",
    "think about the steps for performing PCA. Assume that we are performing PCA on\n",
    "some dataset $\\boldsymbol X$ for $M$ principal components. \n",
    "We then need to perform the following steps, which we break into parts:\n",
    "\n",
    "1. Data normalization (`normalize`).\n",
    "2. Find eigenvalues and corresponding eigenvectors for the covariance matrix $\\boldsymbol S$.\n",
    "   Sort by the largest eigenvalues and the corresponding eigenvectors (`eig`).\n",
    "\n",
    "After these steps, we can then compute the projection and reconstruction of the data onto the spaced spanned by the top $M$ eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def normalize(X):\n",
    "    \"\"\"Normalize the given dataset X\n",
    "    Args:\n",
    "        X: ndarray, dataset\n",
    "    \n",
    "    Returns:\n",
    "        (Xbar, mean, std): ndarray, Xbar is the normalized dataset\n",
    "        with mean 0 and standard deviation 1; mean and std are the \n",
    "        mean and standard deviation respectively.\n",
    "    \n",
    "    Note:\n",
    "        You will encounter dimensions where the standard deviation is\n",
    "        zero, for those when you do normalization the normalized data\n",
    "        will be NaN. Handle this by setting using `std = 1` for those \n",
    "        dimensions when doing normalization.\n",
    "    \"\"\"\n",
    "    mu = np.mean(axis=0) # EDIT THIS\n",
    "    std = np.std(X, axis=0)\n",
    "    std_filled = std.copy()\n",
    "    std_filled[std==0] = 1.\n",
    "    Xbar = (X-mu)/std_filled                # EDIT THIS\n",
    "    return Xbar, mu, std\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def eig(S):\n",
    "    \"\"\"Compute the eigenvalues and corresponding eigenvectors \n",
    "        for the covariance matrix S.\n",
    "    Args:\n",
    "        S: ndarray, covariance matrix\n",
    "    \n",
    "    Returns:\n",
    "        (eigvals, eigvecs): ndarray, the eigenvalues and eigenvectors\n",
    "\n",
    "    Note:\n",
    "        the eigenvals and eigenvecs SHOULD BE sorted in descending\n",
    "        order of the eigen values\n",
    "        \n",
    "        Hint: take a look at np.argsort for how to sort in numpy.\n",
    "    \"\"\"\n",
    "    eig_val,eig_vec=np.linalg.eig(S)\n",
    "    idx=eig_val.argsort()[::-1]\n",
    "    return eig_val[idx],eig_vec[:,idx] # EDIT THIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def projection_matrix(B):\n",
    "    \"\"\"Compute the projection matrix onto the space spanned by `B`\n",
    "    Args:\n",
    "        B: ndarray of dimension (D, M), the basis for the subspace\n",
    "    \n",
    "    Returns:\n",
    "        P: the projection matrix\n",
    "    \"\"\"\n",
    "    P = B@(np.linalg.inv(B.T@B))@B.T # EDIT THIS\n",
    "    return P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, with the help of the functions you have implemented above, let's implement PCA! __When you implement `PCA`, do take advantage of the functions that you have implemented above.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def PCA(X, num_components):\n",
    "    \"\"\"\n",
    "    Args:\n",
    "        X: ndarray of size (N, D), where D is the dimension of the data,\n",
    "           and N is the number of datapoints\n",
    "        num_components: the number of principal components to use.\n",
    "    Returns:\n",
    "        X_reconstruct: ndarray of the reconstruction\n",
    "        of X from the first `num_components` principal components.\n",
    "    \"\"\"\n",
    "    # Compute the data covariance matrix S\n",
    "    S = np.cov(X.T)\n",
    "\n",
    "    # Next find eigenvalues and corresponding eigenvectors for S by implementing eig().\n",
    "    eig_vals, eig_vecs = eig(S)\n",
    "    \n",
    "    # Reconstruct the images from the lowerdimensional representation\n",
    "    # To do this, we first need to find the projection_matrix (which you implemented earlier)\n",
    "    # which projects our input data onto the vector space spanned by the eigenvectors\n",
    "    P = projection_matrix(eig_vecs[:,:num_components]) # projection matrix\n",
    "    \n",
    "    # Then for each data point x_i in the dataset X \n",
    "    #   we can project the original x_i onto the eigenbasis.\n",
    "    X_reconstruct = (P@X.T).T\n",
    "    return X_reconstruct\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Some preprocessing of the data\n",
    "NUM_DATAPOINTS = 1000\n",
    "X = (MNIST.data.reshape(-1, 28 * 28)[:NUM_DATAPOINTS]) / 255.\n",
    "Xbar, mu, std = normalize(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The greater number of of principal components we use, the smaller will our reconstruction\n",
    "error be. Now, let's answer the following question: \n",
    "\n",
    "\n",
    "> How many principal components do we need\n",
    "> in order to reach a Mean Squared Error (MSE) of less than $100$ for our dataset?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mse(predict, actual):\n",
    "    return np.square(predict - actual).sum(axis=1).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "loss = []\n",
    "reconstructions = []\n",
    "for num_component in range(1, 100):\n",
    "    reconst = PCA(Xbar, num_component)\n",
    "    error = mse(reconst, Xbar)\n",
    "    reconstructions.append(reconst)\n",
    "    # print('n = {:d}, reconstruction_error = {:f}'.format(num_component, error))\n",
    "    loss.append((num_component, error))\n",
    "\n",
    "reconstructions = np.asarray(reconstructions)\n",
    "reconstructions = reconstructions * std + mu # \"unnormalize\" the reconstructed image\n",
    "loss = np.asarray(loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also put these numbers into perspective by plotting them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(loss[:,0], loss[:,1]);\n",
    "ax.axhline(100, linestyle='--', color='r', linewidth=2)\n",
    "ax.xaxis.set_ticks(np.arange(1, 100, 5));\n",
    "ax.set(xlabel='num_components', ylabel='MSE', title='MSE vs number of principal components');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But _numbers don't tell us everything_! Just what does it mean _qualitatively_ for the loss to decrease from around\n",
    "$450.0$ to less than $100.0$?\n",
    "\n",
    "Let's find out! In the next cell, we draw the original eight as the leftmost image. Then we show the reconstruction of the image on the right, in descending number of principal components used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "@interact(image_idx=(0, 1000))\n",
    "def show_num_components_reconst(image_idx):\n",
    "    fig, ax = plt.subplots(figsize=(20., 20.))\n",
    "    actual = X[image_idx]\n",
    "    x = np.concatenate([actual[np.newaxis, :], reconstructions[:, image_idx]])\n",
    "    ax.imshow(np.hstack(x.reshape(-1, 28, 28)[np.arange(10)]),\n",
    "              cmap='gray');\n",
    "    ax.axvline(28, color='orange', linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also browse throught the reconstructions for other digits. Once again, `interact` becomes handy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "@interact(i=(0, 10))\n",
    "def show_pca_digits(i=1):\n",
    "    plt.figure(figsize=(4,4))\n",
    "    actual_sample = X[i].reshape(28,28)\n",
    "    reconst_sample = (reconst[i, :] * std + mu).reshape(28, 28)\n",
    "    plt.imshow(np.hstack([actual_sample, reconst_sample]), cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. PCA for high-dimensional datasets\n",
    "\n",
    "Sometimes, the dimensionality of our dataset may be larger than the number of data points we\n",
    "have. Then it might be inefficient to perform PCA with the implementation above. Instead,\n",
    "as mentioned in the lectures, we can implement PCA in a more efficient manner, which we\n",
    "call PCA for high-dimensional data (PCA_high_dim)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the __normalized__ data matrix $\\boldsymbol{\\bar{X}}$ of size $N \\times D$ where $D > N$. To do PCA we perform the following steps:\n",
    "\n",
    "* We solve the following eigenvalue/eigenvector equation for the matrix $\\frac{1}{N} \\boldsymbol{\\bar{X}} \\boldsymbol{\\bar{X}}^T$, i.e. we solve for $\\lambda_i$, $\\boldsymbol c_i$ in\n",
    "$$\\frac{1}{N} \\boldsymbol{\\bar{X}} \\boldsymbol{\\bar{X}}^T \\boldsymbol c_i = \\lambda_i \\boldsymbol c_i.$$\n",
    "\n",
    "* We want to recover original eigenvectors $\\boldsymbol b_i$ of the data covariance matrix $\\boldsymbol S = \\frac{1}{N} \\boldsymbol{\\bar{X}^T} \\boldsymbol{\\bar{X}}$.\n",
    "\n",
    "* Left-multiply the eigenvectors $\\boldsymbol c_i$ by $\\boldsymbol{\\bar{X}}^T$ yields\n",
    "$$\\frac{1}{N} \\boldsymbol{\\bar{X}}^T \\boldsymbol{\\bar{X}} \\boldsymbol{\\bar{X}}^T \\boldsymbol c_i = \\lambda_i \\boldsymbol{\\bar{X}}^T \\boldsymbol c_i$$ and we recover $\\boldsymbol b_i=\\boldsymbol{\\bar{X}}^T \\boldsymbol c_i$ as eigenvector of $\\boldsymbol S$ with the eigenvalue $\\lambda_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "def PCA_high_dim(X, num_components):\n",
    "    \"\"\"Compute PCA for small sample size. \n",
    "    Args:\n",
    "        X: ndarray of size (N, D), where D is the dimension of the data,\n",
    "           and N is the number of data points in the training set. You may assume the input \n",
    "           has been normalized.\n",
    "        num_components: the number of principal components to use.\n",
    "    Returns:\n",
    "        X_reconstruct: (N, D) ndarray. the reconstruction\n",
    "        of X from the first `num_components` principal components.\n",
    "    \"\"\"\n",
    "    N, D = X.shape\n",
    "    M = (1/N)*(X@X.T) # EDIT THIS, compute the matrix \\frac{1}{N}XX^T.\n",
    "    eig_vals, eig_vecs = eig(M) # EDIT THIS, compute the eigenvalues. \n",
    "    U = (X.T@eig_vecs)[:,:num_components] # EDIT THIS. Compute the eigenvectors for the original PCA problem.\n",
    "    # Similar to what you would do in PCA, compute the projection matrix,\n",
    "    # then perform the projection.\n",
    "    P = projection_matrix(U) # projection matrix\n",
    "    X_reconstruct = (P@X.T).T # EDIT THIS.\n",
    "    return X_reconstruct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the same dataset, `PCA_high_dim` and `PCA` should give the same output. \n",
    "Assuming we have implemented `PCA` correctly, we can then use `PCA` to test the correctness\n",
    "of `PCA_high_dim`.\n",
    "\n",
    "We can use this __invariant__\n",
    "to test our implementation of `PCA_high_dim`, assuming that we have correctly implemented `PCA`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_almost_equal(PCA(Xbar, 2), PCA_high_dim(Xbar, 2))\n",
    "# In fact, you can generate random input dataset to verify your implementation.\n",
    "print('correct')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare the running time between `PCA` and `PCA_high_dim`.\n",
    "\n",
    "__Tips__ for running benchmarks or computationally expensive code:\n",
    "\n",
    "When you have some computation that takes up a non-negligible amount of time. Try separating\n",
    "the code that produces output from the code that analyzes the result (e.g. plot the results, comput statistics of the results). In this way, you don't have to recompute when you want to produce more analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def time(f, repeat=100):\n",
    "    times = []\n",
    "    for _ in range(repeat):\n",
    "        start = timeit.default_timer()\n",
    "        f()\n",
    "        stop = timeit.default_timer()\n",
    "        times.append(stop-start)\n",
    "    return np.mean(times), np.std(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "times_mm0 = []\n",
    "times_mm1 = []\n",
    "\n",
    "for datasetsize in np.arange(4, 784, step=20):\n",
    "    XX = Xbar[:datasetsize]\n",
    "    mu, sigma = time(lambda : XX.T @ XX)\n",
    "    times_mm0.append((datasetsize, mu, sigma))\n",
    "    \n",
    "    mu, sigma = time(lambda : XX @ XX.T)\n",
    "    times_mm1.append((datasetsize, mu, sigma))\n",
    "    \n",
    "times_mm0 = np.asarray(times_mm0)\n",
    "times_mm1 = np.asarray(times_mm1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set(xlabel='size of dataset', ylabel='running time')\n",
    "bar = ax.errorbar(times_mm0[:, 0], times_mm0[:, 1], times_mm0[:, 2], label=\"$X^T X$ (PCA)\", linewidth=2)\n",
    "ax.errorbar(times_mm1[:, 0], times_mm1[:, 1], times_mm1[:, 2], label=\"$X X^T$ (PCA_high_dim)\", linewidth=2)\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first benchmark the time taken to compute $\\boldsymbol X^T\\boldsymbol X$ and $\\boldsymbol X\\boldsymbol X^T$. Jupyter's magic command `%time` is quite handy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we benchmark `PCA`, `PCA_high_dim`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "times0 = []\n",
    "times1 = []\n",
    "\n",
    "for datasetsize in np.arange(4, 784, step=100):\n",
    "    XX = Xbar[:datasetsize]\n",
    "    npc = 2\n",
    "    mu, sigma = time(lambda : XX.T @ XX)\n",
    "    times0.append((datasetsize, mu, sigma))\n",
    "    \n",
    "    mu, sigma = time(lambda : XX @ XX.T)\n",
    "    times1.append((datasetsize, mu, sigma))\n",
    "    \n",
    "times0 = np.asarray(times0)\n",
    "times1 = np.asarray(times1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, use the `time` magic command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%time Xbar.T @ Xbar\n",
    "%time Xbar @ Xbar.T\n",
    "pass # Put this here, so that our output does not show the result of computing `Xbar @ Xbar.T`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare the running time for __PCA__ and __PCA\\_high\\_dim__ directly. Spend some time and think about what this plot means. We mentioned in lectures that `PCA_high_dim` are advantageous when\n",
    "we have dataset size $N$ < data dimension $D$. Although our plot for the two running times does not intersect exactly at $N = D$, it does show the trend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set(xlabel='number of datapoints', ylabel='run time')\n",
    "ax.errorbar(times0[:, 0], times0[:, 1], times0[:, 2], label=\"PCA\", linewidth=2)\n",
    "ax.errorbar(times1[:, 0], times1[:, 1], times1[:, 2], label=\"PCA_high_dim\", linewidth=2)\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, with the magic command `time`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%time PCA(Xbar, 2)\n",
    "%time PCA_high_dim(Xbar, 2)\n",
    "pass"
   ]
  }
 ],
 "metadata": {
  "coursera": {
   "course_slug": "mathematics-machine-learning-pca",
   "graded_item_id": "CXC11",
   "launcher_item_id": "ub5A7"
  },
  "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}