4_Neural_Network_Tutorial_Backpropagation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 707,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Pre-requisites\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "# To clear print buffer\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing code from previous tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initializing weight matrices from layer sizes\n",
    "def initializeWeights(layers):\n",
    "    weights = [np.random.randn(o, i+1) for i, o in zip(layers[:-1], layers[1:])]\n",
    "    return weights\n",
    "\n",
    "# Add a bias term to every data point in the input\n",
    "def addBiasTerms(X):\n",
    "        # Make the input an np.array()\n",
    "        X = np.array(X)\n",
    "        \n",
    "        # Forcing 1D vectors to be 2D matrices of 1xlength dimensions\n",
    "        if X.ndim==1:\n",
    "            X = np.reshape(X, (1, len(X)))\n",
    "        \n",
    "        # Inserting bias terms\n",
    "        X = np.insert(X, 0, 1, axis=1)\n",
    "        \n",
    "        return X\n",
    "\n",
    "# Sigmoid function\n",
    "def sigmoid(a):\n",
    "    return 1/(1 + np.exp(-a))\n",
    "\n",
    "# Forward Propagation of outputs\n",
    "def forwardProp(X, weights):\n",
    "    # Initializing an empty list of outputs\n",
    "    outputs = []\n",
    "    \n",
    "    # Assigning a name to reuse as inputs\n",
    "    inputs = X\n",
    "    \n",
    "    # For each layer\n",
    "    for w in weights:\n",
    "        # Add bias term to input\n",
    "        inputs = addBiasTerms(inputs)\n",
    "        \n",
    "        # Y = Sigmoid ( X .* W^T )\n",
    "        outputs.append(sigmoid(np.dot(inputs, w.T)))\n",
    "        \n",
    "        # Input of next layer is output of this layer\n",
    "        inputs = outputs[-1]\n",
    "        \n",
    "    return outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training Neural Networks\n",
    "\n",
    "$$ Y^{(l)}_{n{\\times}o_{l}} = Sigmoid\\;(\\;X^{(l)}_{n{\\times}i_{l}} \\; .* \\; W^{(l)}{^{T}}_{i_{l}{\\times}o_{l}}) \\;\\;\\;\\;\\;\\;-------------(1)$$\n",
    "\n",
    "Neural networks are advantageous when we are able to compute that $W$ which satisfies $Y = Sigmoid(X\\cdot*W)$, for given $X$ and $Y$ (in supervised training).\n",
    "\n",
    "But, since there are so many weights (for bigger networks), it is time-intensive to algebraically solve the above equation. (Something like $W = X^{-1} \\;.*\\; Sigmoid^{-1}(Y)$...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set W to minimize cost (computationally intensive)\n",
    "\n",
    "A quicker way to compute W would be to randomly initialize it, and keep updating its value in such a way as to decrease the cost of the neural network.\n",
    "\n",
    "Define the cost as the mean squared error of the output of the neural network:\n",
    "\n",
    "$$error = yPred-Y$$\n",
    "\n",
    "Here, $yPred$ = ``forwardProp``$(X)$, and $Y$ is the desired output value from the neural network.\n",
    "\n",
    "$$Cost \\; J = \\frac{1}{2} \\sum \\limits_{n} \\frac{ {\\left( error \\right)}^2 }{n} = \\frac{1}{2} \\sum \\limits_{n} \\frac{ {\\left( yPred-Y \\right)}^2 }{n}$$\n",
    "\n",
    "Once we have initialized W, we need to change it such that J is minimized.\n",
    "\n",
    "The best way to minimize J w.r.t. W, is to partially derive J w.r.t. W and equate it to 0: $\\frac{{\\partial}J}{{\\partial}W} = 0$. But, this is computationally intensive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 433,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute COST (J) of Neural Network\n",
    "def nnCost(weights, X, Y):\n",
    "    # Calculate yPred\n",
    "    yPred = forwardProp(X, weights)[-1]\n",
    "    \n",
    "    # Compute J\n",
    "    J = 0.5*np.sum((yPred-Y)**2)/len(Y)\n",
    "    \n",
    "    return J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 434,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initialize network\n",
    "layers = [2, 2, 1]\n",
    "weights = initializeWeights(layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 435,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Declare input and desired output for AND gate\n",
    "X = np.array([[0,0], [0,1], [1,0], [1,1]])\n",
    "Y = np.array([[0], [0], [0], [1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 436,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.284231765606\n"
     ]
    }
   ],
   "source": [
    "# Cost\n",
    "J = nnCost(weights, X, Y)\n",
    "print(J)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomly initialize W, change it to decrease cost (more feasible)\n",
    "\n",
    "Instead, we initialize $W$ by randomly sampling from a standard normal distribution, and then keep changing $W$ so as to decrease the cost $J$.\n",
    "\n",
    "But what value to change $W$ by? To find out, let us focus on the weights of one of the neurons in the last layer, $W^{(L)}_{[k]}$, differentiate $J$ by it to see what we get:\n",
    "\n",
    "$$\\frac{ {\\partial}J} {{\\partial}W^{(L)}_{[k]} }=\\frac{\\partial}{{\\partial}W^{(L)}_{[k]}}\\left(\\frac{1}{2}\\sum\\limits_{n}{\\frac{ {\\left( yPred-Y \\right)}^2 }{n} }\\right)=\\frac{1}{2*n}\\sum\\limits_{n} \\left( \\frac{\\partial} {{\\partial}W^{(L)}_{[k]}} (yPred-Y)^2 \\right)=\\frac{1}{n}\\sum\\limits_{n} \\left( (yPred-Y) * \\frac {{\\partial} \\; yPred} { {\\partial}W^{(L)}_{[k]} } \\right)$$\n",
    "\n",
    "$$\\Rightarrow \\frac{ {\\partial}J} {{\\partial}W^{(L)}_{[k]} } = \\frac{1}{n}\\sum\\limits_{n} \\left( (error) * \\frac {{\\partial} \\; yPred} { {\\partial}W^{(L)}_{[k]} }  \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above equation tells us how $J$ changes by changing $W^{(L)}_{[k]}$. Approximating it for numerical analysis:\n",
    "\n",
    "$${\\Delta}J ={{\\Delta}W^{(L)}_{[k]}} * \\left[ \\frac{1}{n}\\sum\\limits_{n} \\left( (error) * \\frac {{\\partial} \\; yPred} { {\\partial}W^{(L)}_{[k]} } \\right) \\right] \\;\\;\\;\\;\\;\\;-------------(2)$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change $W^{(L)}_{[k]}$ so that $J$ always decreases\n",
    "\n",
    "If we ensure that ${\\Delta}W^{(L)}_{[k]}$ is equal to $-\\left[ \\frac{1}{n}\\sum\\limits_{n} \\left( (error) * \\frac {{\\partial} \\; yPred} { {\\partial}W^{(L)}_{[k]} } \\right) \\right]$, we see that ${\\Delta}J = {\\Delta}W^{(L)}_{[k]}*\\left(-\\left[{\\Delta}W^{(L)}_{[k]}\\right]\\right) = -\\left[{\\Delta}W^{(L)}_{[k]}\\right]^{2} \\Rightarrow$ negative! \n",
    "\n",
    "Thus, we decide to change $W^{(L)}_{[k]}$ by that amount which ensures $J$ always decreases!\n",
    "\n",
    "$${\\Delta}W^{(L)}_{[k]} = -\\left[ \\frac{1}{n}\\sum\\limits_{n} \\left( (error) * \\frac {{\\partial} \\; yPred} { {\\partial}W^{(L)}_{[k]} } \\right) \\right] \\;\\;\\;\\;\\;\\;-------------(3)$$ \n",
    "\n",
    "So, for each weight in the last layer, that ${\\Delta}W^{(L)}_{[k]}$ which shall (for sure) decrease J can be computed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent\n",
    "\n",
    "If we update each weight as $W^{(L)}_{[k]} \\leftarrow W^{(L)}_{[k]} + {\\Delta}W^{(L)}_{[k]}$, it is guaranteed that with the new weights, the neural network shall produce outputs that are closer to the desired output.\n",
    "\n",
    "This is how to train a neural network - randomly initialize $W$, iteratively change $W$ according to eq (3).\n",
    "\n",
    "**This is called Gradient Descent.**\n",
    "\n",
    "One way to think about this is - assuming the graph of $J$ vs. $W$ is like an upturned hill, we are slowly descending down the hill by changing $W$, to the point where $J$ is minimum.\n",
    "\n",
    "J is (sort of) a quadratic function on W, so we can assume it's (sort of) like an upturned hill."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing ${\\Delta}W^{(L)}$ of last layer\n",
    "\n",
    "To compute ${\\Delta}W$, we need to compute $error$ and $\\frac{{\\partial}\\;yPred}{{\\partial}W^{(L)}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Computing error\n",
    "\n",
    "$ error = yPred - Y = $ ``forwardProp``$(X) - Y \\;\\;\\;\\;\\;\\;-------------(4)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, suppose we want to compute those $W$'s in a 3-neuron network that are able to perform AND logic on two inputs.\n",
    "\n",
    "Here, for $X = \\left[\\begin{array}{c}(0,0)\\\\(0,1)\\\\(1,0)\\\\(1,1)\\end{array}\\right]$, $Y = \\left[\\begin{array}{c}0\\\\0\\\\0\\\\1\\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 686,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "weights:\n",
      "1\n",
      "(2, 3)\n",
      "[[-0.87271574  0.35621485  0.95252276]\n",
      " [-0.61981924 -1.49164222  0.55011796]]\n",
      "2\n",
      "(1, 3)\n",
      "[[-1.57656753 -1.10359895 -0.34594249]]\n"
     ]
    }
   ],
   "source": [
    "# Initialize network\n",
    "layers = [2, 2, 1]\n",
    "weights = initializeWeights(layers)\n",
    "\n",
    "print(\"weights:\")\n",
    "for i in range(len(weights)):\n",
    "    print(i+1); print(weights[i].shape); print(weights[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our weights have been randomly initialized. Let us see what yPred they give:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 687,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Declare input and desired output for AND gate\n",
    "X = np.array([[0,0], [0,1], [1,0], [1,1]])\n",
    "Y = np.array([[0], [0], [0], [1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 688,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "outputs\n",
      "[array([[ 0.29468953,  0.34982256],\n",
      "       [ 0.51994117,  0.48258173],\n",
      "       [ 0.37367081,  0.10798781],\n",
      "       [ 0.60731071,  0.17345395]]), array([[ 0.11682925],\n",
      "       [ 0.08969868],\n",
      "       [ 0.11646832],\n",
      "       [ 0.09056134]])]\n"
     ]
    }
   ],
   "source": [
    "# Calculate outputs at each layer by forward propagation\n",
    "outputs = forwardProp(X, weights)\n",
    "print(\"outputs\"); print(outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 689,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 1)\n",
      "[[ 0.11682925]\n",
      " [ 0.08969868]\n",
      " [ 0.11646832]\n",
      " [ 0.09056134]]\n"
     ]
    }
   ],
   "source": [
    "# Calculate yPred as the last output from forward propagation\n",
    "yPred = outputs[-1]\n",
    "print(yPred.shape); print(yPred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 690,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 1)\n",
      "[[ 0.11682925]\n",
      " [ 0.08969868]\n",
      " [ 0.11646832]\n",
      " [-0.90943866]]\n"
     ]
    }
   ],
   "source": [
    "# Error = yPred - Y\n",
    "error = yPred - Y\n",
    "print(error.shape); print(error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Computing $\\frac{{\\partial}\\;yPred}{{\\partial}W^{(L)}_{[k]}}$\n",
    "\n",
    "From eq. (1), $yPred$ can be written as:\n",
    "\n",
    "$$yPred = Sigmoid(X^{(L)}\\;.*\\;W^{(L)}{^{T}})$$\n",
    "\n",
    "So,\n",
    "\n",
    "$$\\frac{{\\partial}\\;yPred}{{\\partial}W^{(L)}_{[k]}} = \\frac{{\\partial}}{{\\partial}W^{(L)}_{[k]}}\\left(Sigmoid\\left(X^{(L)}.*W^{(L)}{^{T}}\\right)\\right) = Sigmoid^{'}\\left(X^{(L)}.*W^{(L)}{^{T}}\\right)*\\left(\\frac{{\\partial}}{{\\partial}W^{(L)}_{[k]}}\\left((X^{(L)}.*W^{(L)}{^{T}})\\right)\\right)$$\n",
    "\n",
    "Here, $yPred$ is an $o_{L}$-dimensional number, and $W^{(L)}_{[k]}$ only affect the $k$-th dimension of $yPred$, i.e. $yPred_{[k]}$. So,\n",
    "\n",
    "$$\\frac{{\\partial}\\;yPred_{[k]}}{{\\partial}W^{(L)}_{[k]}} = Sigmoid^{'}\\left(X^{(L)}.*W^{(L)}_{[k]}\\right)*\\left(\\frac{{\\partial}}{{\\partial}W^{(L)}_{[k]}}\\left((X^{(L)}.*W^{(L)}_{[k]})\\right)\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### - Computing $Sigmoid^{'}\\left(X^{(L)}.*W^{(L)}_{[k]}\\right)$\n",
    "\n",
    "It can be verified that $Sigmoid^{'}(a) = Sigmoid(a)*(1-Sigmoid(a))$. Thus, $Sigmoid^{'}(X^{(L)}.*W^{(L)}_{[k]}{^{T}}) = yPred_{[k]}*(1 - yPred_{[k]})$. So,\n",
    "\n",
    "$$\\frac{{\\partial}\\;yPred_{[k]}}{{\\partial}W^{(L)}_{[k]}} = \\left(yPred_{[k]}*(1 - yPred_{[k]})*\\left(\\frac{{\\partial}(X^{(L)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L)}_{[k]}}\\right)\\right)$$\n",
    "\n",
    "$${\\Delta}W^{(L)}_{[k]} = -\\left[\\frac{1}{n}\\sum\\limits_{n}\\left(error_{[k]}*yPred_{[k]}*(1 - yPred_{[k]})*\\left(\\frac{{\\partial}(X^{(L)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L)}_{[k]}}\\right)\\right)\\right]  \\;\\;\\;\\;\\;\\;-------------(5)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### - Computing $\\frac{{\\partial}}{{\\partial}W^{(L)}_{[k]}}((X^{(L)}.*W^{(L)}_{[k]}))$\n",
    "\n",
    "It can be seen that $\\frac{{\\partial}}{{\\partial}W^{(L)}_{[k]}}((X^{(L)}.*W^{(L)}_{[k]})) = X^{(L)}$\n",
    "\n",
    "We also know that $X^{(L)} = \\left[ \\begin{array}{c} 1 & Y^{(L-1)} \\end{array} \\right]_{n{\\times}i_{L}}$, and $Y^{(L-1)}$ have been computed during Forward Propagation. So,\n",
    "\n",
    "$$\\frac{{\\partial}\\;yPred_{[k]}}{{\\partial}W^{(L)}} = (yPred_{[k]}*(1-yPred_{[k]}))*X^{(L)} $$\n",
    "\n",
    "$${\\Delta}W^{(L)}_{[k]} = -\\left[\\frac{1}{n}\\sum\\limits_{n}\\left(error_{[k]}*yPred_{[k]}*(1 - yPred_{[k]})*X^{(L)}\\right)\\right]\\;\\;\\;\\;\\;\\;-------------(6)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining terms to simplify computation\n",
    "\n",
    "Here, dimension of $error$, $yPred$, and $(1-yPred)$ is $n{\\times}o_{L}$, while that of $X^{(L)}$ is $n{\\times}i_{L}$. A little thought has to be given towards how those quantities are multiplied.\n",
    "\n",
    "First of all, we can combine the mentioned three into one and call it $\\delta$.\n",
    "\n",
    "$${\\delta}_{n{\\times}o_{L}} = error_{n{\\times}o_{L}}*yPred_{n{\\times}o_{L}}*(1-yPred)_{n{\\times}o_{L}} \\;\\;\\;\\;\\;\\;-----(7)$$\n",
    "\n",
    "$${\\Delta}W^{(L)}_{[k]} = -\\left[\\frac{1}{n}\\sum\\limits_{n}\\left({\\delta}_{[k]}*X^{(L)}\\right)\\right]$$\n",
    "\n",
    "We can now combine calculations of all dimensions into this matrix operation: (We will figure out the matrix dimensions below)\n",
    "\n",
    "$${\\Delta}W^{(L)} = -\\left[\\frac{1}{n}\\sum\\limits_{n}\\left({\\delta}*X^{(L)}\\right)\\right]$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way of figuring out how $\\delta$ and $X^{(L)}$ are combined is to see that the dimension of ${\\Delta}W$ is $o_{L}{\\times}i_{L}$, dimension of $\\delta$ is $n{\\times}o_{L}$, and the dimension of $X^{(L)}$ is $n{\\times}i_{L}$.\n",
    "\n",
    "Clearly, the $\\sum\\limits_{n}\\left({\\delta}*X^{(L)}\\right)$ term, when considered for all the weights, is equal to $\\delta^{T}_{o_{L}{\\times}n}\\;.*\\;X^{(L)}_{n{\\times}i_{L}}$, the summation over $n$ being taken care of by the dot product, and the output dimension ${o_{L}{\\times}i_{L}}$ matches that of $W^{(L)}$.\n",
    "\n",
    "Hence, using matrix operations, ${\\Delta}W^{(L)}$ can be found as:\n",
    "\n",
    "$${\\Delta}W^{(L)}_{{o_{L}{\\times}i_{L}}} = -\\frac{1}{n}\\left({\\delta}^{T}{_{o_{L}{\\times}n}}\\;.*\\;X^{(L)}_{n{\\times}i_{L}}\\right) \\;\\;\\;\\;\\;\\;-------------(8)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 691,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 1)\n",
      "[[ 0.01205446]\n",
      " [ 0.00732415]\n",
      " [ 0.01198499]\n",
      " [-0.07490136]]\n"
     ]
    }
   ],
   "source": [
    "# Calculate delta for the last layer\n",
    "delta = np.multiply(np.multiply(error, yPred), 1-yPred)\n",
    "print(delta.shape); print(delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 692,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 3)\n",
      "[[ 1.          0.29468953  0.34982256]\n",
      " [ 1.          0.51994117  0.48258173]\n",
      " [ 1.          0.37367081  0.10798781]\n",
      " [ 1.          0.60731071  0.17345395]]\n"
     ]
    }
   ],
   "source": [
    "# Find input to the last layer\n",
    "xL = addBiasTerms(outputs[-2])\n",
    "print(xL.shape); print(xL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 693,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 3)\n",
      "[[ 0.01088444  0.00841238  0.00098657]]\n"
     ]
    }
   ],
   "source": [
    "# Find deltaW for last layer\n",
    "deltaW = -np.dot(delta.T, xL)/len(Y)\n",
    "print(deltaW.shape); print(deltaW)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 694,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "old weights:\n",
      "1\n",
      "(2, 3)\n",
      "[[-0.87271574  0.35621485  0.95252276]\n",
      " [-0.61981924 -1.49164222  0.55011796]]\n",
      "2\n",
      "(1, 3)\n",
      "[[-1.57656753 -1.10359895 -0.34594249]]\n",
      "new weights:\n",
      "1\n",
      "(2, 3)\n",
      "[[-0.87271574  0.35621485  0.95252276]\n",
      " [-0.61981924 -1.49164222  0.55011796]]\n",
      "2\n",
      "(1, 3)\n",
      "[[-1.5656831  -1.09518657 -0.34495592]]\n",
      "old cost:\n",
      "0.107792308277\n",
      "new cost:\n",
      "0.107601673739\n"
     ]
    }
   ],
   "source": [
    "# Checking cost of neural network before and after change in W^{L}\n",
    "newWeights = [np.array(w) for w in weights]\n",
    "newWeights[-1] += deltaW\n",
    "\n",
    "print(\"old weights:\")\n",
    "for i in range(len(weights)):\n",
    "    print(i+1); print(weights[i].shape); print(weights[i])\n",
    "\n",
    "print(\"new weights:\")\n",
    "for i in range(len(newWeights)):\n",
    "    print(i+1); print(newWeights[i].shape); print(newWeights[i])\n",
    "\n",
    "print(\"old cost:\"); print(nnCost(weights, X, Y))\n",
    "print(\"new cost:\"); print(nnCost(newWeights, X, Y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Congratulations! You've just learned how to back propagate!**\n",
    "(1 layer only)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Back-propagation through layers\n",
    "\n",
    "For the last layer, according to eq. (5),\n",
    "$${\\Delta}W^{(L)}_{[k]} = -\\frac{1}{n}\\sum\\limits_{n}\\left(error_{[k]}*yPred_{[k]}*(1 - yPred_{[k]})*\\left(\\frac{{\\partial}(X^{(L)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L)}_{[k]}}\\right)\\right) = -\\frac{1}{n}\\sum\\limits_{n}\\left(\\delta^{(L)}_{[k]}*\\frac{{\\partial}(X^{(L)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L)}_{[k]}}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing for Layer L-1\n",
    "\n",
    "If we go back one more layer to find out ${\\Delta}W$ for the $p^{th}$ neuron in the $(L-1)^{th}$ layer, backpropagated from the $k^{th}$ neuron in the $L^{th}$ layer, noting that $X^{L} = Y^{L-1}$:\n",
    "\n",
    "$${\\Delta}W^{(L-1)}_{[p]from[k]} = -\\frac{1}{n}\\sum\\limits_{n}\\left(\\delta^{(L)}_{[k]}*\\frac{{\\partial}(X^{(L)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L-1)}_{[p]}}\\right) = -\\frac{1}{n}\\sum\\limits_{n}\\left(\\delta^{(L)}_{[k]}*\\frac{{\\partial}(Y^{(L-1)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L-1)}_{[p]}}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, $Y^{(L-1)}$ is the collected output of the penultimate layer, i.e. the collected output of all neurons in the penultimate layer. $W^{(L-1)}_{[p]}$ is the weight matrix of the $p^{th}$ neuron in the penultimate layer. So,\n",
    "\n",
    "$$\\frac{{\\partial}(Y^{(L-1)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L-1)}_{[p]}} = \\frac{{\\partial}((Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[0]}) + Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[1]}) + ... + Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[p]}) + ... + Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[i_{L-1}-1]})).*W^{(L)}_{[k]})}{{\\partial}W^{(L-1)}_{[p]}}$$\n",
    "\n",
    "We know that change in $W^{(L-1)}_{[p]}$ does not affect $W^{(L)}$ or any $W^{(L-1)}$ weight matrix other than $W^{(L-1)}_{[p]}$. So:\n",
    "\n",
    "$$\\frac{{\\partial}(Y^{(L-1)}.*W^{(L)}_{[k]})}{{\\partial}W^{(L-1)}_{[p]}} = \\frac{{\\partial}(Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[p]}))}{{\\partial}W^{(L-1)}_{[p]}}$$\n",
    "\n",
    "$${\\Delta}W^{(L-1)}_{[p]from[k]} = -\\frac{1}{n}\\sum\\limits_{n}\\left(\\delta^{(L)}_{[k]}*W^{(L)}_{[k]}*\\frac{{\\partial}\\;(Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[p]}))}{{\\partial}W^{(L-1)}_{[p]}}\\right)$$\n",
    "\n",
    "(Ignoring dimensions for now)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know how this goes now.\n",
    "\n",
    "$$\\frac{{\\partial}\\;(Sigmoid(X^{(L-1)}.*W^{(L-1)}_{[p]}))}{{\\partial}W^{(L-1)}_{[p]}} = Sigmoid^{'}(X^{(L-1)}.*W^{(L-1)}_{[p]})*\\frac{{\\partial}(X^{(L-1)}.*W^{(L-1)}_{[p]})}{{\\partial}W^{(L-1)}_{[p]}} = Y^{(L-1)}_{[p]}*(1 - Y^{(L-1)}_{[p]}))*X^{(L-1)}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus,\n",
    "\n",
    "$${\\Delta}W^{(L-1)}_{[p]from[k]} = -\\left[\\frac{1}{n}\\sum\\limits_{n}(\\delta^{(L)}_{[k]}*W^{(L)}_{[k]}*(Y^{(L-1)}_{[p]}*(1 - Y^{(L-1)}_{[p]})*X^{(L-1)})\\right]$$\n",
    "\n",
    "We need to take care of the dimensions. Here, there are two parts: $\\delta^{(L)}_{[k]}*W^{(L)}_{[k]}$, which is only concerned with the $L^{th}$ layer, and $Y^{(L-1)}_{[p]}*(1 - Y^{(L-1)}_{[p]}))*X^{(L-1)}$, which is only concerned with the $(L-1)^{th}$ layer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) Back-propagation Error\n",
    "\n",
    "We can observe here that the terms $\\delta^{(L)}_{[k]}$ and $W^{(L)}_{[k]}$ are back-propagated from the $k^{th}$ neuron of the final layer. Let's combine them and call it the back-propagated error:\n",
    "$$bpError^{(L-1)}_{[k]}{_{n{\\times}i_{L}}} = \\delta^{(L)}_{[k]}*W^{(L)}_{[k]}$$\n",
    "\n",
    "We know that $\\delta^{(L)}{_{n{\\times}o_{L}}}*W^{(L)}{_{o_{L}{\\times}{i_{L}}}}$ is a matrix of dimensions $n{\\times}i_{L} = n{\\times}o_{(L-1)}$, which is the sum of the backprop errors from each neuron in the final layer. Thus,\n",
    "\n",
    "$$bpError^{(L-1)}{_{n{\\times}i_{L}}} = \\delta^{(L)}*W^{(L)} \\;\\;\\;\\;\\;\\;--------------(9)$$\n",
    "\n",
    "We see that for a neuron in the $(L-1)^{th}$ layer, the total error back-propagated to it is the sum of the back-propagated errors from each of the neurons connected to it in the $L^{th}$ layer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, instead of ${\\Delta}W^{(L-1)}_{[p]from[k]}$ from the $k^{th}$ neuron, we can directly consider ${\\Delta}W^{(L-1)}_{[p]}$:\n",
    "\n",
    "$${\\Delta}W^{(L-1)}_{[p]} = -\\left[\\frac{1}{n}\\sum\\limits_{n}(bpError^{(L-1)}_{[p]}*(Y^{(L-1)}_{[p]}*(1 - Y^{(L-1)}_{[p]})*X^{(L-1)})\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) The $Y*(1-Y)*X$ term\n",
    "\n",
    "We can convert the $Y*(1-Y)*X$ term into a matrix operation, with summation over $n$ inherently taken care. Directly considering $Y$ instead of $Y_{[p]}$:\n",
    "\n",
    "$$Y^{(L-1)}*(1 - Y^{(L-1)}))*X^{(L-1)} == (Y^{(L-1)}.*(1 - Y^{(L-1)}))^{T}{_{o_{(L-1)}{\\times}n}} * X^{(L-1)}{_{n{\\times}i_{(L-1)}}}$$\n",
    "\n",
    "We can see that the resultant matrix has the same dimensions as $W^{(L-1)}$ : $o_{(L-1)}{\\times}i_{(L-1)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combining the two\n",
    "\n",
    "To combine $bpError$ and the $Y*(1-Y)*X$ terms, for consistency in dimensions, we need to first dot multiply $bpError_{n{\\times}o_{(L-1)}}$ with $Y^{(L-1)}_{n{\\times}o_{(L-1)}}.*(1 - Y^{(L-1)})_{n{\\times}o_{(L-1)}}$, and then multiply the transpose of that with X.\n",
    "\n",
    "Thus,\n",
    "\n",
    "$${\\Delta}W^{(L-1)}_{o_{(L-1)}{\\times}i_{(L-1)}} = -\\left[\\frac{1}{n}((bpError^{(L-1)}.*Y^{(L-1)}.*(1 - Y^{(L-1)}))^{T} _{o_{(L-1)}{\\times}n}* X^{(L-1)}_{n{\\times}i_{(L-1)}}\\right] \\;\\;\\;\\;\\;\\;--------------(10)$$\n",
    "\n",
    "(Summation across $n$ is taken care of within the matrix multiplication)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simplifying to matrix operation of any layer $l$\n",
    "\n",
    "Just as we had done for the final layer, from equation 9:\n",
    "\n",
    "<center>$bpError^{(l)}_{n{\\times}o_{l}} = \\delta^{(l+1)}_{n{\\times}o_{l+1}}*W^{(l+1)}_{o_{l+1}{\\times}o_{l}}$\n",
    "\n",
    "If we compare equation (10) with equation (6), we can generalize \"error\" there as Backpropagation Error, and the formula for ${\\delta}$ as:\n",
    "\n",
    "$${\\delta}^{(l)}_{n{\\times}o_{l}} = {bpError^{(l)}_{n{\\times}o_{l}}} .* {Y^{(l)}_{n{\\times}o_{l}}} .* (1-Y^{(l)})_{n{\\times}o_{l}}$$\n",
    "\n",
    "Thus,\n",
    "\n",
    "$${\\Delta}W^{(l)}_{{o_{l}{\\times}i_{l}}} = -\\frac{1}{n}\\left({\\delta^{(l)}}^{T}{_{o_{l}{\\times}n}}\\;.*\\;X^{(l)}_{n{\\times}i_{l}}\\right) \\;\\;\\;\\;\\;\\;-------------(11)$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 695,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# IMPLEMENTING BACK-PROPAGATION\n",
    "def backProp(weights, X, Y):\n",
    "    # Forward propagate to find outputs\n",
    "    outputs = forwardProp(X, weights)\n",
    "    \n",
    "    # For the last layer, bpError = error = yPred - Y\n",
    "    bpError = outputs[-1] - Y\n",
    "    \n",
    "    # Back-propagating from the last layer to the first\n",
    "    for l, w in enumerate(reversed(weights)):\n",
    "        \n",
    "        # Find yPred for this layer\n",
    "        yPred = outputs[-l-1]\n",
    "        \n",
    "        # Calculate delta for this layer using bpError from next layer\n",
    "        delta = np.multiply(np.multiply(bpError, yPred), 1-yPred)\n",
    "        \n",
    "        # Find input to the layer, by adding bias to the output of the previous layer\n",
    "        # Take care, l goes from 0 to 1, while the weights are in reverse order\n",
    "        if l==len(weights)-1: # If 1st layer has been reached\n",
    "            xL = addBiasTerms(X)\n",
    "        else:\n",
    "            xL = addBiasTerms(outputs[-l-2])\n",
    "        \n",
    "        # Calculate deltaW for this layer\n",
    "        deltaW = -np.dot(delta.T, xL)/len(Y)\n",
    "        \n",
    "        # Calculate bpError for previous layer to be back-propagated\n",
    "        bpError = np.dot(delta, w)\n",
    "        \n",
    "        # Ignore bias term in bpError\n",
    "        bpError = bpError[:,1:]\n",
    "        \n",
    "        # Change weights of the current layer (W <- W + deltaW)\n",
    "        w += deltaW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 698,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "old weights:\n",
      "1\n",
      "(2, 3)\n",
      "[[-0.87271574  0.35621485  0.95252276]\n",
      " [-0.61981924 -1.49164222  0.55011796]]\n",
      "2\n",
      "(1, 3)\n",
      "[[-1.57656753 -1.10359895 -0.34594249]]\n",
      "old cost:\n",
      "0.107792308277\n"
     ]
    }
   ],
   "source": [
    "# To check with the single back-propagation step done before,\n",
    "# back up the current weights\n",
    "oldWeights = [np.array(w) for w in weights]\n",
    "print(\"old weights:\")\n",
    "for i in range(len(oldWeights)):\n",
    "    print(i+1); print(oldWeights[i].shape); print(oldWeights[i])\n",
    "\n",
    "print(\"old cost:\"); print(nnCost(oldWeights, X, Y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define a function to compute the accuracy of our model, irrespective of the number of neuron in the output layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 699,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Evaluate the accuracy of weights for input X and desired outptut Y\n",
    "def evaluate(weights, X, Y):\n",
    "    yPreds = forwardProp(X, weights)[-1]\n",
    "    # Check if maximum probability is from that neuron corresponding to desired class,\n",
    "    # AND check if that maximum probability is greater than 0.5\n",
    "    yes = sum( int( ( np.argmax(yPreds[i]) == np.argmax(Y[i]) ) and \n",
    "                    ( (yPreds[i][np.argmax(yPreds[i])]>0.5) == (Y[i][np.argmax(Y[i])]>0.5) ) )\n",
    "              for i in range(len(Y)) )\n",
    "    print(str(yes)+\" out of \"+str(len(Y))+\" : \"+str(float(yes/len(Y))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the results of back-propagation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 722,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "950\n",
      "new cost:\n",
      "0.0113971310862\n",
      "new accuracy: \n",
      "4 out of 4 : 1.0\n",
      "[[ 0.03022141]\n",
      " [ 0.13740936]\n",
      " [ 0.13683374]\n",
      " [ 0.7705247 ]]\n"
     ]
    }
   ],
   "source": [
    "# BACK-PROPAGATE, checking old & new weights and costs\n",
    "\n",
    "# Re-initialize to old weights\n",
    "weights = [np.array(w) for w in oldWeights]\n",
    "\n",
    "#print(\"old weights:\")\n",
    "#for i in range(len(weights)):\n",
    "#    print(i+1); print(weights[i].shape); print(weights[i])\n",
    "\n",
    "print(\"old cost: \"); print(nnCost(weights, X, Y))\n",
    "print(\"old accuracy: \"); print(evaluate(weights, X, Y))\n",
    "for i in range(1000):\n",
    "    # Back propagate\n",
    "    backProp(weights, X, Y)\n",
    "\n",
    "    #print(\"new weights:\")\n",
    "    #for i in range(len(weights)):\n",
    "    #    print(i+1); print(weights[i].shape); print(weights[i])\n",
    "    \n",
    "    if i%50==0:\n",
    "        time.sleep(1)\n",
    "        clear_output()\n",
    "        print(i)\n",
    "        print(\"new cost:\"); print(nnCost(weights, X, Y))\n",
    "        print(\"new accuracy: \"); evaluate(weights, X, Y)\n",
    "        print(forwardProp(X, weights)[-1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 718,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Revert back to original weights (if needed)\n",
    "weights = [np.array(w) for w in oldWeights]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training\n",
    "\n",
    "Keep calling backProp() again and again until the cost decreases so much that we reach our desired accuracy.\n",
    "\n",
    "You can observe the cost of the function going down with iterations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problems\n",
    "\n",
    "### - Not reaching desired accuracy fast enough\n",
    "\n",
    "It takes too many iterations of the backProp algorithm for the network to reach the desired output.\n",
    "\n",
    "One of the simplest ways of solving this problem is by adding a Learning Rate (described below) to the back-propagation algorithm.\n",
    "\n",
    "### - Taking too long to compute one iteration\n",
    "\n",
    "Within one iteration, the multiplication and summing operations take too long because there are too many data points feeded into the network.\n",
    "\n",
    "This problem is tackled using Stochastic Gradient Descent (talked about in the next tutorial). The above algorithm is running Batch Gradient Descent. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning Rate\n",
    "\n",
    "Usually, it is desired that we change the amount with which we back propagate, so that we can train our network to reach the desired accuracy faster. So we multiply ${\\Delta}W$ with a factor to control this.\n",
    "\n",
    "$$W \\leftarrow W + \\eta*{\\Delta}W$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $\\eta$ is large, then we take bigger steps to the assumed minimum. If $\\eta$ is small, we take smaller steps.\n",
    "\n",
    "Remember that we are not actually travelling on the gradient, we are only approximating the direction using a ${\\Delta}W$ instead of a ${\\delta}W$. So we don't always point in the direction of the minimum, we could undershoot or overshoot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $\\eta$ is too small, we might take too long to get to the minimum.\n",
    "\n",
    "If $\\eta$ is too big, we might start climbing back up the hill and our cost would keep increasing instead of decreasing!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to ensure that we get the best learning rate is to start at, say, 1,\n",
    "- increase $\\eta$ by 5% if the cost is decreasing\n",
    "- decrease $\\eta$ to 50% if the cost is increasing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Different ways to manipulate learning rate\n",
    "\n",
    "There are various methods available that leverage the variability of learning rate, to produce results that \"converge\" (reach a minimum) faster. The following list includes those with even more complicated methods of trying to converge faster:\n",
    "\n",
    "<center>![Optimizers](images/optimizers.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen, Stochastic Gradient Descent (SGD) itself performs slower than all the other methods, and the one that we are using (Batch Gradient Descent) is even slower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is an implementation of backProp with provision for learning rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 485,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# IMPLEMENTING BACK-PROPAGATION WITH LEARNING RATE\n",
    "# Added eta, the learning rate, as an input\n",
    "def backProp(weights, X, Y, learningRate):\n",
    "    # Forward propagate to find outputs\n",
    "    outputs = forwardProp(X, weights)\n",
    "    \n",
    "    # For the last layer, bpError = error = yPred - Y\n",
    "    bpError = outputs[-1] - Y\n",
    "    \n",
    "    # Back-propagating from the last layer to the first\n",
    "    for l, w in enumerate(reversed(weights)):\n",
    "        \n",
    "        # Find yPred for this layer\n",
    "        yPred = outputs[-l-1]\n",
    "        \n",
    "        # Calculate delta for this layer using bpError from next layer\n",
    "        delta = np.multiply(np.multiply(bpError, yPred), 1-yPred)\n",
    "        \n",
    "        # Find input to the layer, by adding bias to the output of the previous layer\n",
    "        # Take care, l goes from 0 to 1, while the weights are in reverse order\n",
    "        if l==len(weights)-1: # If 1st layer has been reached\n",
    "            xL = addBiasTerms(X)\n",
    "        else:\n",
    "            xL = addBiasTerms(outputs[-l-2])\n",
    "        \n",
    "        # Calculate deltaW for this layer\n",
    "        deltaW = -np.dot(delta.T, xL)/len(Y)\n",
    "        \n",
    "        # Calculate bpError for previous layer to be back-propagated\n",
    "        bpError = np.dot(delta, w)\n",
    "        \n",
    "        # Ignore bias term in bpError\n",
    "        bpError = bpError[:,1:]\n",
    "        \n",
    "        # Change weights of the current layer (W <- W + eta*deltaW)\n",
    "        w += learningRate*deltaW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this back-propagation code, it is better to launch another function that calls it iteratively until we reach the desired accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shall look at training schemes and experiments in the next tutorial."
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}