P02-C05-regularization-scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overfitting and regularization\n",
    "\n",
    "In [the last tutorial](./P02-C03-softmax-regression-scratch.ipynb), we introduced the task of multiclass classification. We showed how you can tackle this problem with a linear model called logistic regression. Owing to some amount of randomness, you might get slightly different results, but when I ran the notebook, the model achieved 88.1% accuracy on the training data and actually did slightly (but not significantly) better on the test data than on the traning data. \n",
    " \n",
    "Not every algorithm that performs well on training data will also perform test data. Take, for example, a trivial algorithm that memorizes its inputs and stores the associated labels. This model would have 100% accuracy on training data but would have no way of making any prediction at all on previously unseen data. \n",
    "\n",
    "The goal of supervised learning is to produce models that *generalize* to previously unseen data. When a model achieves low error on training data but performs much worse on test data, we say that the model has *overfit*. This means that the model has caught on to idiosyncratic features of the training data (e.g. one \"2\" happened to have a white pixel in top-right corner), but hasn't really picked up on general patterns. \n",
    "\n",
    "We can express this more formally. The quantity we really care about is the test error $e$. Because this quantity reflects the error of our model when generalized previously unseen data, we commonly call it the *generalization error*. When we have simple models and abundant data, we expect the generalization error to resemble the training error. When we work with more complex models and fewer examples, we expect the training error to go down but the generalization gap to grow. Fixing the size of the dataset, the following graph should give you some intuition about what we generally expect to see.\n",
    "\n",
    "![](img/regularization-overfitting.png)\n",
    "\n",
    "What precisely constitutes model complexity is a complex matter. Many factors govern whether a model will generalize well. For example a model with more parameters might be considered more complex. A model whose parameters can take a wider range of values might be more comples. Often with neural networks, we think of a model that takes more training steps as more complex, and one subject to *early stopping* as less complex. \n",
    "\n",
    "It can be dificult to compare the complexity among members of very different model classes (say decision trees versus neural networks). Researchers in the field of statistical learning theory have developed a large body of mathematical analysis that formulizes the notion of model complexity and provides guarantees on the generalization error for simple classes of models. *We won't get into this theory but may delve deeper in a future chapter*.\n",
    "\n",
    "To give you some intuition in this chapter, we'll focus on a few factors tend to influence the generalizability of a model class:\n",
    "1. **The number of tunable parameters.** When the number of tunable parameters, sometimes denoted the number of degress of freedom, is large, models tend to be more susceptible to overfitting.\n",
    "2. **The values taken by the parameters.** When weights can take a wider range of values, models can be more susceptible to over fitting.\n",
    "3. **The number of training examples.** It's trivially easy to overfit a dataset containing only one or two examples even if your model is simple. But overfitting a dataset with millions of examples requires an extremely flexible model.\n",
    "\n",
    "\n",
    "When classified handwritten digits before, we didn't overfit because our 60,000 training examples far out numbered the  $784 \\times 10 = 7,840$ weights plus $10$ bias terms gave us far fewer parameters than training examples. Let's see how things can go wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import mxnet as mx\n",
    "from mxnet import nd, autograd\n",
    "import numpy as np\n",
    "ctx = mx.cpu()\n",
    "mx.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the MNIST dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mnist = mx.test_utils.get_mnist()\n",
    "num_examples = 1000\n",
    "batch_size = 64\n",
    "train_data = mx.io.NDArrayIter(\n",
    "    mnist[\"train_data\"][:num_examples], \n",
    "    mnist[\"train_label\"][:num_examples], \n",
    "    batch_size, shuffle=True)\n",
    "test_data = mx.io.NDArrayIter(\n",
    "    mnist[\"test_data\"][:num_examples], \n",
    "    mnist[\"test_label\"][:num_examples], \n",
    "    batch_size, shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Allocate model parameters and define model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = nd.random_normal(shape=(784,10))\n",
    "b = nd.random_normal(shape=10)\n",
    "\n",
    "params = [W, b]\n",
    "\n",
    "for param in params:\n",
    "    param.attach_grad()\n",
    "    \n",
    "def net(X):\n",
    "    y_linear = nd.dot(X, W) + b\n",
    "    yhat = nd.softmax(y_linear, axis=1)\n",
    "    return yhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define loss function and optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cross_entropy(yhat, y):\n",
    "    return - nd.sum(y * nd.log(yhat), axis=0, exclude=True)\n",
    "\n",
    "def SGD(params, lr):    \n",
    "    for param in params:\n",
    "        param[:] = param - lr * param.grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write evaluation loop to calculate accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def evaluate_accuracy(data_iterator, net):\n",
    "    numerator = 0.\n",
    "    denominator = 0.\n",
    "    data_iterator.reset()\n",
    "    for i, batch in enumerate(data_iterator):\n",
    "        data = batch.data[0].as_in_context(ctx).reshape((-1,784))\n",
    "        label = batch.label[0].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 training loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed epoch 100. Loss: 0.379656352769, Train_acc 0.885742, Test_acc 0.689453\n",
      "Completed epoch 200. Loss: 0.152990181215, Train_acc 0.938477, Test_acc 0.728516\n",
      "Completed epoch 300. Loss: 0.0901987688058, Train_acc 0.973633, Test_acc 0.728516\n",
      "Completed epoch 400. Loss: 0.0613622239244, Train_acc 0.986328, Test_acc 0.738281\n",
      "Completed epoch 500. Loss: 0.0450707157253, Train_acc 0.995117, Test_acc 0.734375\n",
      "Completed epoch 600. Loss: 0.0343953531825, Train_acc 0.999023, Test_acc 0.737305\n",
      "Completed epoch 700. Loss: 0.0272877791812, Train_acc 1.0, Test_acc 0.74707\n",
      "Completed epoch 800. Loss: 0.0225696076896, Train_acc 1.0, Test_acc 0.748047\n",
      "Completed epoch 900. Loss: 0.0192749091022, Train_acc 1.0, Test_acc 0.749023\n",
      "Completed epoch 1000. Loss: 0.0168556058273, Train_acc 1.0, Test_acc 0.751953\n"
     ]
    }
   ],
   "source": [
    "epochs = 1000\n",
    "moving_loss = 0.\n",
    "\n",
    "for e in range(epochs):\n",
    "    train_data.reset()\n",
    "    for i, batch in enumerate(train_data):\n",
    "        data = batch.data[0].as_in_context(ctx).reshape((-1,784))\n",
    "        label = batch.label[0].as_in_context(ctx)\n",
    "        label_one_hot = nd.one_hot(label, 10)\n",
    "        with autograd.record():\n",
    "            output = net(data)\n",
    "            loss = cross_entropy(output, label_one_hot)\n",
    "        loss.backward()\n",
    "        SGD(params, .001)\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",
    "    if e % 100 == 99:\n",
    "        print(\"Completed epoch %s. Loss: %s, Train_acc %s, Test_acc %s\" % \n",
    "              (e+1, moving_loss, train_accuracy, test_accuracy))       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What Happened?\n",
    "\n",
    "By the 700th epoch, our model achieves 100% accuracy on the training data. However, it only classifies 75% of the test examples accurately. This is a clear case of overfitting. At a high level, let's reason about what went wrong. Because we have 7450 parameters and only 1000 data points, there are actually many settings of the parameters that could produce 100% accuracy on training data. \n",
    "\n",
    "To get some intuition imagine that we wanted to fit a dataset with 2 dimensional data and 2 data points. Our model has three degrees of freedom, and thus for any dataset can find an arbitrary number of separators that will perfectly caslsify our training points. Note below that we can produce completely orthogonal separators that both classify our training data perfectly. Even if it seems perposterous that they could both describe our training data well.\n",
    "\n",
    "![](img/overfitting-low-data.png)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def l2_penalty(params):\n",
    "    penalty = nd.zeros(shape=1)\n",
    "    for param in params:\n",
    "        penalty = penalty + nd.sum(param ** 2)\n",
    "    return penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regularization\n",
    "\n",
    "Now that we've characterized the problem of overfitting, we can begin talking about some solutions. \n",
    "Broadly speaking the family of techniques geared towards mitigating overfitting are referred to as *regularization*.\n",
    "The core idea is this: when a model is overfitting, its training error is substantially lower than its test error. We're already doing as well as we possibly can on the training data, but our test data performance leaves something to be desired. Typically, regularization techniques attempt to trade of of our training performance in exchange for lowering our test error. \n",
    "\n",
    "There are several straight forward techniques we might employ. Given the intuition from the previous chart, we might attempt to make our model less complex. One way to do this would be to lower the number of free parameters. For example, we could throw away some subset of the our input features (and thus the correpsonding parameters) that we thought were least informative. \n",
    "\n",
    "![](img/regularization-overfitting.png)\n",
    "\n",
    "Another approach is to limit the values that our weights might take. One common approach is to force the weights to take small values. \n",
    "[give more intuition with example of polynomial curve fitting]\n",
    "We can accomplish this by changing our optimization objective to penalize the value of our weights. \n",
    "The most popular regularizer is the $\\ell^2_2$ norm. For linear models, $\\ell^2_2$ regularization as the additional benefit that it makes the solution unique, even when our model is overparametrized.\n",
    "$$\\sum_{i}(\\hat{y}-y)^2 + \\lambda || \\textbf{w} ||^2_2$$\n",
    "Here, $||\\textbf{w}||$ is the $\\ell^2_2$ norm and $\\lambda$ is a hyper-parameter that determines how aggressively we want to push the weights towards 0.\n",
    "\n",
    "In code, we can express the $\\ell^2_2$ penalty succinctly:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Re-initializing the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for param in params:\n",
    "    param[:] = nd.random_normal(shape=param.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training L2-regularized logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed epoch 100. Loss: 421.122502903, Train_acc 0.898438, Test_acc 0.733398\n",
      "Completed epoch 200. Loss: 220.921968033, Train_acc 0.963867, Test_acc 0.786133\n",
      "Completed epoch 300. Loss: 121.568079227, Train_acc 0.981445, Test_acc 0.810547\n",
      "Completed epoch 400. Loss: 70.5325320933, Train_acc 0.989258, Test_acc 0.824219\n",
      "Completed epoch 500. Loss: 44.0794594871, Train_acc 0.991211, Test_acc 0.831055\n",
      "Completed epoch 600. Loss: 30.3448870191, Train_acc 0.992188, Test_acc 0.834961\n",
      "Completed epoch 700. Loss: 23.2058101207, Train_acc 0.991211, Test_acc 0.834961\n",
      "Completed epoch 800. Loss: 19.4902396174, Train_acc 0.991211, Test_acc 0.836914\n",
      "Completed epoch 900. Loss: 17.5540405403, Train_acc 0.991211, Test_acc 0.835938\n",
      "Completed epoch 1000. Loss: 16.5440581832, Train_acc 0.991211, Test_acc 0.838867\n"
     ]
    }
   ],
   "source": [
    "epochs = 1000\n",
    "moving_loss = 0.\n",
    "l2_strength = .1\n",
    "\n",
    "for e in range(epochs):\n",
    "    train_data.reset()\n",
    "    for i, batch in enumerate(train_data):\n",
    "        data = batch.data[0].as_in_context(ctx).reshape((-1,784))\n",
    "        label = batch.label[0].as_in_context(ctx)\n",
    "        label_one_hot = nd.one_hot(label, 10)\n",
    "        with autograd.record():\n",
    "            output = net(data)\n",
    "            loss = nd.sum(cross_entropy(output, label_one_hot)) + l2_strength * l2_penalty(params)\n",
    "        loss.backward()\n",
    "        SGD(params, .001)\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",
    "    if e % 100 == 99:\n",
    "        print(\"Completed epoch %s. Loss: %s, Train_acc %s, Test_acc %s\" % \n",
    "              (e+1, moving_loss, train_accuracy, test_accuracy))       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "By adding $L_2$ regularization we were able to increase the performance on test data from 75% accuracy to 83% accuracy. That's a 32% reduction in error. In a lot of applications, this big an improvement can make the difference between a viable product and useless system. Note that L2 regularization is just one of many ways of controlling capacity. Basically we assumed that small weight values are good. But there are many more ways to constrain the values of the weights:\n",
    "\n",
    "* We could require that the total sum of the weights is small. That is what $L_1$ regularization does via the penalty $\\sum_i |w_i|$. \n",
    "* We could require that the largest weight is not too large. This is what $L_\\infty$ regularization does via the penalty $\\max_i |w_i|$. \n",
    "* We could require that the number of nonzero weights is small, i.e. that the weight vectors are *sparse*. This is what the $L_0$ penalty does, i.e. $\\sum_i I\\{w_i \\neq 0\\}$. This penalty is quite difficult to deal with explicitly since it is nonsmooth. There is a lot of research that shows how to solve this problem approximatley using an $L_1$ penalty. \n",
    "\n",
    "![](img/regularization.png)\n",
    "\n",
    "From left to right: $L_2$ regularization, which constrains the parameters to a ball, $L_1$ regularization, which constrains the parameters to a diamond (for lack of a better name, this is often referred to as an $L_1$-ball), and $L_\\infty$ regularization, which constrains the parameters to a hypercube. \n",
    "\n",
    "All of this raises the question of **why** regularization is any good. After all, choice is good and giving our model more flexibility *ought* to be better (e.g. there are plenty of papers which show improvements on ImageNet using deeper networks). What is happening is somewhat more subtle. Allowing for many different parameter values allows our model to cherry pick a combination that is *just right* for all the training data it sees, without really learning the underlying mechanism. Since our observations are likely noisy, this means that we are trying to approximate the errors at least as much as we're learning what the relation between data and labels actually is. There is an entire field of statistics devoted to this issue - Computational Learning Theory. For now, a few simple rules of thumb suffice:\n",
    "\n",
    "* Fewer parameters tend to be better than more parameters.\n",
    "* Better engineering for a specific problem that takes the actual problem into account will lead to better models, due to the prior knowledge that data scientists have about the problem at hand.\n",
    "* $L_2$ is easier to optimize for than $L_1$. In particular, many optimizers will not work well out of the box for $L_1$. Using the latter requires something called *proximal operators*.\n",
    "* Dropout and other methods to make the model robust to perturbations in the data often work better than off-the-shelf $L_2$ regularization.\n",
    "\n",
    "We conclude with an [XKCD Cartoon](https://xkcd.com/882/) which captures the entire situation much better than the entire paragraph above. \n",
    "\n",
    "![](https://imgs.xkcd.com/comics/significant.png)"
   ]
  },
  {
   "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 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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}