Session 2 - Features_II.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](http://datascience.uci.edu/wp-content/uploads/sites/2/2014/09/data_science_logo_with_image1.png 'UCI_data_science')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Predictive Modeling with Python \n",
    "## Session #2: Feature Engineering II\n",
    "Author: [Eric Nalisnick](http://www.ics.uci.edu/~enalisni/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Schedule for Today\n",
    "\n",
    "|Start Time | Session |\n",
    "|-----------|---------|\n",
    "|8:30am     | Check In|\n",
    "|9:00am     | Feature Engineering I |\n",
    "|10:30am    | Break|\n",
    "|10:45am    | **Feature Engineering II** |\n",
    "|12:30pm    | Lunch |\n",
    "|1:00pm    | Feature Engineering III |\n",
    "|2:30pm    | Break |\n",
    "|2:45pm    | Ensembling |\n",
    "|5:00pm    | End |\n",
    "\n",
    "### Goals of this Lesson\n",
    "- Gradient Descent for PCA\n",
    "- Nonlinear Dimensionality Reduction\n",
    "    - Autoencoder: Model and Learning\n",
    "    - Autoencoding Images\n",
    "    - Denoising Autoencoder\n",
    "\n",
    "### References \n",
    "- [Ch. 14 'Autoencoders' of the *Deep Learning* book](http://www.deeplearningbook.org/contents/autoencoders.html)\n",
    "- *Science* article on Autoencoders: [*Reducing the Dimensionality of Data with Neural Networks* by Hinton and  Salakhutdinov](https://www.cs.toronto.edu/~hinton/science.pdf)\n",
    "- [SciKit-Learn's documentation on neural networks](http://scikit-learn.org/stable/modules/neural_networks_supervised.html)\n",
    "- [SciKit-Learn's documentation on dimensionality reduction](http://scikit-learn.org/stable/modules/decomposition.html#decompositions)\n",
    "\n",
    "## 0.  Preliminaries\n",
    "First we need to import Numpy, Pandas, MatPlotLib..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import time\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we need functions for shuffling the data and calculating classification errrors.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### function for shuffling the data and labels\n",
    "def shuffle_in_unison(features, labels):\n",
    "    rng_state = np.random.get_state()\n",
    "    np.random.shuffle(features)\n",
    "    np.random.set_state(rng_state)\n",
    "    np.random.shuffle(labels)\n",
    "    \n",
    "### calculate classification errors\n",
    "# return a percentage: (number misclassified)/(total number of datapoints)\n",
    "def calc_classification_error(predictions, class_labels):\n",
    "    n = predictions.size\n",
    "    num_of_errors = 0.\n",
    "    for idx in xrange(n):\n",
    "        if (predictions[idx] >= 0.5 and class_labels[idx]==0) or (predictions[idx] < 0.5 and class_labels[idx]==1):\n",
    "            num_of_errors += 1\n",
    "    return num_of_errors/n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.1  Load the dataset of handwritten digits\n",
    "We are going to use the [MNIST dataset](https://en.wikipedia.org/wiki/MNIST_database) throughout this session.  Let's load the data...  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load the 70,000 x 784 matrix\n",
    "from sklearn.datasets import fetch_mldata\n",
    "mnist = fetch_mldata('MNIST original').data\n",
    "\n",
    "# reduce to 5k instances\n",
    "np.random.shuffle(mnist)\n",
    "mnist = mnist[:5000,:]/255.\n",
    "print \"Dataset size: %d x %d\"%(mnist.shape)\n",
    "\n",
    "# subplot containing first image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "digit = mnist[1,:]\n",
    "ax1.imshow(np.reshape(digit, (28, 28)), cmap='Greys_r')\n",
    "\n",
    "# subplot containing second image\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "digit = mnist[2,:]\n",
    "ax2.imshow(np.reshape(digit, (28, 28)), cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 1 Gradient Descent for PCA\n",
    "\n",
    "Recall the [Principal Component Analysis](https://en.wikipedia.org/wiki/Principal_component_analysis) model we covered in the last session.  Again, the goal of PCA is for a given datapoint $\\mathbf{x}_{i}$, find a lower-dimensional representation $\\mathbf{h}_{i}$ such that $\\mathbf{x}_{i}$ can be 'predicted' from $\\mathbf{h}_{i}$ using a linear transformation.  Again, the loss function can be written as: $$ \\mathcal{L}_{\\text{PCA}} = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{x}_{i}\\mathbf{W}\\mathbf{W}^{T})^{2}.$$   \n",
    "Instead of using the closed-form solution we discussed in the previous session, here we'll use gradient descent.  The reason for doing this will become clear later in the session, as we move on to cover a non-linear version of PCA.  To run gradient descent, we of course need the derivative of the loss w.r.t. the parameters, which are in this case, the transformation matrix $\\mathbf{W}$:\n",
    "$$ \\nabla_{\\mathbf{W}} \\mathcal{L}_{\\text{PCA}} = -4\\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{\\tilde x}_{i})^{T}\\mathbf{h}_{i} $$\n",
    "\n",
    "\n",
    "Now let's run our stochastic gradient PCA on the MNIST dataset...\n",
    "#### <span style=\"color:red\">Caution: Running the following PCA code could take several minutes or more, depending on your computer's processing power.</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set the random number generator for reproducability\n",
    "np.random.seed(49)\n",
    "\n",
    "# define the dimensionality of the hidden rep.\n",
    "n_components = 200\n",
    "\n",
    "# Randomly initialize the Weight matrix\n",
    "W = np.random.uniform(low=-4 * np.sqrt(6. / (n_components + mnist.shape[1])),\\\n",
    "                      high=4 * np.sqrt(6. / (n_components + mnist.shape[1])), size=(mnist.shape[1], n_components))\n",
    "# Initialize the step-size\n",
    "alpha = 1e-3\n",
    "# Initialize the gradient\n",
    "grad = np.infty\n",
    "# Set the tolerance \n",
    "tol = 1e-8\n",
    "# Initialize error\n",
    "old_error = 0\n",
    "error = [np.infty]\n",
    "batch_size = 250\n",
    "\n",
    "### train with stochastic gradients\n",
    "start_time = time.time()\n",
    "\n",
    "iter_idx = 1\n",
    "# loop until gradient updates become small\n",
    "while (alpha*np.linalg.norm(grad) > tol) and (iter_idx < 300):\n",
    "    for batch_idx in xrange(mnist.shape[0]/batch_size):\n",
    "        x = mnist[batch_idx*batch_size:(batch_idx+1)*batch_size, :]\n",
    "        h = np.dot(x, W)\n",
    "        x_recon = np.dot(h, W.T)\n",
    "        \n",
    "        # compute gradient\n",
    "        diff = x - x_recon\n",
    "        grad = (-4./batch_size)*np.dot(diff.T, h)\n",
    "    \n",
    "        # update parameters\n",
    "        W = W - alpha*grad\n",
    "    \n",
    "    # track the error\n",
    "    if iter_idx % 25 == 0:\n",
    "        old_error = error[-1]\n",
    "        diff = mnist - np.dot(np.dot(mnist, W), W.T)\n",
    "        recon_error = np.mean( np.sum(diff**2, 1) )\n",
    "        error.append(recon_error)\n",
    "        print \"Epoch %d, Reconstruction Error: %.3f\" %(iter_idx, recon_error)\n",
    "    \n",
    "    iter_idx += 1\n",
    "end_time = time.time()\n",
    "\n",
    "print\n",
    "print \"Training ended after %i iterations, taking a total of %.2f seconds.\" %(iter_idx, end_time-start_time)\n",
    "print \"Final Reconstruction Error: %.2f\" %(error[-1])\n",
    "reduced_mnist = np.dot(mnist, W)\n",
    "print \"Dataset is now of size: %d x %d\"%(reduced_mnist.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Let's visualize a reconstruction..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_idx = 2\n",
    "reconstructed_img = np.dot(reduced_mnist[img_idx,:], W.T)\n",
    "original_img = mnist[img_idx,:]\n",
    "\n",
    "# subplot for original image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.imshow(np.reshape(original_img, (28, 28)), cmap='Greys_r')\n",
    "ax1.set_title(\"Original Painting\")\n",
    "\n",
    "# subplot for reconstruction\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.imshow(np.reshape(reconstructed_img, (28, 28)), cmap='Greys_r')\n",
    "ax2.set_title(\"Reconstruction\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We can again visualize the transformation matrix $\\mathbf{W}^{T}$.  It's rows act as 'filters' or 'feature detectors'.  However, without the orthogonality constraint, we've loss the identifiably of the components..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# two components to show\n",
    "comp1 = 0\n",
    "comp2 = 150\n",
    "\n",
    "# subplot \n",
    "ax1 = plt.subplot(1,2,1)\n",
    "filter1 = W[:, comp1]\n",
    "ax1.imshow(np.reshape(filter1, (28, 28)), cmap='Greys_r')\n",
    "\n",
    "# subplot \n",
    "ax2 = plt.subplot(1,2,2)\n",
    "filter2 = W[:, comp2]\n",
    "ax2.imshow(np.reshape(filter2, (28, 28)), cmap='Greys_r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.  Nonlinear Dimensionality Reduction with Autoencoders\n",
    "\n",
    "In the last session (and section) we learned about *Principal Component Analysis*, a technique that finds some linear projection that reduces the dimensionality of the data while preserving its variance.  We looked at it as a form of unsupervised linear regression, where we predict the data itself instead of some associated value (i.e. a label).  In this section, we will move on to a *nonlinear* dimensionality reduction technique called an [Autoencoder](https://en.wikipedia.org/wiki/Autoencoder) and derive it's optimization procedure. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Defining the Autoencoder Model\n",
    "\n",
    "Recall that PCA is comprised of a linear projection step followed by application of the inverse projection.  An Autoencoder is the same model but with a non-linear transformation placed on the hidden representation.  To reiterate, our goal is: for a datapoint $\\mathbf{x}_{i}$, find a lower-dimensional representation $\\mathbf{h}_{i}$ such that $\\mathbf{x}_{i}$ can be 'predicted' from $\\mathbf{h}_{i}$---but this time, not necessarily with a linear transformation.  In math, this statement can be written as $$\\mathbf{\\tilde x}_{i} = \\mathbf{h}_{i} \\mathbf{W}^{T} \\text{ where } \\mathbf{h}_{i} = f(\\mathbf{x}_{i} \\mathbf{W}). $$  $\\mathbf{W}$ is a $D \\times K$ matrix of parameters that need to be learned--much like the $\\beta$ vector in regression models.  $D$ is the dimensionality of the original data, and $K$ is the dimensionality of the compressed representation $\\mathbf{h}_{i}$.  Lastly, we have the new component, the transformation function $f$.  There are many possible function to choose for $f$; yet we'll use a framilar one, the logistic function $$f(z) = \\frac{1}{1+\\exp(-z)}.$$  The graphic below depicts the autoencoder's computation path: \n",
    "\n",
    "![autoencoder_pipeline](./graphics/ae_pipeline.png)\n",
    "\n",
    "#### Optimization\n",
    "Having defined the Autoencoder model, we look to write learning as an optimization process.  Recall that we wish to make a reconstruction of the data, denoted $\\mathbf{\\tilde x}_{i}$, as close as possible to the original input: $$\\mathcal{L}_{\\text{AE}} = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{\\tilde x}_{i})^{2}.$$  We can make a substitution for $\\mathbf{\\tilde x}_{i}$ from the equation above: $$ = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{h}_{i}\\mathbf{W}^{T})^{2}.$$  And we can make another substitution for $\\mathbf{h}_{i}$, bringing us to the final form of the loss function: $$ = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - f(\\mathbf{x}_{i}\\mathbf{W})\\mathbf{W}^{T})^{2}.$$ \n",
    "\n",
    "\n",
    "## <span style=\"color:red\">STUDENT ACTIVITY (15 mins)</span>   \n",
    "Derive an expression for the gradient: $$ \\nabla_{W}\\mathcal{L}_{\\text{AE}} = ? $$ \n",
    "Take $f$ to be the logistic function, which has a derivative of $f'(z) = f(z)(1-f(z))$.  Those functions are provided for you below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def logistic(x):\n",
    "    return 1./(1+np.exp(-x))\n",
    "\n",
    "def logistic_derivative(x):\n",
    "    z = logistic(x)\n",
    "    return np.multiply(z, 1-z)\n",
    "\n",
    "def compute_gradient(x, x_recon, h, a):\n",
    "    # parameters:\n",
    "    # x: the original data\n",
    "    # x_recon: the reconstruction of x\n",
    "    # h: the hidden units (after application of f)\n",
    "    # a: the pre-activations (before the application of f)\n",
    "\n",
    "    return #TODO\n",
    "   \n",
    "np.random.seed(39)\n",
    "\n",
    "# dummy variables for testing\n",
    "x = np.random.normal(size=(5,3))\n",
    "x_recon = x + np.random.normal(size=x.shape)\n",
    "W = np.random.normal(size=(x.shape[1], 2))\n",
    "a = np.dot(x, W)\n",
    "h = logistic(a)\n",
    "compute_gradient(x, x_recon, h, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should print \n",
    "`\n",
    "array([[ 4.70101821,  2.26494039],\n",
    "       [ 2.86585042,  0.0731302 ],\n",
    "       [ 0.79869215,  0.15570277]])\n",
    "`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Autoencoder (AE) Overview\n",
    "\n",
    "_*Data*_\n",
    "\n",
    "We observe $\\mathbf{x}_{i}$ where\n",
    "\\begin{eqnarray*}\n",
    "\\mathbf{x}_{i} = (x_{i,1}, \\dots, x_{i,D}) &:& \\mbox{set of $D$ explanatory variables (aka features).  No labels.} \n",
    "\\end{eqnarray*}\n",
    "\n",
    "_* Parameters*_\n",
    "\n",
    "$\\mathbf{W}$: Matrix with dimensionality $D \\times K$, where $D$ is the dimensionality of the original data and $K$ the dimensionality of the new features. The matrix encodes the transformation between the original and new feature spaces.\n",
    "\n",
    "_*Error Function*_\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\mathcal{L} = \\sum_{i=1}^{N} ( \\mathbf{x}_{i} - f(\\mathbf{x}_{i} \\mathbf{W}) \\mathbf{W}^{T})^{2}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "### 2.2 Autoencoder Implementation\n",
    "\n",
    "Now let's train an Autoencoder..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set the random number generator for reproducability\n",
    "np.random.seed(39)\n",
    "\n",
    "# define the dimensionality of the hidden rep.\n",
    "n_components = 200\n",
    "\n",
    "# Randomly initialize the transformation matrix\n",
    "W = np.random.uniform(low=-4 * np.sqrt(6. / (n_components + mnist.shape[1])),\\\n",
    "                      high=4 * np.sqrt(6. / (n_components + mnist.shape[1])), size=(mnist.shape[1], n_components))\n",
    "\n",
    "# Initialize the step-size\n",
    "alpha = .01\n",
    "# Initialize the gradient\n",
    "grad = np.infty\n",
    "# Initialize error\n",
    "old_error = 0\n",
    "error = [np.infty]\n",
    "batch_size = 250\n",
    "\n",
    "### train with stochastic gradients\n",
    "start_time = time.time()\n",
    "\n",
    "iter_idx = 1\n",
    "# loop until gradient updates become small\n",
    "while (alpha*np.linalg.norm(grad) > tol) and (iter_idx < 300):\n",
    "    for batch_idx in xrange(mnist.shape[0]/batch_size):\n",
    "        x = mnist[batch_idx*batch_size:(batch_idx+1)*batch_size, :]\n",
    "        pre_act = np.dot(x, W) \n",
    "        h = logistic(pre_act)\n",
    "        x_recon = np.dot(h, W.T)\n",
    "        \n",
    "        # compute gradient\n",
    "        grad = compute_gradient(x, x_recon, h, pre_act)\n",
    "    \n",
    "        # update parameters\n",
    "        W = W - alpha/batch_size * grad\n",
    "    \n",
    "    # track the error\n",
    "    if iter_idx % 25 == 0:\n",
    "        old_error = error[-1]\n",
    "        \n",
    "        diff = mnist - np.dot(logistic(np.dot(mnist, W)), W.T)\n",
    "        recon_error = np.mean( np.sum(diff**2, 1) )\n",
    "        error.append(recon_error)\n",
    "        print \"Epoch %d, Reconstruction Error: %.3f\" %(iter_idx, recon_error)\n",
    "    \n",
    "    iter_idx += 1\n",
    "end_time = time.time()\n",
    "\n",
    "print\n",
    "print \"Training ended after %i iterations, taking a total of %.2f seconds.\" %(iter_idx, end_time-start_time)\n",
    "print \"Final Reconstruction Error: %.2f\" %(error[-1])\n",
    "reduced_mnist = np.dot(mnist, W)\n",
    "print \"Dataset is now of size: %d x %d\"%(reduced_mnist.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_idx = 2\n",
    "reconstructed_img = np.dot(logistic(reduced_mnist[img_idx,:]), W.T)\n",
    "original_img = mnist[img_idx,:]\n",
    "\n",
    "# subplot for original image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.imshow(np.reshape(original_img, (28, 28)), cmap='Greys_r')\n",
    "ax1.set_title(\"Original Digit\")\n",
    "\n",
    "# subplot for reconstruction\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.imshow(np.reshape(reconstructed_img, (28, 28)), cmap='Greys_r')\n",
    "ax2.set_title(\"Reconstruction\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# two components to show\n",
    "comp1 = 0\n",
    "comp2 = 150\n",
    "\n",
    "# subplot \n",
    "ax1 = plt.subplot(1,2,1)\n",
    "filter1 = W[:, comp1]\n",
    "ax1.imshow(np.reshape(filter1, (28, 28)), cmap='Greys_r')\n",
    "\n",
    "# subplot \n",
    "ax2 = plt.subplot(1,2,2)\n",
    "filter2 = W[:, comp2]\n",
    "ax2.imshow(np.reshape(filter2, (28, 28)), cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3  SciKit Learn Version\n",
    "We can hack the Scikit-Learn Regression neural network into an Autoencoder by feeding it the data back as the labels..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPRegressor\n",
    "\n",
    "# set the random number generator for reproducability\n",
    "np.random.seed(39)\n",
    "\n",
    "# define the dimensionality of the hidden rep.\n",
    "n_components = 200\n",
    "\n",
    "# define model\n",
    "ae = MLPRegressor(hidden_layer_sizes=(n_components,), activation='logistic')\n",
    "\n",
    "### train Autoencoder\n",
    "start_time = time.time()\n",
    "ae.fit(mnist, mnist)\n",
    "end_time = time.time()\n",
    "\n",
    "recon_error = np.mean(np.sum((mnist - ae.predict(mnist))**2, 1))\n",
    "W = ae.coefs_[0]\n",
    "b = ae.intercepts_[0]\n",
    "reduced_mnist = logistic(np.dot(mnist, W) + b)\n",
    "\n",
    "print\n",
    "print \"Training ended after a total of %.2f seconds.\" %(end_time-start_time)\n",
    "print \"Final Reconstruction Error: %.2f\" %(recon_error)\n",
    "print \"Dataset is now of size: %d x %d\"%(reduced_mnist.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_idx = 5\n",
    "reconstructed_img = np.dot(reduced_mnist[img_idx,:], ae.coefs_[1]) + ae.intercepts_[1]\n",
    "original_img = mnist[img_idx,:]\n",
    "\n",
    "# subplot for original image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.imshow(np.reshape(original_img, (28, 28)), cmap='Greys_r')\n",
    "ax1.set_title(\"Original Digit\")\n",
    "\n",
    "# subplot for reconstruction\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.imshow(np.reshape(reconstructed_img, (28, 28)), cmap='Greys_r')\n",
    "ax2.set_title(\"Reconstruction\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# two components to show\n",
    "comp1 = 0\n",
    "comp2 = 150\n",
    "\n",
    "# subplot \n",
    "ax1 = plt.subplot(1,2,1)\n",
    "filter1 = W[:, comp1]\n",
    "ax1.imshow(np.reshape(filter1, (28, 28)), cmap='Greys_r')\n",
    "\n",
    "# subplot \n",
    "ax2 = plt.subplot(1,2,2)\n",
    "filter2 = W[:, comp2]\n",
    "ax2.imshow(np.reshape(filter2, (28, 28)), cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4  Denoising Autoencoder (DAE)\n",
    "Lastly, we are going to examine an extension to the Autoencoder called a *Denoising Autoencoder* (DAE).  It has the following loss fuction: $$\\mathcal{L}_{\\text{DAE}} = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - f((\\hat{\\boldsymbol{\\zeta}} \\odot \\mathbf{x}_{i})\\mathbf{W})\\mathbf{W}^{T})^{2} \\ \\text{ where } \\hat{\\boldsymbol{\\zeta}} \\sim \\text{Bernoulli}(p).$$  In words, what we're doing is drawning a Bernoulli (i.e. binary) matrix the same size as the input features, and feeding a *corrupted* version of $\\mathbf{x}_{i}$.  The Autoencoder, then, must try to recreate the original data from a lossy representation.  This has the effect of forcing the Autoencoder to use features that better *generalize*. \n",
    "\n",
    "Let's make the simple change that implements a DAE below..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set the random number generator for reproducability\n",
    "np.random.seed(39)\n",
    "\n",
    "# define the dimensionality of the hidden rep.\n",
    "n_components = 200\n",
    "\n",
    "# Randomly initialize the Beta vector\n",
    "W = np.random.uniform(low=-4 * np.sqrt(6. / (n_components + mnist.shape[1])),\\\n",
    "                      high=4 * np.sqrt(6. / (n_components + mnist.shape[1])), size=(mnist.shape[1], n_components))\n",
    "\n",
    "# Initialize the step-size\n",
    "alpha = .01\n",
    "# Initialize the gradient\n",
    "grad = np.infty\n",
    "# Set the tolerance \n",
    "tol = 1e-8\n",
    "# Initialize error\n",
    "old_error = 0\n",
    "error = [np.infty]\n",
    "batch_size = 250\n",
    "\n",
    "### train with stochastic gradients\n",
    "start_time = time.time()\n",
    "\n",
    "iter_idx = 1\n",
    "# loop until gradient updates become small\n",
    "while (alpha*np.linalg.norm(grad) > tol) and (iter_idx < 300):\n",
    "    for batch_idx in xrange(mnist.shape[0]/batch_size):\n",
    "        x = mnist[batch_idx*batch_size:(batch_idx+1)*batch_size, :]\n",
    "        \n",
    "        # add noise to features\n",
    "        x_corrupt = np.multiply(x, np.random.binomial(n=1, p=.8, size=x.shape))\n",
    "        \n",
    "        pre_act = np.dot(x_corrupt, W) \n",
    "        h = logistic(pre_act)\n",
    "        x_recon = np.dot(h, W.T)\n",
    "        \n",
    "        # compute gradient\n",
    "        diff = x - x_recon\n",
    "        grad = -2.*(np.dot(diff.T, h) + np.dot(np.multiply(np.dot(diff, W), logistic_derivative(pre_act)).T, x_corrupt).T)\n",
    "        # NOTICE: during the 'backward pass', use the uncorrupted features\n",
    "    \n",
    "        # update parameters\n",
    "        W = W - alpha/batch_size * grad\n",
    "    \n",
    "    # track the error\n",
    "    if iter_idx % 25 == 0:\n",
    "        old_error = error[-1]\n",
    "        \n",
    "        diff = mnist - np.dot(logistic(np.dot(mnist, W)), W.T)\n",
    "        recon_error = np.mean( np.sum(diff**2, 1) )\n",
    "        error.append(recon_error)\n",
    "        print \"Epoch %d, Reconstruction Error: %.3f\" %(iter_idx, recon_error)\n",
    "    \n",
    "    iter_idx += 1\n",
    "end_time = time.time()\n",
    "\n",
    "print\n",
    "print \"Training ended after %i iterations, taking a total of %.2f seconds.\" %(iter_idx, end_time-start_time)\n",
    "print \"Final Reconstruction Error: %.2f\" %(error[-1])\n",
    "reduced_mnist = np.dot(mnist, W)\n",
    "print \"Dataset is now of size: %d x %d\"%(reduced_mnist.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_idx = 5\n",
    "reconstructed_img = np.dot(logistic(reduced_mnist[img_idx,:]), W.T)\n",
    "original_img = mnist[img_idx,:]\n",
    "\n",
    "# subplot for original image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.imshow(np.reshape(original_img, (28, 28)), cmap='Greys_r')\n",
    "ax1.set_title(\"Original Painting\")\n",
    "\n",
    "# subplot for reconstruction\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.imshow(np.reshape(reconstructed_img, (28, 28)), cmap='Greys_r')\n",
    "ax2.set_title(\"Reconstruction\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# two components to show\n",
    "comp1 = 0\n",
    "comp2 = 150\n",
    "\n",
    "# subplot \n",
    "ax1 = plt.subplot(1,2,1)\n",
    "filter1 = W[:, comp1]\n",
    "ax1.imshow(np.reshape(filter1, (28, 28)), cmap='Greys_r')\n",
    "\n",
    "# subplot \n",
    "ax2 = plt.subplot(1,2,2)\n",
    "filter2 = W[:, comp2]\n",
    "ax2.imshow(np.reshape(filter2, (28, 28)), cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When training larger autoencoders, you'll see filters that look like these...\n",
    "\n",
    "Regular Autoencoder:\n",
    "![AE Without Noise](http://deeplearning.net/tutorial/_images/filters_corruption_0.png)\n",
    "\n",
    "Denoising Autoencoder:\n",
    "![DAE](http://deeplearning.net/tutorial/_images/filters_corruption_30.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:red\">STUDENT ACTIVITY (until end of session)</span> \n",
    "Your task is to reproduce the faces experiment from the previous session but using an Autoencoder instead of PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_olivetti_faces\n",
    "\n",
    "faces_dataset = fetch_olivetti_faces(shuffle=True)\n",
    "faces = faces_dataset.data # 400 flattened 64x64 images\n",
    "person_ids = faces_dataset.target # denotes the identity of person (40 total)\n",
    "\n",
    "print \"Dataset size: %d x %d\" %(faces.shape)\n",
    "print \"And the images look like this...\"\n",
    "plt.imshow(np.reshape(faces[200,:], (64, 64)), cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset contains 400 64x64 pixel images of 40 people each exhibiting 10 facial expressions.  The images are in gray-scale, not color, and therefore flattened vectors contain 4096 dimensions.\n",
    "\n",
    "### <span style=\"color:red\">Subtask 1: Run (Regular) Autoencoder</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# train Autoencoder model on 'faces'\n",
    "\n",
    "###########################\n",
    "\n",
    "print \"Training took a total of %.2f seconds.\" %(end_time-start_time)\n",
    "print \"Final reconstruction error: %.2f%%\" %(recon_error) \n",
    "print \"Dataset is now of size: %d x %d\"%(faces_reduced.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 2: Reconstruct an image</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# Use learned transformation matrix to project back to the original 4096-dimensional space\n",
    "# Remember you need to use np.reshape() \n",
    "\n",
    "###########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 3: Train a Denoising Autoencoder</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "\n",
    "###########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 4: Generate a 2D scatter plot from both models</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# Run AE for 2 components\n",
    "\n",
    "# Generate plot\n",
    "\n",
    "# Bonus: color the scatter plot according to the person_ids to see if any structure can be seen\n",
    "\n",
    "###########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 5: Train a denoising version of PCA and test its performance</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# Run PCA but add noise to the input first\n",
    "\n",
    "###########################"
   ]
  }
 ],
 "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
}