Backpropagation by hand

I said a few years ago that after Andrew Ng’s class I never wanted to derive backprop by hand again. In 2011 I lost an entire Saturday on the neural network assignment in Octave to an off-by-one in the bias term, and when it finally worked I promised myself that from then on I’d let a library work out the derivatives. Well. The neural nets section of the data science certificate I’ve been doing wants it done from scratch, in R this time, and it turns out that seven years later I’d forgotten nearly all of it. So here it is again, written down properly so I don’t have to do it a third time.

The whole thing is the chain rule. Everybody says that, and it’s true, but it never helped me much, because the hard part isn’t the calculus, it’s the bookkeeping: which number gets multiplied by which, which index is the row, and whether you use the value before or after the activation. My math degree covered the chain rule in the first semester. Nobody covered keeping track of 40 partial derivatives in a for loop at 11pm.

So let’s keep it tiny. Two inputs, two hidden units, one output, sigmoid everywhere, one training example. Small enough that you can do every step on paper and check it against R.

The setup

sigmoid <- function(z) 1 / (1 + exp(-z))

x  <- c(1.0, 0.5)    # one training example
y  <- 1              # its target

W1 <- matrix(c(0.4, -0.6,
               0.3,  0.8), nrow = 2, byrow = TRUE)  # row j = weights into hidden unit j
W2 <- c(0.7, -0.5)                                  # hidden -> output

I’ve left out the biases to keep the arithmetic short. They behave exactly like a weight on an input that’s always 1, so once you can do weights you can do biases. (You do need them for anything real, more on that below.)

Forward

z1   <- as.vector(W1 %*% x)   # 0.10, 0.70
h    <- sigmoid(z1)           # 0.525, 0.668
z2   <- sum(W2 * h)           # 0.0334
o    <- sigmoid(z2)           # 0.508
loss <- 0.5 * (o - y)^2       # 0.121

Hidden unit 1 gets 0.4 × 1.0 + (−0.6) × 0.5 = 0.1, and unit 2 gets 0.3 × 1.0 + 0.8 × 0.5 = 0.7. Squash both and you get 0.525 and 0.668. The output unit gets 0.7 × 0.525 − 0.5 × 0.668 = 0.0334, which is close to zero, so the sigmoid gives almost exactly 0.5. We wanted 1, so the loss is ½(0.508 − 1)² = 0.121. A net with random-ish weights is basically shrugging, which is what you’d expect.

Backward

The question backprop answers is: if I nudge this one weight, how much does the loss move? You want that number, the partial derivative, for every weight, and the trick is that you can get all of them in one sweep from right to left, reusing what you computed for the layer after.

Start at the output. The loss is ½(o − y)², so ∂loss/∂o = o − y = −0.492. The output is sigmoid(z2), and the sigmoid has the nice property that its derivative is s(1 − s), so ∂o/∂z2 = 0.508 × 0.492 = 0.250. Multiply them and you get what everyone calls the delta of the output unit, how much the loss cares about the input to that unit:

delta2 <- (o - y) * o * (1 - o)    # -0.123

Now the hidden-to-output weights. z2 = W2[1] × h1 + W2[2] × h2, so ∂z2/∂W2[1] is just h1. One more multiply:

gW2 <- delta2 * h                  # -0.0645, -0.0821

Both gradients are negative, which means increasing either weight would reduce the loss. That makes sense: both hidden units are positive and we want the output to go up.

The hidden layer is where I always got lost, so slowly. The loss reaches hidden unit 1 only through z2, and z2 depends on h1 through the weight W2[1] = 0.7. So ∂loss/∂h1 = delta2 × 0.7. Then h1 = sigmoid(z1[1]), so multiply by h1(1 − h1) to get to z1[1]. That’s the hidden delta:

delta1 <- delta2 * W2 * h * (1 - h)   # -0.0215, 0.0136
gW1    <- outer(delta1, x)            # 2x2, row j = gradients into hidden unit j

Unit 1’s delta is −0.123 × 0.7 × 0.525 × 0.475 = −0.0215. Unit 2’s is positive, because its outgoing weight is negative: making h2 bigger makes the output smaller, which is the wrong direction. And the gradient for a weight into a hidden unit is that unit’s delta times whatever came in on that wire, so the weight from x1 to h1 gets −0.0215 × 1.0, and from x2 to h1 gets −0.0215 × 0.5. That’s the outer().

That’s the entire algorithm. With more layers you just repeat the hidden step: a unit’s delta is the weighted sum of the deltas of the units it feeds, times its own sigmoid derivative. In matrix form it’s t(W) %*% delta * h * (1 - h), which is the line I’d seen in a hundred tutorials and never really understood until I’d done the scalar version on paper.

Did I get it right?

Don’t trust your derivatives, check them numerically. Nudge one weight up and down by a tiny epsilon, compute the loss both times, and the slope between them should match your gradient:

loss_fn <- function(W1, W2) {
  h <- sigmoid(as.vector(W1 %*% x))
  o <- sigmoid(sum(W2 * h))
  0.5 * (o - y)^2
}

eps <- 1e-5
Wp <- W1; Wp[1, 1] <- Wp[1, 1] + eps
Wm <- W1; Wm[1, 1] <- Wm[1, 1] - eps
(loss_fn(Wp, W2) - loss_fn(Wm, W2)) / (2 * eps)   # -0.02145

−0.02145, against −0.021450 from backprop. They agree to about ten decimal places. If they hadn’t, I’d know there was a bug, and more usefully, because you can check each weight separately, I’d know which layer it was in. Ng had us run a gradient check, and in 2011 it told me my gradients were wrong and nothing about where, because I only ever ran it on the whole network at once. Checking one weight at a time is the part I missed, and it would have found that bias term in minutes.

The bugs I hit this time, in order: using z instead of h in the sigmoid derivative (so z1 * (1 - z1), which is nonsense but runs fine and produces numbers); forgetting the transpose and getting a 2×2 where I wanted a vector, which R happily recycled instead of complaining; and a sign error in the update. Every one of them gave a loss that went down a bit and then stalled, which looks exactly like “the learning rate is wrong”. The gradient check caught all three in a few minutes.

One step

Gradient descent is then one line per weight:

lr <- 0.5
W2 <- W2 - lr * gW2
W1 <- W1 - lr * gW1

Run the forward pass again and the output goes from 0.508 to 0.520, and the loss from 0.121 to 0.115. A small step in the right direction, which is all one step ever is.

Making it learn something

One example is a toy, so the next thing I did was train it on XOR, the four-row truth table that a single-layer perceptron can’t learn (that’s the Minsky and Papert result that supposedly froze the field for a decade, though I’ve read enough arguments about that history to not want to repeat it with any confidence). This is where you need biases. Without them every hidden unit’s boundary has to pass through the origin, and you can’t carve up XOR that way.

I looped over all four rows, summed the gradients, and took a step per epoch, with weights starting uniformly between −1 and 1. With two hidden units, most runs got the loss under 0.01 in roughly 600 to 1,600 epochs. But a fair number of runs, somewhere between a quarter and a third of the random starts I tried, just sat there, outputting 0.5 or so for two of the rows forever. That’s the network finding a local minimum, and it’s a real thing even for a problem this small. Bump it to three hidden units and every start I tried got there. I’m not sure how much intuition from XOR carries over to real networks, but it was the first time “more parameters makes the loss surface easier” felt like something I’d seen rather than something I’d read.

The R is about 40 lines in total. It’s slow, R loops always are, and nobody should train anything real this way. I’m going back to letting a library do the derivatives, again. But I’ve now actually watched every number move, which I couldn’t honestly claim in 2011.