Neural Networks: The Backward Pass

Imagine someone tells you to minimise some random function \(g(x)\) you know nothing about. The only thing you can ask is, for an \(x\) you pick, what is \(g(x)\)? Most people will start with some random \(x\) and nudge it in different directions to see which direction produces a smaller \(g(x)\) and hope that if they keep nudging \(x\) in that direction, \(g(x)\) will keep going down.

Russell and Norvig explain a similar algorithm in their book Artificial Intelligence–A Modern Approach with the statement, "trying to find the top of Mount Everest in a thick fog while suffering from amnesia" [1]. Except, in our case, we are trying to find the base-camp in a thick fog with amnesia, i. e., going down. We can't see too far, we don't remember where we came from, but we have a sense of downwards direction. We just know that if we keep stepping (making nudges) in the downwards direction, we shall eventually reach a local minimum. This idea is the very key to finding the correct weights and biases. Calculus students know this rate of change as the derivative and most algebra students know it as the slope. We never get into algorithms to minimise the simple derivatives or slopes because we can just set them equal to zero and solve. But, in order to illustrate the algorithm, let's consider a simple polynomial as an example.

Take this function \(g(x) = x^2\), we know that the derivative is \(2x\) and \(2x = 0 \iff x = 0\) hence \(x = 0 \Rightarrow g(x) = 0\). We can do this because \(g\) has such a nice and simple definition. What \(g\) it did not have this nice and simple definition like our network \(f\). We could do what we just talked about, start at some random \(x_0\) and make small \(\eta\) sized nudges to it using the negative derivative: \(x_{n+1} = x_n - \eta g'(x_n)\) as shown in figure 2. The small nudge to the initial \(x_0\) gives us \(x_1\) and then we calculate \(-\eta g'(x_1)\) to figure out the nudge we want to make to \(x_1\) and the series continues till we find some \(x_i\) for which \(g(x_i)\) is small enough.

For \(\vec{p} = [a, b]\), take \(g(\vec{p}) = g(a, b) = a^2 + b^2\), imagine this to be our error function where \(a\) is analogous to \(\mathcal{W}\) and \(b\) is analogous to \(\mathbf{b}\). Let's write the gradient descent algorithm,

1. Pick a random \(\vec{p}_{0}\).
2. Let \(\vec{p}_{i+1} = \vec{p}_{i} - \eta \nabla g(\vec{p}_i)\).
3. Repeat step 2 for all \(i\).

Can you implement the gradient descent for the function \(g\) in your favourite programming language? Here, I'll do it in Python, but first let's write the gradient of \(g\).

\[ \nabla g(\vec{p}_i) = \Bigg[\frac{\partial g(\vec{p})}{\partial p_1}, \frac{\partial g(\vec{p})}{\partial p_2}\Bigg] = \Bigg[\frac{\partial g(a,b)}{\partial a}, \frac{\partial g(a,b)}{\partial b}\Bigg] = \Bigg[\frac{\partial (a^2 + b^2)}{\partial a}, \frac{\partial (a^2 + b^2)}{\partial b}\Bigg] = [2a, 2b] \]

def g(a, b, get_gr=False):
  val, nabla_g = a**2 + b**2, [2 * a, 2 * b]
  error = (0 - a)**2 + (0 - b)**2
  return nabla_g, error if get_gr else val

def GD(a=10.0, b=10.0, eta=0.3, print_steps=True):  # Gradient Descent
  step, step_error_sum = 0, 0
  for step in range(7):
    gradient, step_error = g(a, b, get_gr=True)
    step_error_sum += step_error
    a = a - (eta * gradient[0])
    b = b - (eta * gradient[1])
    if print_steps:
      to_print = 'Step: {:>2}, (a, b): ({:>1.2}, {:>1.2}) Error: {:>3.5}'
      print(to_print.format(step, a, b, step_error))
  return step_error_sum / step

GD()

We know that \(g(a, b) = 0\) for \((a,b)=(0,0)\). As we can see in the printed output (also in figure 3), at each gradient descent step, the parameters \((a, b)\) get closer to \((0,0)\) alongside the error approaching \(0\). This means that we converge towards \((a,b)\) that minimise \(g(a,b)\) while descending on the gradient.

Step:  0, (a, b): (4.0, 4.0)     Error: 200.0
Step:  1, (a, b): (1.6, 1.6)     Error: 32.0
Step:  2, (a, b): (0.64, 0.64)   Error: 5.12
Step:  3, (a, b): (0.26, 0.26)   Error: 0.8192
Step:  4, (a, b): (0.1, 0.1)     Error: 0.13107
Step:  5, (a, b): (0.041, 0.041) Error: 0.020972
Step:  6, (a, b): (0.016, 0.016) Error: 0.0033554

Let's start with the \(\frac{\partial E}{\partial W^{(L)}}\).

We have a composite function: \((\vec{y} - \vec{a}^{(L)})^2\). We'll need the chain rule to move on. The application of chain rule here needs a little care. You might already be able to take polynomial derivatives using the chain rule without thinking much of it; but, in order to truly understand its application here, the reader should not only be able to apply the chain rule to the polynomials but also be aware of the steps they are taking in terms of its notation. Here is a nice example that walks us through the chain rule in Leibniz's notation⁶, which is what we'll be using. We state the chain rule.

\[ \frac{\operatorname{d} f \circ g (x)}{\operatorname{d} x} = \frac{\operatorname{d} f(g(x))}{\operatorname{d} x} = \frac{\operatorname{d} f(g(x))}{\operatorname{d} g(x)} \times \frac{\operatorname{d} g(x)}{\operatorname{d} x} = \frac{\operatorname{d} f}{\operatorname{d} g} \times \frac{\operatorname{d} g}{\operatorname{d} x} \]

Following is an intuitive way to think about it [3].

\[\begin{array}{c} \text{The change in }f\circ g\text{ caused by}\\ \text{a small unit change in }x \end{array}=\begin{array}{c} \text{The change in }f\text{ caused by}\\ \text{a small unit change in }g \end{array}\times\begin{array}{c} \text{The change in }g\text{ caused by}\\ \text{a small unit change in }x. \end{array}\]

Therefore,

Similarly for \(\frac{\partial E}{\partial \vec{b}^{(L)}}\), the same derivation will apply except at the end we'll have,

\[ \frac{\partial (\vec{b}^{(L)})}{\partial \vec{b}^{(L)}} = 1 \Rightarrow \frac{\partial E}{\partial \vec{b}^{(L)}} = -\frac{\partial (\vec{y} - \vec{a}^{(L)})^2}{\partial (\vec{y}-\vec{a}^{(L)})} \frac{\partial (\sigma(\vec{z}^{(L)}))}{\partial \vec{z}^{(L)}} \frac{\partial \vec{z}^{(L)}}{\partial \vec{b}^{(L)}} = -2(\vec{y}-\vec{a}^{(L)})\sigma'(\vec{z}^{(L)})(1) \]

You'll often see the term \(-2(\vec{y} - \vec{a}^{(L)}) \sigma'(\vec{z}^{(L)})\) written as \(\delta^{(L)}\). This is known as the delta rule or error in the layer \(l\). We have the three grand equations (each arranged correctly per dimension) as a result⁷:

But what about all the rest of the partial derivatives for layers \(l < L\)? Well, we just keep applying the chain rule in the above derivation instead of stopping at \(\partial (W^{L}\vec{a}^{(L-1)})\). Let's do it for \(l = L-1, \frac{\partial E}{\partial W^{(L-1)}}\).

We now generalise (correcting dimensions) through induction, for all \(l = L - i\),

The only thing keeping us from running the implementation of the stochastic gradient descent was the lack of a way to calculate the gradient of the error function \(E\). That is no longer true! We have a plan of attack. We do a forward pass using some \(\vec{x} \in \mathcal{B}_i\) and then we use equations (1), (2) and (3) to calculate the layer \(L\) partial derivative in the gradient. From there, we use the \(\delta^{(L)}\) to recursively keep calculating \(\delta^{(l)}\) to be used in equations (4), (5) and (6) in order to calculate the rest of the partial derivatives backwards. This backwards pass is why we call this algorithm Backpropagation. Let's implement backpropagation for \(\vec{x} \in \mathcal{B}_i\).

def backpropagation(self, x, y):
  outputs, activations = self.forward_pass(x)
  gradient = self.backward_pass(outputs, activations, y)
  return gradient, self.error(y, activations[-1])

def backward_pass(self, outputs, activations, y):
  gradient_W = np.empty(shape=self.L - 1, dtype=np.object)
  gradient_b = np.empty(shape=self.L - 1, dtype=np.object)
  z, a = outputs[-1], activations[-1]  # z^L, a^L
  delta = -2 * (y - a) * self.sigmoid(z, derivative=True)  # delta^L eq 2
  delta = delta.reshape((1, len(delta)))
  for l in range(self.L - 1, 0, -1):
    a_prev = activations[l - 1]
    a_prev = a_prev.reshape((len(a_prev), 1)).T
    pC_w = np.dot(delta.T, a_prev)  # eq 1 or 4
    pC_b = delta.flatten()  # eq 3 or 6
    gradient_W[l - 1], gradient_b[l - 1] = pC_w, pC_b
    if l == 1:
      break
    z, a = outputs[l - 2], activations[l - 1]
    delta = np.dot(delta, self.W[l - 1]) * self.sigmoid(z, derivative=True)  # eq 5
  gradient = np.empty(shape=2, dtype=np.object)
  gradient[0], gradient[1] = gradient_W, gradient_b
  return gradient

Let's first test network.py by building and training the network in figure 4. We produce the training and testing examples (skipping validation). The training examples are what we pass as X, y. The testing examples are the ones we use to test the accuracy of a trained network. We don't train the network on this set for the sake of an unbiased measure of performance.

\[ \text{Accuracy} = ACC = \frac{\text{Number of Correctly Classified Examples in Test Set}}{\text{Number of Incorrectly Classified Examples in Test Set}} \]

from network import Network
import numpy as np

N, n = 700, 10
X = np.array([list(map(int, '{0:010b}'.format(e))) for e in range(2**n)])
y = np.array([int(e % 2 == 0) for e in range(2**n)])
X, y, X_test, y_test = X[:N], y[:N], X[N:], y[N:]

Since we have the data, we build and train the figure 4 network.

np.random.seed(0)
net = Network(X, y, structure=[10, 3, 2], epochs=1000, bt_size=8, eta=0.3)

# train
net.tango()

# test
predictions = np.array([net.predict(x.flatten()) for x in X_test])
acc = np.sum(predictions == y_test) / len(y_test)
print('Network Accuracy: {}'.format(acc))

We get the output:

* Epoch:    0, Error: 0.48236
* Epoch:    1, Error: 0.3234
* Epoch:    2, Error: 0.16235
* Epoch:    3, Error: 0.075001
* Epoch:    4, Error: 0.040232
* Epoch:    5, Error: 0.025256
* Epoch:    6, Error: 0.017793
* Epoch:    7, Error: 0.013523
* Epoch:    8, Error: 0.010801
* Epoch:    9, Error: 0.0089463
[snip]
* Epoch:  998, Error: 0.00004031
* Epoch:  999, Error: 0.00004026
Network Accuracy: 1.0

Usually if you see an accuracy of 1, that should make you suspicious of having over-fit. However, in this case, the parity classification of a binary number is pretty easy, i. e., we can just look at the least significant bit. In-fact, let's plot our two weight matrices to see if that is what the network learned.

Looking at figure 5 and sure enough, training assigned the heaviest weight values to the least significant bit. This pattern shows that sometimes by looking at the weights, even humans can learn about new patterns or simply perform dimensionality reduction of the input \(\vec{x}\) (feature selection).

The benchmark problem that all neural network tutorials solve is the classification of the handwritten digits. At the time of writing this article, the data is available for free at http://yann.lecun.com/exdb/mnist. If you are anything like me, you'll open that link and say to yourself, "Okay, wow, now what?". To deal with that feeling, I wrote a little downloader class in Python that will download, read and reshape the data from that site into the format we want.

from mnsit_handwritten_digits import MNSIT
from network import Network
import numpy as np

# MNSIT Data
mn = MNSIT(path='../mnsit_data/')
mn.plot_image(999, source='training')
train_X, test_X, train_y, test_y = mn.get_data()

You'll need Matplotlib for plotting. The line mn.plot_image(999, source='training') will plot the \(999^{th}\) image (figure 6) in the training set.

MNSIT images are normalised to dimensions \(28 \times 28 = 784\). We'll create a network with 3 layers:

1. \(28^2 = 784\) neurons in the input layer to hold the activations caused by the pixel values in the flattened image.
2. 32 neurons in the hidden layer, for no particular reason.
3. 10 neurons in the output layer to hold the probabilities of each class 0 through 9.

# Network
np.random.seed(0)
net = Network(train_X,
              train_y,
              structure=[784, 32, 10],
              epochs=300,
              bt_size=256)
# net.load_weights_biases(path='weights_biases.npy')

# train
net.tango()
# net.save_weights_biases(path='weights_biases.npy')

# test
predictions = np.array([net.predict(x.flatten()) for x in test_X])
acc = np.sum(predictions == test_y) / len(test_y)
print('Network Accuracy: {}'.format(acc))

Output after running,

* Epoch:    0, Error: 1.1664
* Epoch:    1, Error: 0.837
* Epoch:    2, Error: 0.74615
* Epoch:    3, Error: 0.66387
* Epoch:    4, Error: 0.61246
* Epoch:    5, Error: 0.57026
* Epoch:    6, Error: 0.54528
* Epoch:    7, Error: 0.52077
* Epoch:    8, Error: 0.49434
* Epoch:    9, Error: 0.47249
* Epoch:   10, Error: 0.44921
[snip]
* Epoch:   49, Error: 0.21924
[snip]
Network Accuracy: 0.9259

The training time here will be longer depending on your machine. With the network above, I was able to achieve an accuracy of \(93\%\). You may load my converged set of \((\mathcal{W}, \mathbf{b})\) by uncommenting the line after the network instantiation. Make sure you give the correct path to the file weights_biases.npy (and you should probably comment net.tango()).