class13_PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Data Mining:<br>Statistical Modeling and Learning from Data\n",
    "\n",
    "## Dr. Ciro Cattuto<br>Dr. Laetitia Gauvin<br>Dr. André Panisson\n",
    "\n",
    "### Exercises - Principal Component Analysis\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MNIST dataset contains 70000 images of hand-written digits. Each image is a matrix of 28x28 grayscale pixels. The value of each pixel goes from 0 to 255.\n",
    "\n",
    "There is more information available in the MNIST dataset repository:\n",
    "http://yann.lecun.com/exdb/mnist/\n",
    "\n",
    "Let's first fetch the dataset from the internet (which may take a while, note the asterisk [*]):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_mldata\n",
    "mnist = fetch_mldata(\"MNIST Original\", data_home='.') # fetch dataset from internet\n",
    "X_all = mnist.data\n",
    "y_all = mnist.target\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset is ordered, and we will get just part of it to speed up this exercise.\n",
    "\n",
    "For this, we need to shuffle the entire data and get the first instances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.utils import shuffle\n",
    "\n",
    "# shuffle dataset (which is sorted!)\n",
    "X_all, y_all = shuffle(X_all, y_all, random_state=42)\n",
    "\n",
    "# take only the first instances, to shorten runtime\n",
    "X = X_all[:1500]\n",
    "y = y_all[:1500]\n",
    "\n",
    "# take the other instances for testing\n",
    "X_test = X_all[1500:]\n",
    "y_test = y_all[1500:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the sklearn documentation:\n",
    "\n",
    "> Many elements used in the objective function of a learning algorithm (such as the RBF kernel of Support Vector Machines or the l1 and l2 regularizers of linear models) assume that all features are centered around zero and have variance in the same order. If a feature has a variance that is orders of magnitude larger that others, it might dominate the objective function and make the estimator unable to learn from other features correctly as expected.\n",
    "\n",
    "We will then rescale the training data, since this is important to improve the performance of the classification algorithm that we will use later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rescale grayscale from -1 to 1\n",
    "X = X/255.0*2 - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at some of the instances in the dataset we just loaded:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc(\"image\", cmap=\"binary\")\n",
    "plt.figure(figsize=(8,4))\n",
    "for i in xrange(10):\n",
    "    plt.subplot(2,5,i+1)\n",
    "    plt.imshow(X[i].reshape(28,28))\n",
    "    plt.title(y[i])\n",
    "    plt.xticks(())\n",
    "    plt.yticks(())\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA (Principal Component Analysis) converts a dataset of possibly correlated features into a dataset of linearly uncorrelated features called principal components. is equivalent to **Singular Value Decomposition** (SVD):\n",
    "\n",
    "$$\\mathbf{X} = \\mathbf{U} \\boldsymbol{\\Sigma} \\mathbf{V}^T$$\n",
    "\n",
    "The value returned from `transform` is in fact the matrix $\\mathbf{U}$. The parameter `n_components` is the number of components to keep. if n_components is not set, all components are kept; it means that the result of `transform` is a matrix $n \\times n$ (if $d \\geq n$) or $n \\times d$ (if $d < n$).\n",
    "\n",
    "The values in the attribute `explained_variance_` are related to $\\boldsymbol{\\Sigma}$ in the SVD decomposition. Finally, the values in the attribute `components_` are related to $\\mathbf{V}^T$.\n",
    "\n",
    "The optional parameter `whiten=True` make it possible to project the data onto the singular space while scaling each component to unit variance. This is often useful if the models down-stream make strong assumptions on the isotropy of the signal: this is for example the case for Support Vector Machines with the RBF kernel and the K-Means clustering algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA()\n",
    "\n",
    "pca.fit(X)\n",
    "\n",
    "U = pca.transform(X)\n",
    "S = pca.explained_variance_\n",
    "V = pca.components_\n",
    "\n",
    "print (U.shape)\n",
    "print (S.shape)\n",
    "print (V.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It does not make much sense to look at the transformed images, they will look like noise to us.\n",
    "\n",
    "Instead, let's have a look at the most important directions on which the dataset was projected. These directions can be accessed in the $\\mathbf{V}^T$ matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "\n",
    "for i in xrange(10): # loop over all means\n",
    "    plt.subplot(2,5,i+1)\n",
    "    plt.imshow(V[i].reshape(28,28))\n",
    "    plt.xticks(())\n",
    "    plt.yticks(())\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The later directions (here, from the 100-th on) show noise and small variations between different, but very similar training instances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,4))\n",
    "\n",
    "for i in xrange(10): # loop over all means\n",
    "    plt.subplot(2,5,i+1)\n",
    "    plt.imshow(V[100+i].reshape(28,28))\n",
    "    plt.xticks(())\n",
    "    plt.yticks(())\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of \"interesting\" dimensions can be seen from the importance of the found directions.\n",
    "\n",
    "We can find this information by plotting the `explained_variance_` attribute from the `PCA` variable. This is related to the $\\Sigma$ matrix in the following way:\n",
    "\n",
    "$$\\mathbf{ev} = \\frac{1}{N} \\Sigma^2 $$\n",
    "\n",
    "In order to find the most relevant directions, we can use the `explained_variance_` to get 90% of the dataset variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(pca.explained_variance_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev_cumsum = np.cumsum(pca.explained_variance_)/(pca.explained_variance_).sum()\n",
    "ev_at90 = ev_cumsum[ev_cumsum<0.9].shape[0]\n",
    "print (ev_at90)\n",
    "\n",
    "plt.plot(ev_cumsum)\n",
    "plt.vlines(ev_at90, 0, 1, linestyles='dashed')\n",
    "plt.hlines(0.9, 0, 500, linestyles='dashed');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the intrinsic dimensionality is not higher than 80, even though the dataset has 784 dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By setting the parameter `n_components` to $\\mathbf{k}$, we limit the number of components to $\\mathbf{k}$. At this point, the decomposition is a **Truncated SVD**.\n",
    "\n",
    "As consequence, the value returned from `transform` is in fact a matrix $\\mathbf{U} \\in \\mathbb{R}^{n \\times k}$, where n is the number of samples and k is the number of components. The number of values in the attribute `explained_variance_`, related to $\\boldsymbol{\\Sigma}$ in the SVD decomposition, is $\\mathbf{k}$. Finally, the values in the attribute `components_` are a matrix $\\mathbf{V}^T \\in \\mathbb{R}^{k \\times m}$, where $m$ is the number of features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(80)\n",
    "\n",
    "pca.fit(X)\n",
    "\n",
    "U = pca.transform(X)\n",
    "S = pca.explained_variance_\n",
    "V = pca.components_\n",
    "\n",
    "print (U.shape)\n",
    "print (S.shape)\n",
    "print (V.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "From now on, we will try to train a SVM classificator to classify the numbers.\n",
    "\n",
    "For this, create a SVM classifier with **rbf** kernel and default parameters and assign it to a variable with name **`svm`**.\n",
    "Then, create a list named **`gammas`** with 10 values for the rbf coefficient (gamma) starting from $10^{-4}$ to $10^{0}$, spaced evenly on a log scale. For each value of gamma, set the model parameter **`gamma`** to this value, calculate a list of 3-fold cross-validation scores, and add the average score to a list with name **`scores`**. Use **`X`** and **`y`** to train the model.\n",
    "\n",
    "Finally, plot the score for each gamma value.\n",
    "\n",
    "Make note of the time that it takes to search for the best gamma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('Best gamma:', gammas[np.argmax(scores)])\n",
    "print ('Best score:', scores[np.argmax(scores)])\n",
    "\n",
    "semilogx(gammas, scores)\n",
    "plt.xlabel('gamma')\n",
    "plt.ylabel('Score (accuracy)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the score of the best model in the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svm = SVC(kernel='rbf', gamma=gammas[np.argmax(scores)])\n",
    "svm.fit(X, y)\n",
    "\n",
    "X_test_norm = X_test/255.0*2 - 1\n",
    "print (svm.score(X_test_norm, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next exercise, reduce the **`X`** data to a matrix **`U`** using PCA with 80 components.\n",
    "\n",
    "Next, reproduce exactly the same steps as the previous exercise, but use **`U`** instead of **`X`** to train the model.\n",
    "\n",
    "Make note of the time that it takes to search for the best gamma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('Best gamma:', gammas[np.argmax(scores)])\n",
    "print ('Best score:', scores[np.argmax(scores)])\n",
    "\n",
    "semilogx(gammas, scores)\n",
    "plt.xlabel('gamma')\n",
    "plt.ylabel('Score (accuracy)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the score of the best model that uses PCA in the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svm = SVC(kernel='rbf', gamma=gammas[np.argmax(scores)])\n",
    "svm.fit(U, y)\n",
    "\n",
    "X_test_norm = X_test/255.0*2 - 1\n",
    "U_test = pca.transform(X_test_norm)\n",
    "print (svm.score(U_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, find the best value also for the C parameter and check its score in the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}