02_neural_networks_python/01_introduction_mlp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MNIST handwritten digits classification with MLPs, by hand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Bethany Lusch adapting notebooks by Tanwi Mallick, Prasanna Balaprakash and Taylor Childers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goal of notebook:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we'll train a multi-layer perceptron model (a basic kind of neural network) to classify handwritten digits. We'll build up the code by hand. Next week, we show how this can be done using existing Python libraries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Machine learning task:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial works through a supervised learning problem, specifically classification.\n",
    "\n",
    "Imagine you are making a machine for the post office that will automatically sort mail by zip code. The MNIST dataset contains thousands of examples of handwritten numbers, with each digit labeled 0-9. We will use deep learning to create a function that classifies each image of one number as a digit 0-9.\n",
    "<img src=\"images/mnist_task.png\"  align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, the needed imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import tensorflow as tf\n",
    "\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MNIST data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll load the MNIST handwritten digits data set. The first time we may have to download the data, which can take a while.\n",
    "\n",
    "<img src=\"images/MnistExamples.png\"  align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MNIST is a popular dataset, so we can download it via the TensorFlow library. Note:\n",
    "- x is for the inputs (images of handwritten digits) and y is for the labels or outputs (digits 0-9)\n",
    "- We are given \"training\" and \"test\" datasets. Training datasets are used to fit the model. Test datasets are saved until the end, when we are satisfied with our model, to estimate how well our model generalizes to new data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(x_train, y_train), (x_test, y_test) = tf.keras.datasets.mnist.load_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do some pre-processing on the images: convert from integer to float32 and normalize the pixels to be within 0 to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = x_train.astype(numpy.float32)\n",
    "x_test  = x_test.astype(numpy.float32)\n",
    "\n",
    "x_train /= 255.\n",
    "x_test  /= 255.\n",
    "\n",
    "print(x_train.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training data (x_train) is a tensor of 60,000 images, each of size 28x28 pixels.\n",
    "\n",
    "For this notebook, we flatten each image to a vector, so x_train is a matrix of size (60000, 28*28)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = x_train.reshape(x_train.shape[0], numpy.prod(x_train[0,:,:].shape))\n",
    "x_test = x_test.reshape(x_test.shape[0], numpy.prod(x_test[0,:,:].shape))\n",
    "\n",
    "print(x_train.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y_train is a 60000-dimensional vector containing the correct classes (\"0\", \"1\", ..., \"9\") for each training sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = y_train.astype(numpy.int32)\n",
    "y_test  = y_test.astype(numpy.int32)\n",
    "\n",
    "print()\n",
    "print('MNIST data loaded: train:',len(x_train),'test:',len(x_test))\n",
    "print('X_train:', x_train.shape)\n",
    "print('y_train:', y_train.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a closer look. Here are the first 10 training digits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pltsize=1\n",
    "plt.figure(figsize=(10*pltsize, pltsize))\n",
    "\n",
    "for i in range(10):\n",
    "    plt.subplot(1,10,i+1)\n",
    "    plt.axis('off')\n",
    "    plt.imshow(numpy.reshape(x_train[i,:], (28, 28)), cmap=\"gray\")\n",
    "    plt.title('Class: '+str(y_train[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin with a simple linear model: linear regression, like last week. \n",
    "We add one complication: each example is a vector (flattened image), so the \"slope\" multiplication becomes a dot product.\n",
    "\n",
    "Note, like before, we consider multiple examples at once. input_images is a matrix where each row is an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_model(W, input_images):\n",
    "    # f(x) = xW returns m-length vector, where m is the number of examples\n",
    "    return numpy.dot(input_images, W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like last week, the loss is mean squared error (MSE):\n",
    "\n",
    "$\\large{MSE = \\frac{1}{n} \\sum_{i=1}^{n} (y - \\hat{y})^{2}}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate(W, input_images, true_labels):\n",
    "    predicted_label = linear_model(W, input_images)\n",
    "    MSE = numpy.mean((true_labels - predicted_label)**2) \n",
    "    return MSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Update weights using gradient decent \n",
    "$\\large{W = W - \\eta \\frac{\\partial J(W)}{\\partial W} }$,  where, $W$ is the network weight, $\\eta$ is the learning rate and $J(W)$ is the objective function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def learn(input_images, true_labels, current_W, learning_rate=0.000001):\n",
    "    # first we need dJW/dW where JW = MSE \n",
    "    n = input_images.shape[0] # get number of examples to average over\n",
    "    label_predictions = linear_model(current_W, input_images)\n",
    "    # calculate gradient: one entry per partial derivative for an entry in vector W\n",
    "    dJW_dW = (2./n) * numpy.dot(input_images.transpose(), label_predictions - true_labels)\n",
    "    # now we update W\n",
    "    new_W = current_W - (learning_rate * dJW_dW)  # gradient update step\n",
    "    return new_W"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, we usually don't use all of the training data to calculate each step. We use a random subset. This makes the steps faster and noisier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# At the moment we take the simple route and use a fixed subset. \n",
    "batch_size = 100\n",
    "\n",
    "x_train_batch = x_train[:batch_size, :]\n",
    "y_train_batch = y_train[:batch_size,numpy.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_features = x_train.shape[1] # this is the number of pixels\n",
    "\n",
    "# Randomly initialize W\n",
    "W = .01 * numpy.random.rand(num_features,1)\n",
    "\n",
    "# now iterate num_iters times, with the step size defined by learning_rate\n",
    "learning_rate = 0.0005  \n",
    "num_iters = 5000\n",
    "losses = numpy.zeros(num_iters,)\n",
    "\n",
    "for i in range(0, num_iters):\n",
    "    # all the magic here\n",
    "    W = learn(x_train_batch, y_train_batch, W, learning_rate)\n",
    "    losses[i] = evaluate(W, x_train_batch, y_train_batch)\n",
    "plt.plot(losses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check results so far"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, this simple linear model f(x) = xW is not very accurate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pltsize=1\n",
    "plt.figure(figsize=(10*pltsize, pltsize))\n",
    "predicted_labels = linear_model(W, x_train[:10,:])\n",
    "\n",
    "for i in range(10):\n",
    "    plt.subplot(1,10,i+1)\n",
    "    plt.axis('off')\n",
    "    plt.imshow(numpy.reshape(x_train[i,:], (28, 28)), cmap=\"gray\")\n",
    "    plt.title('%1.2f' % predicted_labels[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In-class exercise: implement \"accuracy\" - number of images correctly labeled"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Improving the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examples of basic ways to improve:\n",
    "- Add bias term: f(x) = xW + b \n",
    "- Reformulate as classification (output integers, not real numbers), like logistic regression\n",
    "- Minimize something other than mean squared error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sigmoid function encourages outputs of 0 and 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    z = 1/(1 + numpy.exp(-x))\n",
    "    return(z)\n",
    "\n",
    "def classification_model(A, b, input_images):\n",
    "    # f(x) = sigmoid(xA + b) returns m-length vector, where m is the number of examples\n",
    "    return sigmoid(numpy.dot(input_images, W) + b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = numpy.arange(-10, 10, step=.1)\n",
    "plt.plot(x, sigmoid(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To handle multiple classes, it's common to use a one-hot encoding:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one-hot encoding:\n",
    "nb_classes = 10\n",
    "y_train_onehot = tf.keras.utils.to_categorical(y_train, nb_classes)\n",
    "y_test_onehot = tf.keras.utils.to_categorical(y_test, nb_classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pltsize=1\n",
    "plt.figure(figsize=(10*pltsize, pltsize))\n",
    "\n",
    "for i in range(10):\n",
    "    plt.subplot(1,10,i+1)\n",
    "    plt.axis('off')\n",
    "    plt.imshow(numpy.reshape(x_train[i,:], (28, 28)), cmap=\"gray\")\n",
    "    plt.title('Class: '+str(y_train[i]))\n",
    "    print('Training sample',i,': class:',y_train[i], ', one-hot encoded:', y_train_onehot[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To handle multi-class classification, it's common to use softmax instead of sigmoid. It's related but forces the outputs to sum to 1, like a probability distribution. The class with the highest value is the prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def softmax(x):\n",
    "    normalization = numpy.sum(numpy.exp(x),axis=1)\n",
    "    z = numpy.exp(x) / normalization[:,None]\n",
    "    return(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Randomly initialize W and b\n",
    "W = .01 * numpy.random.rand(num_features,10)\n",
    "b = .01 * numpy.random.rand(10,)\n",
    "\n",
    "def classification_model(W, b, input_images):\n",
    "    # f(x) = softmax(xW + b) returns m x 10 matrix, where m is the number of examples and 10 is the number of classes\n",
    "    return softmax(numpy.dot(input_images, W) + b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To demonstrate softmax, we can apply this model with random W & b\n",
    "predicted_labels = classification_model(W, b, x_train[:5,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here the class \"probabilities\" might all be very similar since the model hasn't been trained\n",
    "print(predicted_labels[0,:])\n",
    "print(sum(predicted_labels[0,:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, for multi-class classification problem, it is common to minimize a different \"loss\" function instead of mean squared error, like categorical cross-entropy. You can read more [here](https://gombru.github.io/2018/05/23/cross_entropy_loss/). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above improvements are not enough to classify these images. We move to a nonlinear model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A neural network has multiple layers. A basic layer is $\\sigma(xA + b)$, where $\\sigma$ is a nonlinear \"activation function.\" An example neural network with two layers adds another affine transformation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$f(x) = \\sigma(\\sigma(xW_1 + b_1)W_2 + b_2)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is often drawn as a network like this:\n",
    "\n",
    "\n",
    "<img src=\"images/tiny_network.png\" width=\"300\" hight=\"300\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Writing it in numpy, with sigmoid as the activation function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonlinear_model(W1, W2, b1, b2, input_images):\n",
    "    return sigmoid(numpy.dot(sigmoid(numpy.dot(input_images, W1) + b1), W2) + b2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding another layer:\n",
    "\n",
    "$f(x) = \\sigma(\\sigma(\\sigma(xW_1 + b_1)W_2 + b_2)W_3 + b_3)$\n",
    "\n",
    "\n",
    "<img src=\"images/three_layer_network.png\" width=\"300\" hight=\"300\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonlinear_model(A1, A2, A3, b1, b2, b3, input_images):\n",
    "    return sigmoid(numpy.dot(sigmoid(numpy.dot(sigmoid(numpy.dot(input_images, A1) + b1), A2) + b2), A3) + b3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will learn about other activation functions (nonlinearities) next week."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimizing a nonlinear model: \n",
    "- We still want to use stochastic gradient descent (or a variant), but now the gradients are more complicated\n",
    "- The gradients can be calculated with calculus (chain rule!). To save on computation, we move backward through the layers, saving intermediate results for re-use. This is called **back-propagation**.\n",
    "- Applying the current network to the data is often called the \"forward pass,\" and calculating the gradients is called the \"backward pass.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonlinear neural networks can fit more complicated data than linear models. On the other hand, deep learning training can be tricky. \n",
    "\n",
    "1. Unlike linear regression, the objective function that you're minimizing (some measure of error) is non-convex, so there can be many local optima. As we learned about last week, the learning rate can help you jump into a new area, although too much jumping can be bad. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![multiple local minima](images/nonconvex.png)\n",
    "\n",
    "Image source: [firsttimeprogrammer.blogspot.com](http://firsttimeprogrammer.blogspot.com/2014/09/multivariable-gradient-descent.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some fancier versions of gradient descent optimization algorithms that are more effective, such as:\n",
    "- Adam \n",
    "- RMSprop\n",
    "- Adadelta\n",
    "- Adagrad\n",
    "\n",
    "For far more information, and some cool animations, see https://ruder.io/optimizing-gradient-descent/ or https://distill.pub/2017/momentum/. It will be easier to try them out next week when we are using TensorFlow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Neural networks can be overly flexible/complicated and \"overfit\" your data. This is like what happens if you fit a high-degree polynomial:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img src=\"images/bias_vs_variance.png\" width=\"800\" hight=\"500\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. To improve the generalization of our model on previously unseen data, we employ a technique known as regularization, which constrains our optimization problem in order to discourage complex models. Next week, we'll learn about Dropout for regularization. A simpler form is to add a penalty for large weights (\"L2 regularization\")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img src=\"images/test_data_rule.png\" width=\"800\" hight=\"500\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In-class exercise: split the training data into training & validation, and track validation loss during the training loop. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tip: https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html\n",
    "\n",
    "How do we know when to stop training? For example, you might stop when the validation loss stops improving."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the basics, we can experiment with more complicated networks. Rather than implement these all by hand, we will move to using existing Python packages next week. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are variants of \"universal approximation theorems\" roughly stating that there exists a nonlinear neural network with one hidden layer (possibly very wide) can fit an \"arbitrary\" nice/smooth function arbitrarily well. However, we can make the optimizaiton easier with fancier layers than \"fully connected,\" like convolutional layers, which we will learn about next week. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}