P03-C03-mlp-dropout-scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dropout regularization from scratch\n",
    "\n",
    "[**chapter status - very rough draft**]\n",
    "\n",
    "If you're reading the tutorials in sequence, \n",
    "then you might remember from Part 2 \n",
    "that machine learning models \n",
    "can be susceptible to overfitting. \n",
    "To recap: in machine learning,\n",
    "our goal is to discover general patterns.\n",
    "For example, we might want to learn as association between genetic markers\n",
    "and the development of dementia in adulthood. \n",
    "Our hope would be to uncover a pattern that could be applied successfully to assess risk for the entire population.\n",
    "\n",
    "However, when we train models, we don't have access to the entire population (or current or potential humans).\n",
    "Instead, we can access only a small, finite sample.\n",
    "Even in a large hospital system, we might get hundreds of thousands of medical records. \n",
    "Given such a finite sample size, it's possible to uncover spurious associations \n",
    "that don't hold up for unseen data.\n",
    "\n",
    "Let's consider an extreme pathological case. \n",
    "Imagine that you want to learn to predict\n",
    "which people will repay their loans. \n",
    "A lender hires you as a data scientist \n",
    "to investigate the case and gives you complete files on 100 applicants,\n",
    "of which 5 defaulted on their loans within 3 years. \n",
    "The files might include hundreds of features \n",
    "including income, occupation, credit score, length of employment etcetera.\n",
    "Imagine that they additionally give you video footage of their interview with a lending agent. \n",
    "That might seem like a lot of data! \n",
    "\n",
    "Now suppose that after generating an enormous set of features,\n",
    "you discover that of the 5 applicants who defaults, \n",
    "all 5 were wearing blue shirts during their interviews,\n",
    "while only 40% of general population wore blue shirts. \n",
    "Theirs a good chance that any model you train would pick up on this signal \n",
    "and use it as an important part of its learned pattern.\n",
    "\n",
    "Even if defaulters are no more lkely to wear blue shirts, \n",
    "there's a 1% chance that we'll observe all five defaulters wearing blue shirts.\n",
    "And keeping the sample size low we have hundreds or thousands of features,\n",
    "we may observe a large number of spurious correlations. \n",
    "Given trillions of training examples, these false associations might disappear. \n",
    "But we seldom have that luxury.\n",
    "\n",
    "The phenomena of fitting our training distribution more closely than the real distribution\n",
    "is called *overfitting*, and the techniques used to combat overfitting are called regularization.\n",
    "In the previous chapter, we introduced one classical approach to regularize statistical models. \n",
    "We penalized the size (the $\\ell^2_2$ nor) of the weights, coercing them to take smaller values.\n",
    "In probabilistic terms we might say this imposes a Gaussian prior on the value of the weights. \n",
    "But in more intuitive, functional terms, we can say this encourages the model to spead out its weights among many features and not to depend to much on a small number of potentially spurious associations. \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With great flexibility comes overfitting liability\n",
    "\n",
    "Given many more features than examples, linear models can overfit. \n",
    "But when there are many more examples than features, \n",
    "linear models can usually be counted on not to overfit.\n",
    "Unfortunately this propensity to generalize well comes at a cost. \n",
    "For every feature, a linear model has to assign it either positive or negative weight.\n",
    "Linear models can't take into account nuanced interactions between features.\n",
    "In more formal texts, you'll see this phenomena discussed as the bias-variance tradeoff.\n",
    "Linear models have high bias, (they can only represent a small class of functions),\n",
    "but low variance (they give similar results across different random samples of the data).\n",
    "[**point to more formal discussion of generalziation when chapter exists**]\n",
    "\n",
    "Deep neural networks, however, occupy the opposite end of the bias-variance spectrum.\n",
    "Neural networks are so flexible because they aren't confined to looking at each feature individually.\n",
    "Instead, they can learn complex interactions among groups of features. \n",
    "For example, they might infer that \"Nigeria\" and \"Wester Union\" \n",
    "appearing together in an email indicates spam \n",
    "but that \"Nigeria\" without \"Western Union\" does not connote spam. \n",
    "\n",
    "Even for a small number of features, deep neural networks are capable of overfitting.\n",
    "As one demonstration of the incredible flexibility of neural networks,\n",
    "researchers showed that [neural networks perfectly classify randomly labeled data](https://arxiv.org/abs/1611.03530).\n",
    "Let's think about what means. \n",
    "If the labels are assigned uniformly at random, and there are 10 classes, \n",
    "then no classifier can get better than 10% accuracy on holdout data. \n",
    "Yet even in these situations, when there is no true pattern to be learned, \n",
    "neural networks can perfectly fit the training labels. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dropping out activations\n",
    "\n",
    "In 2012, Professor Geoffrey Hinton and his students including Nitish Srivastava \n",
    "introduced a new idea for how to regularize neural network models. \n",
    "The inuition goes something like this. \n",
    "When a neural network overfits badly to training data,\n",
    "each layer depends too heavily on the exact configuration\n",
    "of features in the previous layer. \n",
    "\n",
    "To prevent the neural network from depending too much on any exact activation pathway,\n",
    "Hinton and Srivastava proposed randomly *dropping out* (i.e. setting to $0$) \n",
    "the hidden nodes in every layer with probability $.5$.\n",
    "Given a network with $n$ nodes we are sampling uniformly at random from the $2^n$ \n",
    "networks in which a subset of the nodes are turned off. \n",
    "\n",
    "![](./img/dropout.png)\n",
    "\n",
    "One intuition here is that because the nodes to drop out are chosen randomly on every pass,\n",
    "the representations in each layer can't depend on the exact values taken by nodes in the previous layer. \n",
    "\n",
    "## Making predictions with dropout models\n",
    "\n",
    "However, when it comes time to make predictions, \n",
    "we want to use the full representational power of our model. \n",
    "In other words, we don't want to drop out activations at test time.\n",
    "One principled way to justify the use of all nodes simultaneously,\n",
    "despite not training in this fasion,\n",
    "is that it's a form of model averaging. \n",
    "At each layer we average the representations of all of the $2^n$ dropout networks.\n",
    "Because each node has a $.5$ probability of being on during training, \n",
    "its vote is scaled by $.5$ when we use all nodes at prediction time\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import mxnet as mx\n",
    "import numpy as np\n",
    "from mxnet import nd, autograd\n",
    "import numpy as np\n",
    "mx.random.seed(1)\n",
    "ctx = mx.gpu()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The MNIST dataset\n",
    "\n",
    "Let's go ahead and grab our data.\n",
    "\n",
    "[**SWITCH TO CIFAR TO GET BETTER FEEL FOR REGULARIZATION**]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mnist = mx.test_utils.get_mnist()\n",
    "batch_size = 64\n",
    "def transform(data, label):\n",
    "    return data.astype(np.float32)/255, label.astype(np.float32)\n",
    "train_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=True, transform=transform),\n",
    "                                      batch_size, shuffle=True)\n",
    "test_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=False, transform=transform),\n",
    "                                     batch_size, shuffle=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W1 = nd.random_normal(shape=(784,256), ctx=ctx) *.01\n",
    "b1 = nd.random_normal(shape=256, ctx=ctx) * .01\n",
    "\n",
    "W2 = nd.random_normal(shape=(256,128), ctx=ctx) *.01\n",
    "b2 = nd.random_normal(shape=128, ctx=ctx) * .01\n",
    "\n",
    "W3 = nd.random_normal(shape=(128,10), ctx=ctx) *.01\n",
    "b3 = nd.random_normal(shape=10, ctx=ctx) *.01\n",
    "\n",
    "params = [W1, b1, W2, b2, W3, b3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, let's allocate space for gradients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for param in params:\n",
    "    param.attach_grad()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activation functions\n",
    "\n",
    "If we compose a multi-layer network but use only linear operations, then our entire network will still be a linear function. That's because $\\hat{y} = X \\cdot W_1 \\cdot W_2 \\cdot W_2 = X \\cdot W_4 $ for $W_4 = W_1 \\cdot W_2 \\cdot W3$. To give our model the capacity to capture nonlinear functions, we'll need to interleave our linear operations with activation functions. In this case, we'll use the rectified linear unit (ReLU):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def relu(X):\n",
    "    return nd.maximum(X, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dropout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dropout(X, drop_probability):\n",
    "    keep_probability = 1 - drop_probability\n",
    "    mask = nd.random_uniform(0, 1.0, X.shape, ctx=X.context) < keep_probability\n",
    "    #############################\n",
    "    #  Avoid division by 0 when scaling\n",
    "    #############################\n",
    "    if keep_probability > 0.0:\n",
    "        scale = (1/keep_probability)\n",
    "    else:\n",
    "        scale = 0.0\n",
    "    return mask * X * scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "[[  0.   1.   2.   3.]\n",
       " [  4.   5.   6.   7.]\n",
       " [  8.   9.  10.  11.]\n",
       " [ 12.  13.  14.  15.]\n",
       " [ 16.  17.  18.  19.]]\n",
       "<NDArray 5x4 @cpu(0)>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = nd.arange(20).reshape((5,4))\n",
    "dropout(A, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "[[  0.   0.   4.   6.]\n",
       " [  8.  10.  12.  14.]\n",
       " [  0.  18.  20.  22.]\n",
       " [ 24.   0.   0.  30.]\n",
       " [ 32.   0.   0.  38.]]\n",
       "<NDArray 5x4 @cpu(0)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dropout(A, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "[[ 0.  0.  0.  0.]\n",
       " [ 0.  0.  0.  0.]\n",
       " [ 0.  0.  0.  0.]\n",
       " [ 0.  0.  0.  0.]\n",
       " [ 0.  0.  0.  0.]]\n",
       "<NDArray 5x4 @cpu(0)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dropout(A, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Softmax output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def softmax(y_linear):\n",
    "    exp = nd.exp(y_linear-nd.max(y_linear))\n",
    "    partition = nd.nansum(exp, axis=0, exclude=True).reshape((-1,1))\n",
    "    return exp / partition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The *softmax* cross-entropy loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def softmax_cross_entropy(yhat_linear, y):\n",
    "    return - nd.nansum(y * nd.log_softmax(yhat_linear), axis=0, exclude=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the model\n",
    "\n",
    "Now we're ready to define our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def net(X, drop_prob=0.0):\n",
    "    #######################\n",
    "    #  Compute the first hidden layer \n",
    "    #######################    \n",
    "    h1_linear = nd.dot(X, W1) + b1\n",
    "    h1 = relu(h1_linear)\n",
    "    h1 = dropout(h1, drop_prob)\n",
    "    \n",
    "    #######################\n",
    "    #  Compute the second hidden layer\n",
    "    #######################\n",
    "    h2_linear = nd.dot(h1, W2) + b2\n",
    "    h2 = relu(h2_linear)\n",
    "    h2 = dropout(h2, drop_prob)\n",
    "    \n",
    "    #######################\n",
    "    #  Compute the output layer.\n",
    "    #  We will omit the softmax function here \n",
    "    #  because it will be applied \n",
    "    #  in the softmax_cross_entropy loss\n",
    "    #######################\n",
    "    yhat_linear = nd.dot(h2, W3) + b3\n",
    "    return yhat_linear"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def SGD(params, lr):    \n",
    "    for param in params:\n",
    "        param[:] = param - lr * param.grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation metric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def evaluate_accuracy(data_iterator, net):\n",
    "    numerator = 0.\n",
    "    denominator = 0.\n",
    "    for i, (data, label) in enumerate(data_iterator):\n",
    "        data = data.as_in_context(ctx).reshape((-1,784))\n",
    "        label = label.as_in_context(ctx)\n",
    "        label_one_hot = nd.one_hot(label, 10)\n",
    "        output = net(data)\n",
    "        predictions = nd.argmax(output, axis=1)\n",
    "        numerator += nd.sum(predictions == label)\n",
    "        denominator += data.shape[0]\n",
    "    return (numerator / denominator).asscalar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Execute the training loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 0. Loss: 0.737156236043, Train_acc 0.850967, Test_acc 0.8539\n",
      "Epoch 1. Loss: 0.394209427167, Train_acc 0.92225, Test_acc 0.923\n",
      "Epoch 2. Loss: 0.296510421107, Train_acc 0.946, Test_acc 0.9452\n",
      "Epoch 3. Loss: 0.232048722573, Train_acc 0.9563, Test_acc 0.9531\n",
      "Epoch 4. Loss: 0.205553142149, Train_acc 0.962967, Test_acc 0.9591\n",
      "Epoch 5. Loss: 0.178349442085, Train_acc 0.969233, Test_acc 0.9641\n",
      "Epoch 6. Loss: 0.175119599567, Train_acc 0.9735, Test_acc 0.967\n",
      "Epoch 7. Loss: 0.157515936016, Train_acc 0.975067, Test_acc 0.9688\n",
      "Epoch 8. Loss: 0.14164880119, Train_acc 0.977933, Test_acc 0.9705\n",
      "Epoch 9. Loss: 0.129475182254, Train_acc 0.9798, Test_acc 0.9729\n"
     ]
    }
   ],
   "source": [
    "epochs = 10\n",
    "moving_loss = 0.\n",
    "learning_rate = .001\n",
    "\n",
    "for e in range(epochs):\n",
    "    for i, (data, label) in enumerate(train_data):\n",
    "        data = data.as_in_context(ctx).reshape((-1,784))\n",
    "        label = label.as_in_context(ctx)\n",
    "        label_one_hot = nd.one_hot(label, 10)\n",
    "        with autograd.record():\n",
    "            ################################\n",
    "            #   Drop out 50% of hidden activations on the forward pass\n",
    "            ################################\n",
    "            output = net(data, drop_prob=.5)\n",
    "            loss = softmax_cross_entropy(output, label_one_hot)\n",
    "        loss.backward()\n",
    "        SGD(params, learning_rate)\n",
    "        \n",
    "        ##########################\n",
    "        #  Keep a moving average of the losses\n",
    "        ##########################\n",
    "        if i == 0:\n",
    "            moving_loss = nd.mean(loss).asscalar()\n",
    "        else:\n",
    "            moving_loss = .99 * moving_loss + .01 * nd.mean(loss).asscalar()\n",
    "            \n",
    "    test_accuracy = evaluate_accuracy(test_data, net)\n",
    "    train_accuracy = evaluate_accuracy(train_data, net)\n",
    "    print(\"Epoch %s. Loss: %s, Train_acc %s, Test_acc %s\" % (e, moving_loss, train_accuracy, test_accuracy)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "Noice. With just two hidden layers containing 256 and 128 hidden nodes, repsectively, we can achieve over 95% accuracy on this task. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For whinges or inquiries, [open an issue on  GitHub.](https://github.com/zackchase/mxnet-the-straight-dope)"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}