02_Lab_2_Principal_component_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Mining and Exploration [INFR11007]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 2: Principal component analysis (PCA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Principal component analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis) is one of the most used techniques for exploratory data analysis and preprocessing.\n",
    "\n",
    "We start this lab off by introducing some fundamental concepts regarding covariance matrices. We then look at several methods for identifying the principal component (PC) directions and scores for an artificial dataset. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by importing the packages we will need throughout the lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Covariance matrices\n",
    "Covariance matrices are by definition [symmetric positive semi-definite](https://en.wikipedia.org/wiki/Positive-definite_matrix#Positive-semidefinite). One property of positive semi-definite matrices is that their eigenvalues are all non-negative. You are provided with the function  `generate_spsd_matrix()` which generates a random symmetric positive semi-definite matrix. The `random_seed` parameter can be used to set the numpy random seed in order to ensure reproducible results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def generate_spsd_matrix(d, random_seed=None):\n",
    "    \"\"\"Reproducible random generation of symmetric \n",
    "    positive semi-definite matrices. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    d : integer\n",
    "        Number of dimensions.\n",
    "    \n",
    "    random_seed : integer\n",
    "        Random seed number.\n",
    "        \n",
    "    Returns \n",
    "    ----------\n",
    "    A : ndarray, shape (n,n)\n",
    "        Symmetric positive definite matrix.\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    if random_seed is not None:\n",
    "        np.random.seed(random_seed)\n",
    "    A = np.random.randn(d,d)\n",
    "    return np.dot(A.T, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1 ==========\n",
    "\n",
    "Write a function which checks if a matrix is positive semi-definite. \n",
    "\n",
    "Generate a random 10 x 10 matrix by using the `generate_spd_matrix()` function and double-check that it is positive semi-definite by using the function you just wrote.\n",
    "\n",
    "*Hint: you might find the [`numpy.linalg.eigh()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html) function useful.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sampling data with a given covariance matrix\n",
    "### Covariance matrix eigendecomposition\n",
    "\n",
    "Suppose that you want to sample (i.e. generate) some data points from a multivariate normal distribution with known mean vector and covariance matrix. We will assume that we have access to the tool required to sample from a standard univariate normal distribution, that is, one with zero mean and unit variance (all you need for that is just the cumulative function of the normal distribution). \n",
    "\n",
    "We know that if a univariate gaussian random variable $x$ has zero mean and unit variance, then $y = \\sigma x + \\mu$ has mean $\\mu$ and variance $\\sigma^2$ . Therefore, we can use this property to sample from an arbitrary univariate normal distribution with mean $\\mu$ and variance $\\sigma^2$.\n",
    "\n",
    "By repeating the process described above $d$ times, we can sample from a $d$-dimensional normal distribution with an arbitrary mean vector and diagonal covariance matrix. But, how can we sample from a multivariate normal distribution with given covariance matrix $\\mathbf{C}$? The answer is through decomposing the covariance matrix. One such method is the [eigen (or spectral) decomposition](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix) which, for a real symmetric matrix (such as covariance matrices), takes the following form:\n",
    "$$\\mathbf{C} = \\mathbf{U} \\mathbf{\\Lambda} \\mathbf{U}^T$$\n",
    "where $\\mathbf{U}$ is an orthogonal matrix containing the (right) eigenvectors of $\\mathbf{C}$ and $\\mathbf{\\Lambda}$ is a diagonal matrix whose entries are the eigenvalues of $\\mathbf{C}$. Now, you might wonder how this decomposition can help us sample from a multivariate distribution.\n",
    "\n",
    "Assume $\\mathbf{x}$ is sampled from a standard multivariate normal distribution (i.e. zero mean vector, identity covariance matrix). We can now create a new random variable $\\mathbf{y}$ by transforming $\\mathbf{x}$ as follows:\n",
    "$$\\mathbf{y} =  \\mathbf{U} \\mathbf{\\Lambda}^{1/2} \\mathbf{x} + \\mathbf{\\mu}$$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2 ==========\n",
    "\n",
    "Assume that $\\mathbf{x}$ is a multivariate random variable with zero mean and unit covariance matrix, and $\\mathbf{C} = \\mathbf{U} \\mathbf{\\Lambda} \\mathbf{U}^T$ is the eigendecomposition of $\\mathbf{C}$.\n",
    "\n",
    "If $\\mathbf{y} =  \\mathbf{U} \\mathbf{\\Lambda}^{1/2} \\mathbf{x} + \\mathbf{\\mu}$, show that  $\\mathbf{y}$ has mean $\\mathbf{\\mu}$ and covariance $\\mathbf{C} $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that the above holds without any assumptions on the probability distribution of $\\mathbf{x}$. As a special case, we focus on the multivariate normal distribution. \n",
    "\n",
    "The procedure of sampling one data point from a multivariate normal distribution with mean $\\mathbf{\\mu}$ and covariance matrix $\\mathbf{C}$ is as follows:\n",
    "\n",
    "1. Decompose the covariance matrix by using its eigendecompostion, i.e. $\\mathbf{C} = \\mathbf{U} \\mathbf{\\Lambda} \\mathbf{U}^T$\n",
    "2. Sample a data point $\\mathbf{x} \\in \\mathbb{R}^d$ from a standard normal distribution (i.e. zero mean, identity covariance matrix).\n",
    "3. Transform $\\mathbf{x}$ by using $\\mathbf{y} =  \\mathbf{U} \\mathbf{\\Lambda}^{1/2} \\mathbf{x} + \\mathbf{\\mu}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 3 ==========\n",
    "\n",
    "Write a function that generates random samples from a multivariate normal distribution with given mean vector and  covariance matrix. You should make use of the [`numpy.random.randn()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.randn.html) function that generates samples from a standard multivariate gaussian distribution and the eigendecomposition of the covariance matrix. For computing the eigendecomposition of a symmetric matrix you should use the [`numpy.linalg.eigh()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html) function.\n",
    "\n",
    "**Important:**\n",
    "\n",
    "Most scientific computing packages like `scikit-learn` follow the $n \\times d$ convention for a data matrix, where $n$ is the number of samples and $d$ the dimensionality, that is, $\\mathbf{X} \\in \\mathbb{R} ^ {n \\times d}$. Please note that many textbooks and the lecture notes for this course follow reverse notation (i.e. $d \\times n$). Here, you should follow the $n \\times d$ convention. If the standard normal data are stored in a matrix $\\mathbf{X}$, then you should compute $\\mathbf{Y}$ by $$\\mathbf{Y} = \\mathbf{X} \\mathbf{\\Lambda}^{\\frac{1}{2}} \\mathbf{U}^T $$\n",
    "\n",
    "Finally, generate a 3 x 3 random covariance matrix `C` by using the `generate_positive_semi_definite()` function.  Use the function you just wrote to generate 1 million random samples with zero mean and covariance matrix `C`. Compute the empirical covariance matrix of the data (you can use [`numpy.cov()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html) by setting `rowvar=False`) and double-check that it is a good estimate of the true covariance matrix `C`.\n",
    "\n",
    "*Tip: don't get stuck on this, move on after 10 mins or so (it's not critical)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note: numpy built-in function\n",
    "Numpy implements the [`numpy.random.multivariate_normal()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.multivariate_normal.html) function that can be used to generate samples from a multivariate normal distribution with given mean vector and covariance matrix. You are encouraged to use such built-in functions whenever available, as they will most likely be highly optimised, and bug-free. Nevertheless, it is some times very useful to know what these functions do under the hood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Covariance matrix structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are provided with the following function `generate_gaussian_data()` that can be used to generate a multivariate gaussian dataset from a given mean and covariance. When the mean and covariance are not defined, they are generated at random. The `random_seed` parameter can be used to ensure reproducible results. The function returns a tuple containing three items; the dataset, the true mean, and the true covariance matrix of the probability distribution the data were sampled from. Execute the cell below to load this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def generate_gaussian_data(n_samples, n_features=None, mu=None, cov=None, random_seed=None):\n",
    "    \"\"\"Generates a multivariate gaussian dataset.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    n_samples : integer\n",
    "        Number of samples.\n",
    "        \n",
    "    n_features : integer\n",
    "        Number of dimensions (features).\n",
    "        \n",
    "    mu : array, optional (default random), shape (n_features,)\n",
    "        Mean vector of normal distribution.\n",
    "        \n",
    "    cov : array, optional (default random), shape (n_features,n_features)\n",
    "        Covariance matrix of normal distribution.\n",
    "    \n",
    "    random_seed : integer\n",
    "        Random seed.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    x : array, shape (n_samples, n_features)\n",
    "        Data matrix arranged in rows (i.e. \n",
    "        columns correspond to features and \n",
    "        rows to observations).\n",
    "    \n",
    "    mu : array, shape (n_features,)\n",
    "        Mean vector of normal distribution.\n",
    "        \n",
    "    cov : array, shape (n_features,n_features)\n",
    "        Covariance matrix of normal distribution.\n",
    "        \n",
    "    Raises\n",
    "    ------\n",
    "    \n",
    "    ValueError when the shapes of mu and C are not compatible\n",
    "    with n_features.\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    if random_seed is not None:\n",
    "        np.random.seed(random_seed)\n",
    "        \n",
    "    if mu is None:\n",
    "        mu = np.random.randn(n_features,)\n",
    "    else:\n",
    "        if n_features is None:\n",
    "            n_features = mu.shape[0]\n",
    "        else:\n",
    "            if mu.shape[0] != n_features:\n",
    "                raise ValueError(\"Shape mismatch between mean and number of features.\")\n",
    "                \n",
    "    if cov is None:\n",
    "        cov = generate_spsd_matrix(n_features, random_seed=random_seed)\n",
    "    else:\n",
    "        if (cov.shape[0] != n_features) or (cov.shape[1] != n_features):\n",
    "            raise ValueError(\"Shape mismatch between covariance and number of features.\")\n",
    "            \n",
    "    x = np.random.multivariate_normal(mu, cov, n_samples)\n",
    "    return (x, mu, cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 4 ==========\n",
    "\n",
    "By using the function provided above, generate a 2-dimensional gaussian dataset consisting of 1000 samples. Use a zero mean and the following covariance matrix:\n",
    "\n",
    "$$ \\mathbf{C} = \\left( \\begin{array}{ccc} \n",
    "1 & 0.3 \\\\\n",
    "0.3 & 2 \\\\ \\end{array} \\right) $$\n",
    "\n",
    "Print the empirical mean, covariance and correlation matrices. You should use numpy built-in functions for this purpose, but you are also free to implement your own functions if you are keen. \n",
    "Finally, use the seaborn [`jointplot()`](http://seaborn.pydata.org/generated/seaborn.jointplot.html) function to produce a joint scatter plot of the two variables. This function also shows the marginal histograms on the top and right hand sides of the plot. Label axes appropriately. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 5 ==========\n",
    "\n",
    "Repeat the same procedure as above, but now modify the covariance matrix such that the **true correlation ceofficient** between the two variables is 0.6 while the variances stay unchanged. \n",
    "\n",
    "Visualise the `jointplot()` of the new dataset and compare to the one from the previous question. Do your observations agree with the modification you just applied to the covariance matrix?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal component analysis (PCA)\n",
    "\n",
    "Let us now dive into PCA. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this section, we will use a 3-dimensional Gaussian dataset. Execute the cell below to generate the dataset and print the true mean and covariance matrix of the distribution the data was sampled from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Generates a 3D dataset and prints true mean and covariance\n",
    "x_3d, mu_true, C_true = generate_gaussian_data(n_samples=1000, n_features=3, random_seed=10)\n",
    "print(\"Dataset consists of {} samples and {} features.\".format(x_3d.shape[0], x_3d.shape[1]))\n",
    "print(\"True mean:\\n{}\".format(mu_true))\n",
    "print(\"True covariance matrix:\\n{}\".format(C_true))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 6 ==========\n",
    "\n",
    "Estimate the empirical mean and covariance matrix from this dataset and save them in variables `mu_est` and `C_est` respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA sklearn implementation\n",
    "\n",
    "Sklearn offers a class implementation of `pca`. Please spend a minute to read the [documentation](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html) of this class. The principal component (PC) directions of a dataset are computed by using the [`fit()`](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html#sklearn.decomposition.PCA.fit) method and stored into the `components_` attribute. The PC scores can be computed by using the [`transform()`](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html#sklearn.decomposition.PCA.transform) method. The amount of variance explained by each of the selected components is stored into the `explained_variance_` attribute. \n",
    "\n",
    "Please note, this method centers the input data (i.e. removes the mean) before computing the PC directions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 7 ==========\n",
    "\n",
    "Initialise a `pca` object and \"fit\" it using the dataset `x_3d`. Print the three PC directions. Store the PC scores for `x_3d` in an array called `pc_scores`.\n",
    "\n",
    "*Hint: according to the documentation the components are stored row-wise, i.e. the `components_ array` has shape `[n_components, n_features]`. You should print the PC directions column-wise, i.e. take the transpose of the `components_` array.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most often, we do not want to compute all PC directions, but only a few (i.e. dimensionality reduction). We can define the desired number of PCs by setting the [`n_components`]() parameter appropriately when we initialise the `pca` object. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 8 ==========\n",
    "\n",
    "Initialise a `pca_new` object with 2 PCs and \"fit\" it using the dataset `x_3d`. Print the two PC directions. \n",
    "\n",
    "*Hint: the 2 PC directions should be the same as the first 2 directions you computed in the previous question. The reason for this is ultimately that PCA by sequential and simultaneous variance maximisation give the same result.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA implementation #1: covariance matrix eigendecomposition\n",
    "\n",
    "Now we want to implement PCA from scratch. The procedure can be summarised as follows:\n",
    "\n",
    "1. Center the data (i.e. remove the mean).\n",
    "2. Compute the empirical covariance matrix.\n",
    "3. Perform the eigendecomposition of the estimated covariance matrix.\n",
    "4. Sort eigenvalues and associated eigenvectors, in eigenvalue descending order. The sorted eigenvectors correspond to the PC directions. If we want to reduce the dimensionality, we select the first `k` eigenvectors corresponding to the `k` largest eigenvalues (`k` < `d`).\n",
    "5. To compute PC scores we project the centered data matrix (i.e. matrix product) onto the PC directions.\n",
    "\n",
    "(Regarding 3 and 4: some algorithms for eigendecompositions allow you to specify the number of eigenvectors and eigenvalues that should be computed. You then do not have to compute the complete eigendecomposition, which is wasteful if you are only interested in a few principle components.)\n",
    "\n",
    "### ========== Question 9 ==========\n",
    "\n",
    "Compute and print all three PC directions in the dataset `x_3d` by using the procedure described above. Hopefully, you will have already performed the first two steps in question 6. Then compute the PC scores for the centered data matrix.\n",
    "\n",
    "As happens very often when writing code, it is likely that there will be a few bugs in your implementation. To double-check that your code is correct, compare the computed PC directions and scores to the ones achieved in question 7 with scikit-learn.\n",
    "\n",
    "*Hint: you might (or might not) find that some of the PC directions/scores you have computed have opposite signs to the ones returned by the sklearn implementation. Do not worry about this, the two solutions are equivalent (why?). To make debugging easier, you are provided with the following function, `solutions_equivalent()` which tests wether two solutions are equivalent, regardless of their signs. Execute the following cell to load this function.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def solutions_equivalent(b1, b2):\n",
    "    \"\"\"Checks whether two PC directions/scores\n",
    "    solutions are equivalent regardless of their .\n",
    "    respective signs.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    s1 : array, shape(n_axes,)\n",
    "        First solution.\n",
    "        \n",
    "    s2 : array, shape(n_axes,)\n",
    "        Second solution.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    True if solutions are equivalent.\n",
    "    \n",
    "    Raises\n",
    "    ------\n",
    "    \n",
    "    ValueError if the two bases do not have\n",
    "    the same dimensionality.\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    s1 = np.asarray(b1)\n",
    "    s2 = np.asarray(b2)\n",
    "    \n",
    "    if s1.shape != s2.shape:\n",
    "        raise ValueError(\"Solutions must have the same dimensionality.\")\n",
    "    \n",
    "    dims_equal = []\n",
    "    for dim in xrange(s1.shape[1]):\n",
    "        if (np.allclose(s1[:,dim],s2[:,dim]) or np.allclose(s1[:,dim],-s2[:,dim])):\n",
    "            dims_equal.append(True)\n",
    "        else:\n",
    "            dims_equal.append(False)\n",
    "    return all(element for element in dims_equal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA implementation #2: data matrix singular value decomposition (SVD)\n",
    "\n",
    "### ========== Question 10 ==========\n",
    "\n",
    "The SVD decomposition of the centered data matrix is: $$\\mathbf{X} =  \\mathbf{V} \\mathbf{S} \\mathbf{U}^T  $$ where the columns $\\mathbf{V} \\in \\mathbb{R}^{n\\times n}$ and $\\mathbf{U} \\in \\mathbb{R}^{d\\times d}$ contain the left- and right-singular vectors of $\\mathbf{X}$, and $\\mathbf{S} \\in \\mathbb{R}^{n\\times d}$ is a diagonal matrix containing its singular values.\n",
    "\n",
    "Compute the PC directions and scores in the dataset `x_3d` by using the [SVD](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html) decomposition of the centered data matrix (cf. lecture notes). Double-check that these are correct by using the `solutions_equivalent()` function.\n",
    "\n",
    "*Hint 1: make sure you center the data before computing the SVD of `x_3d`.*\n",
    "\n",
    "*Hint 2: Remember, we are dealing with `n` x `d` data matrices here, as opposed to `d` x `n` matrices in the lecture notes (that is why we have written $\\mathbf{X} =  \\mathbf{V} \\mathbf{S} \\mathbf{U}^T$ as opposed to $\\mathbf{X} =  \\mathbf{U} \\mathbf{S} \\mathbf{V}^T$ in the lecture notes). Which singular vectors of the centered data matrix correspond to the PC directions/scores in this case?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume now that you want to do dimensionality reduction by using `k` PCs only. One option would be to follow the same procedure as above and discard the last `d-k` PC directions. This, however, can be rather wasteful. There exist [SVD variations](https://en.wikipedia.org/wiki/Low-rank_approximation) that approximate the matrix $\\mathbf{X}$ with a low-rank matrix. \n",
    "\n",
    "Scipy implements the function [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) that computes only the `k` largest singular values/vectors of a given matrix. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 11 ==========\n",
    "\n",
    "Compute and print the first 2 PC directions by using the [`svds`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.svds.html) decomposition of the centered data matrix.\n",
    "\n",
    "*Hint: as opposed to the numpy function `numpy.linalg.svd`, `svds` does not automatically sort the singular values/vectors in singular value descending order. You should do this operation manually.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import svds\n",
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA implementation #3: Gram matrix eigendecomposition\n",
    "\n",
    "We have seen in the lectures that we can compute the PC directions and scores by performing the eigendecomposition of the [Gram matrix](https://en.wikipedia.org/wiki/Gramian_matrix) $\\mathbf{G}$.\n",
    "\n",
    "The procedure of computing the PC directions and scores based on the eigendecomposition of the Gram matrix is as follows:\n",
    "1. Center data.\n",
    "2. Compute Gram matrix.\n",
    "3. Perform the eigendecomposition of the Gram matrix.\n",
    "4. Sort eigenvalues and eigenvectors in decreasing eigenvalue order.\n",
    "5. Pick the first `k` eigenvectors. This is required even if we do not to perform dimensionality reduction, since the eigendecomposition of the Gram matrix will yield `n` eigenvectors.\n",
    "6. The PC scores are given by these eigenvectors scaled by $\\mathbf{\\Sigma} ^ {1/2}$, where $\\mathbf{\\Sigma}$ is the diagonal matrix containing the eigenvalues of the gram matrix.\n",
    "7. To get the PC directions, we project the centered data matrix onto the eigenvectors and scale by $\\mathbf{\\Sigma} ^ {-1/2}$.\n",
    "\n",
    "### ========== Question 12 ==========\n",
    "\n",
    "Compute the PC directions and scores in the dataset `x_3d` by using the eigendecomposition of the Gram matrix. Double-check that these are correct by using the `solutions_equivalent()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 13 ==========\n",
    "\n",
    "Can you think of a particular case where it would be beneficial to compute the PC scores by using the eigendecomposition of the Gram matrix, as opposed to SVD decomposition of the data matrix?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "One particular case where we would use the eigendecomposition of the Gram matrix, is when we only have access to that matrix and not the data matrix. In this case, it is still possible to identify the PC scores and hence do dimensionality reduction, as we will see in Lab 3.\n",
    "\n",
    "In addition, by working directly with the Gram matrix, we can use the kernel trick and identify non-linear manifolds by using almost the same bit of code. This is exactly what [kernel principal component analysis (kPCA)](https://en.wikipedia.org/wiki/Kernel_principal_component_analysis) does, but more on this in the next lab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimensionality reduction with PCA\n",
    "\n",
    "PCA can be used to reduce the dimensionality of a data from $d$ to $k$.\n",
    "\n",
    "### ========== Question 14 ==========\n",
    "\n",
    "Use any implementation of PCA of your choice to project the dataset `x_3d` onto 1, 2, and 3 dimensions. This projection should just take the form $$\\mathbf{X}_{red} = \\mathbf{X} \\mathbf{W}_{1:k}$$ where $\\mathbf{X}$ is the centered data matrix and $\\mathbf{W}$ is the matrix containing the PC directions as columns. The approximation of the data by the $k$ principle components is $$\\mathbf{X}_{recon} = \\mathbf{X}_{red} \\mathbf{W}_{1:k} ^ T + \\mathbf{\\mu}$$.\n",
    "\n",
    "Plot the mean-squared-error and explained variance of the approximation as a function of the number of components used. Label axes appropriately. You should make use of the sklearn [`mean_squared_error()`](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_error.html) and [`explained_variance_score()`](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.explained_variance_score.html) metrics. For `explained_variance_score()`, you should set the `multioutput` parameter to `variance_weighted`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The black points lie on a line (the 1st PC direction). The red points lie in a 2D plane (spanned by the first two PC directions)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 15 ==========\n",
    "\n",
    "How can you compute the variance explained by k principal components by only looking at the eigenvalues of the covariance matrix?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Low-rank SVD approximation for image compression [optional]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In lecture, we have seen that the SVD allows us to find a low rank approximation of the data matrix. We here exemplify the low rank approximation property of the SVD on a image compression task.\n",
    "\n",
    "As you might already know, grey-scale images are represented in the digital world as 2D matrices, whose elements correspond to pixel intensities. We here approximate this matrix by a low-rank approximation through the SVD. If there are correlations between the pixels in the image (which happens to be the case for [natural images](http://www.naturalimagestatistics.net/)), then we should be able to achieve a relatively good reconstruction of the image by using only a few components.\n",
    "\n",
    "Let us first load a sample image from the scipy package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Load sample image\n",
    "from scipy.misc import face\n",
    "img = face(gray=True)\n",
    "print(\"Image array dimensionality: {}\".format(img.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualise the image by using the matplotlib imshow function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Show image \n",
    "sns.set_style(\"white\")\n",
    "plt.figure()\n",
    "plt.imshow(img, cmap=plt.cm.gray)\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 16 [optional] ==========\n",
    "\n",
    "Write a function image_low_rank_approx() that takes as input an image (i.e. 2-dimensional array) and an integer k and reconstructs the image by using a k-rank SVD approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 17 [optional] ==========\n",
    "Perform a low-rank approximation of the image stored in img by using a varying number of ranks (i.e. from 1 to 500) and visualise the approximation. What do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here.*\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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": 2
}