Backpropagation From Scratch: The Engine, Then the Network

Before building the machine, one question has to be settled, because if the answer went the other way there would be nothing to build.

The perceptron failed on XOR because one line cannot separate the four points. The obvious fix is to stack: run the input through one linear layer, then another. Does that help?

No, and the proof is two lines. A linear layer is $\mathbf{h} = W_1\mathbf{x} + \mathbf{b}_1$. Feed it to another, $\mathbf{y} = W_2\mathbf{h} + \mathbf{b}_2$, and substitute:

The composition is $W\mathbf{x} + \mathbf{b}$ with $W = W_2W_1$ and $\mathbf{b} = W_2\mathbf{b}_1 + \mathbf{b}_2$. A stack of linear layers is a single linear layer. Ten of them, a thousand of them: still one line, still unable to do XOR.

It is worth watching that happen rather than believing it:

import numpy as np
rng = np.random.default_rng(0)
W1, b1 = rng.normal(size=(3, 2)), rng.normal(size=3)
W2, b2 = rng.normal(size=(1, 3)), rng.normal(size=1)
x = rng.normal(size=2)

two_layers = W2 @ (W1 @ x + b1) + b2
one_layer  = (W2 @ W1) @ x + (W2 @ b1 + b2)
print(two_layers[0], one_layer[0], abs(two_layers[0] - one_layer[0]))

-4.612963371048  -4.612963371048  0.00e+00

Not approximately equal. Bit-for-bit identical, because it is the same arithmetic rearranged.

So depth buys nothing on its own. What buys something is putting a nonlinear function between the layers — and that is the entire reason activation functions exist. They are not a biological flourish or a normalisation trick. Without one, the second layer is decoration.

Now the mathematics, and it is one rule you already know applied somewhere slightly unfamiliar.

The single-variable chain rule says that if $L$ depends on $c$ and $c$ depends on $x$, then $\frac{dL}{dx} = \frac{dL}{dc} \cdot \frac{dc}{dx}$. Derivatives multiply along a chain.

The part that matters here is what happens when a variable feeds more than one downstream path. If $x$ influences $L$ through $a$ and also through $b$, the contributions add:

Multiply along a path, sum across paths. That is the whole of backpropagation, and every implementation detail in the rest of this chapter — including the += in the code and the zero_grad() call that trips up everyone who writes their first training loop — is a direct consequence of that second word.

Take a concrete circuit of five operations, with $x = 0.5$ and $y = 1.4$:

Note that $x$ appears three times: in $a$, in $b$, and directly in $L$. Do the backward pass on paper, right to left, starting from $\frac{dL}{dL} = 1$:

Hold on to that number. In a few pages a program is going to produce it without being told any of this.

The insight that makes it programmable: every one of those steps was local. To push a gradient through the multiplication node, you needed the incoming gradient and the two stored input values — nothing about the rest of the circuit. Each operation knows how to differentiate itself.

So make a number that remembers what produced it.

class Value:
    """A number that remembers where it came from."""

    def __init__(self, data, _children=(), _op=""):
        self.data = data
        self.grad = 0.0
        self._backward = lambda: None
        self._prev = set(_children)      
        self._op = _op

Four fields. data is the value. grad accumulates $\frac{\partial L}{\partial \text{self}}$. _prev is the set of Values this one was computed from — the edges of the graph. And _backward is a closure that each operation installs: it knows how to push this node's gradient one step back to its inputs.

Every operator follows the same shape: compute the output, record the parents, install the local rule.

    def __add__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data + other.data, (self, other), "+")

        def _backward():
            self.grad += out.grad       
            other.grad += out.grad      

        out._backward = _backward
        return out

    def __mul__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data * other.data, (self, other), "*")

        def _backward():
            self.grad += other.data * out.grad   
            other.grad += self.data * out.grad   

        out._backward = _backward
        return out

    def tanh(self):
        t = math.tanh(self.data)
        out = Value(t, (self,), "tanh")

        def _backward():
            self.grad += (1 - t * t) * out.grad

        out._backward = _backward
        return out

    def relu(self):
        out = Value(self.data if self.data > 0 else 0.0, (self,), "relu")

        def _backward():
            self.grad += (1.0 if out.data > 0 else 0.0) * out.grad

        out._backward = _backward
        return out

Read the four _backward bodies as a table and the flow patterns from the paper derivation are sitting right there:

Every single one uses += and never =. That is the "sum across paths" rule, encoded. A node that feeds two consumers gets called twice, and the two contributions add up on their own.

Then the driver, which is the only part with any global knowledge:

    def backward(self):
        order, seen = [], set()

        def build(v):
            if v in seen:
                return
            seen.add(v)
            for child in v._prev:
                build(child)
            order.append(v)

        build(self)
        self.grad = 1.0
        for v in reversed(order):       
            v._backward()               

build produces a topological ordering of the graph: every node appears after all of its inputs. Walking that list in reverse guarantees that when you call a node's _backward, its own gradient is already complete — every consumer downstream of it has already contributed. Get the order wrong and you push a half-finished gradient backwards, which produces a wrong answer with no error message.

x = Value(0.5)
y = Value(1.4)
a = x * y
b = x + y
c = a * b
d = c.tanh()
L = d + x
L.backward()
print(x.grad, y.grad)

forward:  a=0.7000  b=1.9000  c=1.3300  d=0.8692  L=1.3692
backward: dL/dd=1.0000  dL/dc=0.2444  dL/da=0.4644  dL/db=0.1711
          dL/dx=1.8212   dL/dy=0.4033

1.8212. The same number, from a program that was told the rule for +, the rule for *, the rule for tanh, and nothing about this circuit.

Two independent checks, because "it matches what I derived" is a weak test when the same person did both.

Numerical differentiation. Nudge the input and measure. The centred difference $\frac{L(x+h) - L(x-h)}{2h}$ estimates the derivative without any calculus at all:

dL/dx:  analytic=1.821202805  numeric=1.821202805  |diff|=1.80e-10
dL/dy:  analytic=0.403269235  numeric=0.403269235  |diff|=7.64e-12

Against PyTorch, which has an industrial autodiff engine written by people who do this for a living:

torch dL/dx=1.821202805316   ours=1.821202805316   |diff|=2.22e-16
torch dL/dy=0.403269234753   ours=0.403269234753   |diff|=1.11e-16

Agreement at $2 \times 10^{-16}$, which is machine epsilon for a 64-bit float: the two engines are performing identical arithmetic. Keep the numerical check in your pocket — it is the tool for debugging a new layer's backward pass, and it is the reason a wrong gradient is findable at all.

The same circuit, different inputs. Set $x = 2$ and $y = -3$, which makes $c = 6$:

x=0.5, y=1.4:  dL/dc = 0.244400     three paths into x:  0.6501 + 0.1711 + 1.0000 = 1.8212
x=2.0, y=-3.0: dL/dc = 0.000025     three paths into x:  0.0001 + -0.0001 + 1.0000 = 0.9999

The gradient crossing the $\tanh$ node fell by a factor of 9,945. Everything upstream of it — in a real network, every layer before it — receives essentially nothing. The two paths through the circuit have gone silent; only the direct connection that skips the $\tanh$ still carries signal.

That is the vanishing gradient problem, in one node. Stack forty layers of $\tanh$ and multiply forty such factors together, and the early layers stop learning entirely. It is also, incidentally, an argument for skip connections that you can see here in miniature: the path that bypassed the nonlinearity is the only one that survived.

Every _backward uses +=. That is correct — it is how paths sum. But it has a consequence that catches everybody: gradients accumulate across calls to backward() too. The engine has no idea that your second call is a new training step rather than another path in the same graph.

So a training loop has to clear them:

for step in range(steps):
    ys = [model(x) for x, _ in DATA]
    loss = sum((yp - yt) ** 2 for yp, (_, yt) in zip(ys, DATA))

    for p in model.parameters():   
        p.grad = 0.0

    loss.backward()
    for p in model.parameters():
        p.data -= lr * p.grad

This is optimizer.zero_grad() in PyTorch, and the usual advice is that forgetting it breaks training. So let us delete those two lines and see how broken it is. Same seeds, same everything, 200 steps of XOR:

At the small learning rates, the buggy version wins every single row. It converges when the correct version stalls.

That is not a fluke and it is worth understanding, because it explains why this bug is so hard to catch. If you never clear the gradient, then at step $k$ the parameter is updated by the sum of every gradient computed so far. On a loss that keeps pointing roughly the same way, that sum grows steadily, and the effect is a learning rate that increases on its own. At $\eta = 0.05$, where the correct algorithm is crawling, the runaway step size looks exactly like a fix.

Then look at the bottom three rows. At $\eta = 0.3$ the same mechanism blows the model apart — loss 8.0 is what a model collapsed to a constant $\pm 1$ scores — half of the 16 that four maximally wrong answers would cost — — while the correct version now converges cleanly.

So the honest statement is not "always call zero_grad or your model will not train". It is: without it you are not running gradient descent any more. You are running something whose step size drifts upward at a rate nobody chose, and it will appear to work, sometimes better than the real thing, right up until it does not — at which point you will blame the learning rate, the initialisation, or the data. This is the shape of the worst bugs in machine learning: they do not crash, they change the algorithm into a different algorithm that occasionally scores better.

With the engine done, a neural network is barely any code. A neuron is a dot product, a bias and an activation; a layer is a list of neurons; a network is a list of layers.

class Neuron:
    def __init__(self, nin):
        self.w = [Value(random.uniform(-1, 1)) for _ in range(nin)]
        self.b = Value(0.0)

    def __call__(self, x):
        act = sum((wi * xi for wi, xi in zip(self.w, x)), self.b)
        return act.tanh()

    def parameters(self):
        return self.w + [self.b]


class Layer:
    def __init__(self, nin, nout):
        self.neurons = [Neuron(nin) for _ in range(nout)]

    def __call__(self, x):
        out = [n(x) for n in self.neurons]
        return out[0] if len(out) == 1 else out

    def parameters(self):
        return [p for n in self.neurons for p in n.parameters()]


class MLP:
    def __init__(self, nin, nouts):
        sizes = [nin] + nouts
        self.layers = [Layer(sizes[i], sizes[i + 1]) for i in range(len(nouts))]

    def __call__(self, x):
        for layer in self.layers:
            x = layer(x)
        return x

    def parameters(self):
        return [p for layer in self.layers for p in layer.parameters()]

There is no backward pass in any of that. Not one line. The Value class already knows how to differentiate whatever these classes happen to build, which is the point of having written it first: an autodiff engine does not know it is being used for a neural network.

Now the problem from Chapter 1. Two inputs, two hidden units, one output, nine parameters:

step   1: loss 4.156690
step  10: loss 4.005572
step  50: loss 3.996708
step 100: loss 3.510700
step 200: loss 0.038597

[0, 0] -> -0.9081  (target -1)  ok
[0, 1] -> +0.8934  (target +1)  ok
[1, 0] -> +0.8906  (target +1)  ok
[1, 1] -> -0.9207  (target -1)  ok

Four out of four. The function that no perceptron can compute — proved in Chapter 1 by four inequalities that demanded $b$ be both positive and negative — is computed by nine numbers found automatically.

The satisfying part is not that it works. It is being able to see how, because with two hidden units the intermediate representation is a point in a plane and you can just print it.

Trained to a loss of 0.001241, here is where each input lands after the hidden layer, and what the output neuron does with it:

Look at the first and fourth rows. The inputs $(0,0)$ and $(1,1)$ are diagonally opposite corners of the square — as far apart as two points in this problem can be — and the hidden layer maps them to $(0.82, -0.85)$ and $(0.84, -0.86)$. Almost the same point. The layer has folded the plane so that the two rejected corners land on top of each other, and once they are in the same place, one line separates them from the other two.

And the output neuron is exactly that line. Its learned parameters are $\mathbf{w} = (-3.1153,\ +3.0893)$, $b = +2.7697$, so its decision boundary is

which is a straight line — a perceptron, the same object from Chapter 1, unchanged. It could not solve XOR then and it cannot now. What changed is that it is no longer looking at the input; it is looking at a space the first layer built for it, in which the problem is linearly separable.

That is what a learned representation is, and it is worth being precise because the phrase gets used loosely for the rest of this course, and for the rest of the field. It is not a compression, a summary, or an embedding in any mystical sense. It is a change of coordinates, learned rather than designed, whose only job is to make the next layer's job easy.

There is a theorem here, and it is usually quoted badly.

Cybenko in 1989 and Hornik in 1991 proved that a feedforward network with a single hidden layer and a suitable activation function can approximate any continuous function on a compact set, to any accuracy you like, given enough hidden units.³⁴ It is a genuine and important result: it says the architecture is not the limitation.

Now read what it omits. It does not say how many units — the bound can be astronomically large. It does not say the weights can be found; it asserts existence, and gradient descent from a random start is not an oracle. And it says nothing about behaviour on data you have not seen, which is the second half of Chapter 6.

The gap between "exists" and "findable" is not academic. Here is the same XOR problem, 50 random initialisations each, 1000 steps, only the hidden layer size changed:

With the minimum viable architecture, one run in four never gets there — it settles into a configuration it cannot descend out of, exactly the local minimum that Chapter 3 showed on a one-dimensional surface. Add one unit and the failures nearly vanish, not because the network became more expressive (two units already suffice — 38 runs prove it) but because extra dimensions give the descent more directions to escape through.

And then eight units does slightly worse than four. At a fixed learning rate and step budget, more capacity is not monotonically better. Anyone who tells you the fix for a stuck network is always a bigger network is extrapolating from the middle of that table.

This is the same lesson as the convergence theorem in Chapter 1, and it will be the same lesson in Chapter 10 about scaling laws, in the form that chapter gives it: a prediction of the loss is not a prediction of the capability you are paying for, and the distance between the two is where the engineering lives.

You now have a network that trains. That is a smaller achievement than it feels, because the network you have trains on four examples and is measured on the same four.

Run the same code on a real dataset and a new set of problems appears, none of which is about gradients. The loss goes down for a while and then stops. Or it goes down on the training data and up on everything else. Or it does not move at all from the first step, and the cause turns out to be the range of the initial random weights. Or one unit's input drifted negative on every example in epoch three and it has been dead ever since, silently, taking a chunk of the model's capacity with it.

These are not exotic failures; they are the normal condition of a network that has just been written, and none of them announces itself. The gradient is correct — you checked it against PyTorch to sixteen decimal places — and the model still does not learn.

Chapter 6 is about that: initialisation, normalisation, overfitting and regularisation, and the diagnostic habit of asking which of those is happening before changing anything. It is the difference between a network that runs and a network that works.

The Value class in this chapter descends directly from Andrej Karpathy's micrograd, and his video The spelled-out intro to neural networks and backpropagation: building micrograd is the best three hours you can spend on this material if you want it explained a second way by someone else. His 2016 post Yes you should understand backprop argues the case for writing one yourself and is assigned reading in Stanford's CS224n. The CS231n notes on backpropagation (cs231n.github.io/optimization-2) are the canonical treatment of the flow patterns tabulated above. For the mathematics as calculus on a graph rather than as neural-network folklore, chapter 5.6 of Mathematics for Machine Learning by Deisenroth, Faisal and Ong is unusually clear; and Baydin, Pearlmutter, Radul and Siskind's survey Automatic Differentiation in Machine Learning: a Survey (arXiv:1502.05767) is the reference for the field as a whole, including the forward/reverse trade-off discussed above.