pset2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem Set 2\n",
    "Designed by Kun He, with help from Kate Saenko.\n",
    "Programming part adapted from [Stanford CS231n](http://cs231n.stanford.edu/).\n",
    "\n",
    "## Preamble\n",
    "- **[Note] The provided code is written in Python 2.7 style, some parts are incompatible with 3.x. We're working on a fix right now.**\n",
    "- By now you should have a working Python installation. If not, refer to the instructions in problem set 1.\n",
    "- Clone or download this git repository, which contains the `cs591s2` codebase.\n",
    "- The code for loading the MNIST dataset requires the `scikit-learn` package. If you previously did the Anaconda installation to set up your Python environment (see problem set 1), you should already have it. Otherwise, you can install it following the instructions here: http://scikit-learn.org/stable/install.html\n",
    "- Be sure to show **all your steps** in derivations!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1: Cross-Entropy and Softmax\n",
    "This problem asks you to derive the cross-entropy loss for a multiclass classification problem using maximum likelihood.\n",
    "Consider the multiclass classification problem in which each input is assigned to one of $K$ mutually exclusive classes. The binary target variables $y_k$ ∈ {0, 1} have a \"one-hot\" coding scheme, where the value is 1 for the indicated class and 0 for all others. Assume that we can interpret the network outputs as $h_k(x,\\theta) = p(y_k = 1|x)$, or the probability of the kth class.\n",
    "\n",
    "**Q1.1**: Show that the maximum likelihood estimate of the parameters $\\theta$ can be obtained by minimizing the multiclass **cross-entropy** loss function \n",
    "<p>\n",
    "$L(\\theta)= - \\frac{1}{N}\\sum_{i=1}^{N} \\sum_{k=1}^{K} y_{ik} \\log(h_k(x_i,\\theta))$\n",
    "</p>\n",
    "<p>\n",
    "where $N$ is the number of examples $\\{x_i,y_i\\}$. </p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**[Double click here to add your answer]**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q1.2**: Assume the penultimate layer outputs raw scores (called logits) for each class $z_k(x,\\theta)$ for input $x$. Then the class probabilities can be obtained using **softmax** with **temperature parameter $T>0$**, i.e. \n",
    "<p>\n",
    "$h_k(x,\\theta)=\\frac{\\exp(z_k(x,\\theta)/T)}{\\sum_{j}exp(z_j(x,\\theta)/T)}$. \n",
    "</p> \n",
    "\n",
    "Consider a single pair $(x,y)$ where $y$ is a one-hot vector, and the corresponding loss $L=-\\sum_{k}y_k\\log h_k(x,\\theta)$. Compute the partial derivative $\\frac{\\partial L}{\\partial z_k(x,\\theta)}$. \n",
    "\n",
    "<p>\n",
    "Hint: with the standard softmax ($T=1$) it is known that\n",
    "$\\frac{\\partial L}{\\partial z_k(x,\\theta)}=h_k(x,\\theta)-y_{k}$.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**[Double click here to add your answer]**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2: Simple Regularization Methods\n",
    "In learning neural networks we typically minimize a loss function $\\mathcal{L}(\\theta)$ with respect to the network parameters $\\theta$. It is also important that we *regularize* the network to reduce overfitting. A simple and popular regularization strategy is to penalize some *norm* of $\\theta$.\n",
    "\n",
    "**Q2.1**:  L2 regularization\n",
    "\n",
    "We can penalize the L2 norm of $\\theta$: we modify our objective function to be $\\mathcal{L}(\\theta) + \\lambda \\|\\theta\\|^2$ where $\\lambda$ is the weight of regularization, and let $g=\\frac{\\partial \\mathcal{L}}{\\partial \\theta}$. Please derive the update rule for minimizing this objective using gradient descent with step size $\\eta$. In other words, at time $t+1$, express the new parameters $\\theta_{t+1}$ in terms of the old parameters $\\theta_t$, the gradient $g_t$, $\\eta$, and $\\lambda$.\n",
    "\n",
    "Hint: this update rule has an additional term compared to standard gradient descent, which is called \"weight decay\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**[Double click here to add your answer]**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q2.2**:  L1 regularization\n",
    "\n",
    "Now let's consider L1 regularization: our objective in this case is $\\mathcal{L}(\\theta) + \\lambda \\|\\theta\\|_1$. Derive the update rule. This leads to a different form of weight decay.\n",
    "\n",
    "(Technically this becomes Stochastic *Sub-Gradient* Descent since the L1 norm is not differentiable at 0. But practically it is usually not an issue.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**[Double click here to add your answer]**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3: Backprop in a simple MLP\n",
    "This problem asks you to derive all the steps of the backpropagation algorithm for a simple classification network. Consider a fully-connected neural network, also known as a multi-layer perceptron (MLP), with a single hidden layer and a one-node output layer. The hidden and output nodes use an elementwise sigmoid activation function and the loss layer uses cross-entropy loss:\n",
    "<p>\n",
    "$f(z)=\\frac{1}{1+exp(-z))}$\n",
    "<br>\n",
    "$L(\\hat{y},y)=-yln(\\hat{y}) - (1-y)ln(1-\\hat{y})$\n",
    "</p>\n",
    "<p>\n",
    "The computation graph for an example network is shown below. Note that it has an equal number of nodes in the input and hidden layer (3 each), but, in general, they need not be equal. Also, to make the application of backprop easier, we show the <i>computation graph</i> which shows the dot product and activation functions as their own nodes, rather than the usual graph showing a single node for both.\n",
    "</p>\n",
    "\n",
    "<img src=\"mlpgraph.png\" style=\"height:200px;\">\n",
    "\n",
    "The backpropagation algorithm for an MLP is described in **Goodfellow 6.5.4**, and is reproduced below. For simplicity, we will assume no regularization on the weights, so you can ignore the terms involving $\\Omega$. The forward step is: \n",
    "\n",
    "<img src=\"forward.png\" style=\"width:500px;\">\n",
    "\n",
    "and the backward step is:\n",
    "\n",
    "<img src=\"backward.png\" style=\"width:500px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write down each step of the backward pass explicitly for all layers, i.e. for the loss and $k=2,1$, compute all gradients above, expressing them as a function of variables $x, y, h, W, b$. <i>Hint: you should substitute the updated values for the gradient $g$ in each step and simplify as much as possible.</i>  Specifically, compute the following (we have replaced the superscript notation $u^{(i)}$ with $u^i$):\n",
    "\n",
    "**Q3.1**: $\\nabla_{\\hat{y}}L(\\hat{y},y)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.2**: $\\nabla_{a^2}J$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.3**: $\\nabla_{b^2}J$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.4**: $\\nabla_{W^2}J$ <br><i>Hint: this should be a vector, since $W^2$ is a vector. </i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.5**: $\\nabla_{h^1}J$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.6**: $\\nabla_{b^1}J$, $\\nabla_{W^1}J$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.7** Briefly, explain how would the computational speed of backpropagation be affected if it did not include a forward pass?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4: XOR problem\n",
    "Design a neural network (by hand) to solve the XOR problem. The network should take two binary variables and output 1 when only one of the two input variables is 1, and 0 otherwize. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q4.1** Write the XOR function in terms of the logical functions (gates) $OR(x_1,x_2)$, $AND(x_1,x_2)$, $NAND(x_1,x_2)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Design the network, using the hyperbolic tangent function (tanh, http://reference.wolfram.com/language/ref/Tanh.html) as the activation function in all of the nodes. In the diagram below, we have filled in most of the values of the parameters.\n",
    "Fill in the remaining parameters, keeping in mind that tanh outputs [-1,+1], not [0,1]. Note that we need to appropriately change the AND node to take +1/-1 as inputs. Also, we must add an extra last layer to convert the final output from +1/-1 to 0/1. \n",
    "\n",
    "Hint: assume that tanh outputs -1 for any input $x\\leq -2$, +1 for any input $x\\geq 2$, 0 for $x=0$.\n",
    "\n",
    "<img src=\"xor1.png\" style=\"height:130px;\"><img src=\"xor2.png\" style=\"height:110px;\">\n",
    "\n",
    "**Q4.2**: What are the missing weights $a,b,c,d,e$ of the OR, NAND, AND and CONVERT subnetworks, respectively?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** `[double click here to add a solution]` **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 5 (Programming): Implementing a simple MLP\n",
    "In this problem we will develop a neural network with fully-connected layers, aka Multi-Layer Perceptron (MLP). We will use it to perform classification on the MNIST dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from cs591s2.classifiers.mlp import TwoLayerMLP\n",
    "\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = (10.0, 8.0) # set default size of plots\n",
    "plt.rcParams['image.interpolation'] = 'nearest'\n",
    "plt.rcParams['image.cmap'] = 'gray'\n",
    "\n",
    "# for auto-reloading external modules\n",
    "# see http://stackoverflow.com/questions/1907993/autoreload-of-modules-in-ipython\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "def rel_error(x, y):\n",
    "    \"\"\" returns relative error \"\"\"\n",
    "    return np.max(np.abs(x - y) / (np.maximum(1e-8, np.abs(x) + np.abs(y))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the class `TwoLayerMLP` in the file `cs591s2/classifiers/mlp.py` to implement our network. The network parameters are stored in the instance variable `self.params` where keys are string parameter names and values are numpy arrays. Below, we initialize toy data and a toy model that we will use to develop your implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create a small net and some toy data to check your implementations.\n",
    "# Note that we set the random seed for repeatable experiments.\n",
    "input_size = 4\n",
    "hidden_size = 10\n",
    "num_classes = 3\n",
    "num_inputs = 5\n",
    "\n",
    "def init_toy_model(actv, std=1e-1):\n",
    "    np.random.seed(0)\n",
    "    return TwoLayerMLP(input_size, hidden_size, num_classes, std=std, activation=actv)\n",
    "\n",
    "def init_toy_data():\n",
    "    np.random.seed(1)\n",
    "    X = 10 * np.random.randn(num_inputs, input_size)\n",
    "    y = np.array([0, 1, 2, 2, 1])\n",
    "    return X, y\n",
    "\n",
    "X, y = init_toy_data()\n",
    "print 'X = ', X\n",
    "print\n",
    "print 'y = ', y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.1 Forward pass\n",
    "Our 2-layer MLP uses a softmax output layer (standard, $T=1$) and the multiclass cross-entropy loss to perform classification. Both are defined in Problem 1.\n",
    "\n",
    "Please take a look at method `TwoLayerMLP.loss` in the file `cs591s2/classifiers/mlp.py`. This function takes in the data and weight parameters, and computes the class scores (or the logits $z_k(x,\\theta)$), the loss ($L$), and the gradients on the parameters. \n",
    "\n",
    "- Implement the first part of the function to compute `scores` and `loss`. Afterwards, run the following two test cases.\n",
    "\n",
    "Note 1: **If you're not careful, you could run into numerical underflow/overflow problems with softmax and cross-entropy.** In particular, it involves the [log-sum-exp operation](https://en.wikipedia.org/wiki/LogSumExp) where exponentiated numbers are summed. This can result in underflow/overflow, e.g. getting \"nan\" (stands for Not A Number) for seemingly ordinary numerical operations. Read about the solution in the link.\n",
    "\n",
    "Note 2: You're strongly encouraged to implement in a **vectorized** way, by not using `for` loops over either the example index $i$ or class index $k$. `For` loops can be much slower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "net = init_toy_model('sigmoid')\n",
    "scores = net.loss(X)\n",
    "print '(1) Your scores:'\n",
    "print scores\n",
    "print\n",
    "correct_scores = np.asarray([[-0.20202237, -0.3826274,  -0.37823944],\n",
    "                             [-0.09999725, -0.29461579, -0.298211  ],\n",
    "                             [-0.04333782, -0.29038971, -0.3097813 ],\n",
    "                             [-0.0305597,  -0.26184001, -0.30413287], \n",
    "                             [-0.02637045, -0.24190783, -0.29750761]])\n",
    "\n",
    "# The difference should be very small (< 1e-7)\n",
    "print 'Difference between your scores and correct scores:'\n",
    "print np.sum(np.abs(scores - correct_scores))\n",
    "print \n",
    "\n",
    "loss, _ = net.loss(X, y, reg=0.1)\n",
    "correct_loss = 1.182248\n",
    "\n",
    "# The difference should be very small (< 1e-7)\n",
    "print '(2) Your loss: %f'%(loss)\n",
    "print 'Difference between your loss and correct loss:'\n",
    "print np.sum(np.abs(loss - correct_loss))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.2 Backward pass\n",
    "- Implement the second part to compute gradient of the loss with respect to the variables `W1`, `b1`, `W2`, and `b2`, stored in `grads`.\n",
    "\n",
    "Now debug your backward pass using a numeric gradient check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use numeric gradient checking to check your implementation of the backward pass.\n",
    "# If your implementation is correct, the difference between the numeric and\n",
    "# analytic gradients should be less than 1e-8 for each of W1, W2, b1, and b2.\n",
    "from cs591s2.utils.gradient_check import eval_numerical_gradient\n",
    "\n",
    "loss, grads = net.loss(X, y, reg=0.1)\n",
    "\n",
    "# these should all be very small\n",
    "for param_name in grads:\n",
    "    f = lambda W: net.loss(X, y, reg=0.1)[0]\n",
    "    param_grad_num = eval_numerical_gradient(f, net.params[param_name], verbose=False)\n",
    "    print '%s max relative error: %e' % (param_name, rel_error(param_grad_num, grads[param_name]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.3 Train the Sigmoid network\n",
    "To train the network we will use stochastic gradient descent (SGD), implemented in `TwoLayerNet.train`. Then we train a two-layer network on toy data.\n",
    "\n",
    "- Implement the prediction function `TwoLayerNet.predict`, which is called during training to keep track of training and validation accuracy.\n",
    "\n",
    "You should get the final training loss around 0.1, which is good, but not too great for such a toy problem.  One problem is that the gradient magnitude for W1 (the first layer weights) stays small all the time, and the neural net doesn't get much \"learning signals\". This has to do with the saturation problem of the sigmoid activation function, discussed in **Goodfellow 6.3.2**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "net = init_toy_model('sigmoid', std=1e-1)\n",
    "stats = net.train(X, y, X, y,\n",
    "                  learning_rate=0.5, reg=1e-5,\n",
    "                  num_epochs=100, verbose=False)\n",
    "\n",
    "print 'Final training loss: ', stats['loss_history'][-1]\n",
    "\n",
    "# plot the loss history and gradient magnitudes\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(stats['loss_history'])\n",
    "plt.xlabel('epoch')\n",
    "plt.ylabel('training loss')\n",
    "plt.title('Training Loss history')\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(stats['grad_magnitude_history'])\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('||W1||')\n",
    "plt.ylim(0,1)\n",
    "plt.title('Gradient magnitude history (W1)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.4 Using softplus activation\n",
    "In practice, the Rectified Linear Unit (ReLU) and softplus activations are much more successful than sigmoid. In particular, they don't suffer from the saturation problem. \n",
    "Now please look at method `TwoLayerMLP.loss`, where we have provided partial implementations for the forward and backward computation for the ReLU activation. \n",
    "\n",
    "- Complete the implementation for the **softplus** activation.\n",
    "- Train the network with softplus, and report your final training loss (see below).\n",
    "\n",
    "Note: an implementation of the softplus function for scalar inputs was given in problem 1, but here you need to implement a **numerically stable**, and **vectorized** version. \n",
    "\n",
    "We do the numerical checks as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "net = init_toy_model('softplus', std=1e-1)\n",
    "scores = net.loss(X)\n",
    "print '(1) Your scores:'\n",
    "print scores\n",
    "print\n",
    "correct_scores = np.asarray([[-0.77980255, -1.38043532, -0.8706494 ],\n",
    "                             [-0.19887433, -1.22677463, -0.59603118],\n",
    "                             [-0.39705994, -1.06544023, -0.91229474],\n",
    "                             [-0.09768135, -0.58642231, -0.63485595],\n",
    "                             [-0.03301159, -0.34776212, -0.4272356 ]])\n",
    "\n",
    "# The difference should be very small (< 1e-7)\n",
    "print 'Difference between your scores and correct scores:'\n",
    "print np.sum(np.abs(scores - correct_scores))\n",
    "print \n",
    "\n",
    "loss, grads = net.loss(X, y, reg=0.1)\n",
    "correct_loss = 1.320973\n",
    "\n",
    "# The difference should be very small (< 1e-7)\n",
    "print '(2) Your loss: %f'%(loss)\n",
    "print 'Difference between your loss and correct loss:'\n",
    "print np.sum(np.abs(loss - correct_loss))\n",
    "print\n",
    "    \n",
    "# The differences should all be very small\n",
    "print '(3) checking gradients'\n",
    "for param_name in grads:\n",
    "    f = lambda W: net.loss(X, y, reg=0.1)[0]\n",
    "    param_grad_num = eval_numerical_gradient(f, net.params[param_name], verbose=False)\n",
    "    print '%s max relative error: %e' % (param_name, rel_error(param_grad_num, grads[param_name]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that it's working, let's train the network. Is the neural net able to get stronger learning signals (i.e. gradients) this time? Report your final training loss; it should be better than what you get with sigmoid. \n",
    "\n",
    "You're welcome to try ReLU as well. Should the results be similar to using softplus, or not? Are they actually similar?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "net = init_toy_model('softplus', std=1e-1)\n",
    "stats = net.train(X, y, X, y,\n",
    "                  learning_rate=0.1, reg=1e-5,\n",
    "                  num_epochs=50, verbose=False)\n",
    "\n",
    "print 'Final training loss: ', stats['loss_history'][-1]\n",
    "\n",
    "# plot the loss history\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(stats['loss_history'])\n",
    "plt.xlabel('epoch')\n",
    "plt.ylabel('training loss')\n",
    "plt.title('Training Loss history')\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(stats['grad_magnitude_history'])\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('||W1||')\n",
    "plt.title('Gradient magnitude history (W1)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load MNIST data\n",
    "Now that you have implemented a two-layer network that works on toy data, let's try some real data. The MNIST dataset is a standard benchmark in machine learning. It consists of 70,000 grayscale handwritten digit images, which we split into 50,000 training, 10,000 validation and 10,000 testing. The images are of size 28x28, which are flattened into 784-d vectors.\n",
    "\n",
    "**Note**: the function `get_MNIST_data` requires the `scikit-learn` package. If you previously did anaconda installation to set up your Python environment, you should already have it. Otherwise, you can install it following the instructions here: http://scikit-learn.org/stable/install.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load MNIST\n",
    "from cs591s2.utils.data_utils import get_MNIST_data\n",
    "X_train, y_train, X_val, y_val, X_test, y_test = get_MNIST_data()\n",
    "print 'Train data shape: ', X_train.shape\n",
    "print 'Train labels shape: ', y_train.shape\n",
    "print 'Validation data shape: ', X_val.shape\n",
    "print 'Validation labels shape: ', y_val.shape\n",
    "print 'Test data shape: ', X_test.shape\n",
    "print 'Test labels shape: ', y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.5 Train a network on MNIST\n",
    "We will now train a network on MNIST with 64 hidden units in the hidden layer. We train it using SGD, and decrease the learning rate with an exponential rate over time; this is achieved by multiplying the learning rate with a constant factor `learning_rate_decay` (which is less than 1) after each epoch. In effect, we are using a high learning rate initially, which is good for exploring the solution space, and using lower learning rates later to encourage convergence to a local minimum (or [saddle point](http://www.offconvex.org/2016/03/22/saddlepoints/), which may happen more often).\n",
    "\n",
    "- Train your MNIST network with 3 different activation functions: sigmoid, softplus, and ReLU. Compare their performances.\n",
    "- Which activation function would you choose in practice? Why?\n",
    "\n",
    "We first define some variables and utility functions. The `plot_stats` function plots the hisotories of gradient magnitude, training loss, and accuracies on the training and validation sets. The `visualize_weights` function visualizes the weights learned in the first layer of the network. In most neural networks trained on visual data, the first layer weights typically show some visible structure when visualized. Both functions help you to diagnose the training process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "input_size = 28 * 28\n",
    "hidden_size = 64\n",
    "num_classes = 10\n",
    "\n",
    "# Plot the loss function and train / validation accuracies\n",
    "def plot_stats(stats):\n",
    "    plt.subplot(3, 1, 1)\n",
    "    plt.plot(stats['grad_magnitude_history'])\n",
    "    plt.title('Gradient magnitude history (W1)')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('||W1||')\n",
    "    plt.ylim(0, np.minimum(100,np.max(stats['grad_magnitude_history'])))\n",
    "    plt.subplot(3, 1, 2)\n",
    "    plt.plot(stats['loss_history'])\n",
    "    plt.title('Loss history')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Loss')\n",
    "    plt.ylim(0, 100)\n",
    "    plt.subplot(3, 1, 3)\n",
    "    plt.plot(stats['train_acc_history'], label='train') \n",
    "    plt.plot(stats['val_acc_history'], label='val')\n",
    "    plt.title('Classification accuracy history')\n",
    "    plt.xlabel('Epoch')\n",
    "    plt.ylabel('Clasification accuracy')\n",
    "    plt.show()\n",
    "\n",
    "# Visualize the weights of the network\n",
    "from cs591s2.utils.vis_utils import visualize_grid\n",
    "def show_net_weights(net):\n",
    "    W1 = net.params['W1']\n",
    "    W1 = W1.T.reshape(-1, 28, 28)\n",
    "    plt.imshow(visualize_grid(W1, padding=3).astype('uint8'))\n",
    "    plt.gca().axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.5.1 Sigmoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sigmoid_net = TwoLayerMLP(input_size, hidden_size, num_classes, activation='sigmoid', std=1e-1)\n",
    "\n",
    "# Train the network\n",
    "sigmoid_stats = sigmoid_net.train(X_train, y_train, X_val, y_val,\n",
    "                                  num_epochs=20, batch_size=100,\n",
    "                                  learning_rate=1e-3, \n",
    "                                  learning_rate_decay=0.95,\n",
    "                                  reg=0.5, verbose=True)\n",
    "print\n",
    "\n",
    "# Predict on the training set\n",
    "train_acc = (sigmoid_net.predict(X_train) == y_train).mean()\n",
    "print 'Sigmoid final training accuracy: ', train_acc\n",
    "\n",
    "# Predict on the validation set\n",
    "val_acc = (sigmoid_net.predict(X_val) == y_val).mean()\n",
    "print 'Sigmoid final validation accuracy: ', val_acc\n",
    "\n",
    "# Predict on the test set\n",
    "test_acc = (sigmoid_net.predict(X_test) == y_test).mean()\n",
    "print 'Sigmoid test accuracy: ', test_acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show stats and visualizations\n",
    "plot_stats(sigmoid_stats)\n",
    "show_net_weights(sigmoid_net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.5.2 Softplus\n",
    "\n",
    "Note: **Be careful with potential numerical overflow in your softplus implementation!** The toy example may not have exposed that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "softplus_net = TwoLayerMLP(input_size, hidden_size, num_classes, activation='softplus', std=1e-1)\n",
    "\n",
    "# Train the network\n",
    "softplus_stats = softplus_net.train(X_train, y_train, X_val, y_val,\n",
    "                                    num_epochs=20, batch_size=100, \n",
    "                                    learning_rate=1e-3, \n",
    "                                    learning_rate_decay=0.95,\n",
    "                                    reg=0.5, verbose=True)\n",
    "print\n",
    "\n",
    "# Predict on the training set\n",
    "train_acc = (softplus_net.predict(X_train) == y_train).mean()\n",
    "print 'Softplus final training accuracy: ', train_acc\n",
    "\n",
    "# Predict on the validation set\n",
    "val_acc = (softplus_net.predict(X_val) == y_val).mean()\n",
    "print 'Softplus final validation accuracy: ', val_acc\n",
    "\n",
    "# Predict on the test set\n",
    "test_acc = (softplus_net.predict(X_test) == y_test).mean()\n",
    "print 'Softplus test accuracy: ', test_acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show stats and visualizations\n",
    "plot_stats(softplus_stats)\n",
    "show_net_weights(softplus_net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Q5.5.3 ReLU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "relu_net = TwoLayerMLP(input_size, hidden_size, num_classes, activation='relu', std=1e-1)\n",
    "\n",
    "# Train the network\n",
    "relu_stats = relu_net.train(X_train, y_train, X_val, y_val,\n",
    "                            num_epochs=20, batch_size=100,\n",
    "                            learning_rate=1e-3, \n",
    "                            learning_rate_decay=0.95,\n",
    "                            reg=0.5, verbose=True)\n",
    "print\n",
    "\n",
    "# Predict on the training set\n",
    "train_acc = (relu_net.predict(X_train) == y_train).mean()\n",
    "print 'ReLU final training accuracy: ', train_acc\n",
    "\n",
    "# Predict on the validation set\n",
    "val_acc = (relu_net.predict(X_val) == y_val).mean()\n",
    "print 'ReLU final validation accuracy: ', val_acc\n",
    "\n",
    "# Predict on the test set\n",
    "test_acc = (relu_net.predict(X_test) == y_test).mean()\n",
    "print 'ReLU test accuracy: ', test_acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show stats and visualizations\n",
    "plot_stats(relu_stats)\n",
    "show_net_weights(relu_net)"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}