Part 5 Numerical Gradient Checking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align = 'center'> Neural Networks Demystified </h1>\n",
    "<h2 align = 'center'> Part 5: Numerical Gradient Checking </h2>\n",
    "\n",
    "\n",
    "<h4 align = 'center' > @stephencwelch </h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/pHMzNW8Agq4\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x105339e50>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('pHMzNW8Agq4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last time, we did a bunch of calculus to find the rate of change of our cost, J, with respect to our parameters, W. Although each calculus step was pretty straight forward, it’s still easy to make mistakes. What’s worse, is that our network doesn’t have a good way to tell us that it’s broken – code with incorrectly implemented gradients may appear to be functioning just fine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the most nefarious kind of error when building complex systems. Big, in-your-face errors suck initially, but it’s clear that you must fix this error for your work to succeed. More subtle errors can be more troublesome because they hide in your code and steal hours of your time, slowly degrading performance, while you wonder what the problem is. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A good solution here is to test the gradient computation part of our code, just as developer would unit test new portions of their code. We’ll combine a simple understanding of the derivative with some mild cleverness to perform numerical gradient checking. If our code passes this test, we can be quite confident that we have computed and coded up our gradients correctly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get started, let’s quickly review derivatives. Derivates tell us the slope, or how steep a function is. Once you’re familiar with calculus, it’s easy to take for granted the inner workings of the derivative - we just accept that the derivative of x^2 is 2x by the power rule. However, depending on how mean your calculus teacher was, you may have spent months not being taught the power rule, and instead required to compute derivatives using the definition. Taking derivatives this way is a bit tedious, but still important - it provides us a deeper understanding of what a derivative is, and it’s going to help us solve our current problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The definition of the derivative is really a glorified slope formula. The numerator gives us the change in y values, while the denominator is convenient way to express the change in x values. By including the limit, we are applying the slope formula across an infinitely small region – it’s like zooming in on our function, until it becomes linear. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The definition tells us to zoom in until our x distance is infinitely small, but computers can’t really handle infinitely small numbers, especially when they’re in the bottom parts of fractions - if we try to plug in something too small, we will quickly lose precision. The good news here is that if we plug in something reasonable small, we can still get surprisingly good numerical estimates of the derivative. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll modify our approach slightly by picking a point in the middle of the interval we would like to test, and call the distance we move in each direction epsilon. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s test our method with a simple function, x squared. We’ll choose a reasonable small value for epsilon, and compute the slope of x^2 at a given point by finding the function value just above and just below our test point. We can then compare our result to our symbolic derivative 2x, at the test point. If the numbers match, we’re in business!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "#Import Code from previous videos:\n",
    "from partFour import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "epsilon = 1e-4\n",
    "x = 1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "numericalGradient = (f(x+epsilon)- f(x-epsilon))/(2*epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.9999999999996696, 3.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numericalGradient, 2*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add helper functions to our neural network class: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Neural_Network(object):\n",
    "    def __init__(self):        \n",
    "        #Define Hyperparameters\n",
    "        self.inputLayerSize = 2\n",
    "        self.outputLayerSize = 1\n",
    "        self.hiddenLayerSize = 3\n",
    "        \n",
    "        #Weights (parameters)\n",
    "        self.W1 = np.random.randn(self.inputLayerSize,self.hiddenLayerSize)\n",
    "        self.W2 = np.random.randn(self.hiddenLayerSize,self.outputLayerSize)\n",
    "        \n",
    "    def forward(self, X):\n",
    "        #Propogate inputs though network\n",
    "        self.z2 = np.dot(X, self.W1)\n",
    "        self.a2 = self.sigmoid(self.z2)\n",
    "        self.z3 = np.dot(self.a2, self.W2)\n",
    "        yHat = self.sigmoid(self.z3) \n",
    "        return yHat\n",
    "        \n",
    "    def sigmoid(self, z):\n",
    "        #Apply sigmoid activation function to scalar, vector, or matrix\n",
    "        return 1/(1+np.exp(-z))\n",
    "    \n",
    "    def sigmoidPrime(self,z):\n",
    "        #Gradient of sigmoid\n",
    "        return np.exp(-z)/((1+np.exp(-z))**2)\n",
    "    \n",
    "    def costFunction(self, X, y):\n",
    "        #Compute cost for given X,y, use weights already stored in class.\n",
    "        self.yHat = self.forward(X)\n",
    "        J = 0.5*sum((y-self.yHat)**2)\n",
    "        return J\n",
    "        \n",
    "    def costFunctionPrime(self, X, y):\n",
    "        #Compute derivative with respect to W and W2 for a given X and y:\n",
    "        self.yHat = self.forward(X)\n",
    "        \n",
    "        delta3 = np.multiply(-(y-self.yHat), self.sigmoidPrime(self.z3))\n",
    "        dJdW2 = np.dot(self.a2.T, delta3)\n",
    "        \n",
    "        delta2 = np.dot(delta3, self.W2.T)*self.sigmoidPrime(self.z2)\n",
    "        dJdW1 = np.dot(X.T, delta2)  \n",
    "        \n",
    "        return dJdW1, dJdW2\n",
    "    \n",
    "    #Helper Functions for interacting with other classes:\n",
    "    def getParams(self):\n",
    "        #Get W1 and W2 unrolled into vector:\n",
    "        params = np.concatenate((self.W1.ravel(), self.W2.ravel()))\n",
    "        return params\n",
    "    \n",
    "    def setParams(self, params):\n",
    "        #Set W1 and W2 using single paramater vector.\n",
    "        W1_start = 0\n",
    "        W1_end = self.hiddenLayerSize * self.inputLayerSize\n",
    "        self.W1 = np.reshape(params[W1_start:W1_end], (self.inputLayerSize , self.hiddenLayerSize))\n",
    "        W2_end = W1_end + self.hiddenLayerSize*self.outputLayerSize\n",
    "        self.W2 = np.reshape(params[W1_end:W2_end], (self.hiddenLayerSize, self.outputLayerSize))\n",
    "        \n",
    "    def computeGradients(self, X, y):\n",
    "        dJdW1, dJdW2 = self.costFunctionPrime(X, y)\n",
    "        return np.concatenate((dJdW1.ravel(), dJdW2.ravel()))\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the same approach to numerically evaluate the gradient of our neural network. It’s a little more complicated this time, since we have 9 gradient values, and we’re interested in the gradient of our cost function. We’ll make things simpler by testing one gradient at a time. We’ll “perturb” each weight - adding epsilon to the current value and  computing the cost function, subtracting epsilon from the current value and computing the cost function, and then computing the slope between these two values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def computeNumericalGradient(N, X, y):\n",
    "        paramsInitial = N.getParams()\n",
    "        numgrad = np.zeros(paramsInitial.shape)\n",
    "        perturb = np.zeros(paramsInitial.shape)\n",
    "        e = 1e-4\n",
    "\n",
    "        for p in range(len(paramsInitial)):\n",
    "            #Set perturbation vector\n",
    "            perturb[p] = e\n",
    "            N.setParams(paramsInitial + perturb)\n",
    "            loss2 = N.costFunction(X, y)\n",
    "            \n",
    "            N.setParams(paramsInitial - perturb)\n",
    "            loss1 = N.costFunction(X, y)\n",
    "\n",
    "            #Compute Numerical Gradient\n",
    "            numgrad[p] = (loss2 - loss1) / (2*e)\n",
    "\n",
    "            #Return the value we changed to zero:\n",
    "            perturb[p] = 0\n",
    "            \n",
    "        #Return Params to original value:\n",
    "        N.setParams(paramsInitial)\n",
    "\n",
    "        return numgrad "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll repeat this process across all our weights, and when we’re done we’ll have a numerical gradient vector, with the same number of values as we have weights. It’s this vector we would like to compare to our official gradient calculation. We see that our vectors appear very similar, which is a good sign, but we need to quantify just how similar they are. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "NN = Neural_Network()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.04889799, -0.06551702,  0.02013916,  0.04615717, -0.05295928,\n",
       "        0.01661075, -0.25285186, -0.13721988, -0.10629924])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numgrad = computeNumericalGradient(NN, X, y)\n",
    "numgrad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.04889799, -0.06551702,  0.02013916,  0.04615717, -0.05295928,\n",
       "        0.01661075, -0.25285186, -0.13721988, -0.10629924])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad = NN.computeGradients(X,y)\n",
    "grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A nice way to do this is to divide the norm of the difference by the norm of the sum of the vectors we would like to compare. Typical results should be on the order of 10^-8 or less if you’ve computed your gradient correctly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0134439421141014e-10"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(grad-numgrad)/norm(grad+numgrad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that’s it, we can now check our computations and eliminate gradient errors before they become a problem. Next time we’ll train our Neural Network. "
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}