Forward and Backward Propagation

You have a network with 109,386 weights. It just predicted "3" for an image of an 8, and you would like it to do better. Which of those 109,386 numbers should go up, which should go down, and by how much?

There is an obvious way to find out. Take the first weight, add a tiny amount to it, run the whole network forward again, and see whether the error got smaller. If it did, that weight should go up. Then put it back and try the second weight. Then the third.

Cost the plan. One forward pass through that network is roughly 220,000 arithmetic operations (about 110,000 multiplications and as many additions). You need one per weight, so 109,386 × 220,000 ≈ 24 billion operations — to compute the gradient for a single training image. Multiply by 60,000 images, multiply by 30 passes over the data, and you are asking for about $4 \times 10^{16}$ operations. Weeks of GPU time, for a model that trains in ninety seconds.

Backpropagation computes the exact same 109,386 numbers in about the cost of two forward passes. Not two per weight. Two, total. That factor of fifty thousand is not an optimisation; it is the difference between deep learning existing and not existing. What follows is how it works, worked through with real numbers small enough to check by hand.

Forward propagation: what the network computes

Before you can assign blame you need to know what happened. The forward pass runs data through the layers, and it is worth writing out precisely because the backward pass is its mirror image.

For layer $l$ with weight matrix $W^{(l)}$, bias vector $b^{(l)}$ and activation function $f$:

with $a^{(0)} = x$, the input. Two names to keep straight, because every backprop derivation depends on the distinction: $z$ is the pre-activation (the raw weighted sum) and $a$ is the activation (after the non-linearity).

In code, with a batch of examples stacked as rows:

1import numpy as np23def relu(z):4    return np.maximum(0.0, z)56def forward(x, params):7    """x: (batch, 784). Returns prediction plus every intermediate value."""8    cache = {"a0": x}9    z1 = x  @ params["W1"] + params["b1"]      # (batch, 128)10    a1 = relu(z1)11    z2 = a1 @ params["W2"] + params["b2"]      # (batch, 64)12    a2 = relu(z2)13    z3 = a2 @ params["W3"] + params["b3"]      # (batch, 10) raw class scores14    cache.update(z1=z1, a1=a1, z2=z2, a2=a2, z3=z3)15    return z3, cache

Note the cache. Every intermediate value is kept, not thrown away. This is not sloppiness — the backward pass needs $a^{(l-1)}$ to compute the gradient of $W^{(l)}$, and it needs $z^{(l)}$ to compute the derivative of the activation. That cache is why training uses several times more memory than inference does, and it is why "reduce the batch size" is the first thing anyone suggests when you hit an out-of-memory error.

The forward pass is not just a prediction. It is also the recording of everything the backward pass will need.

The chain rule, which is the whole trick

A network is a composition of functions. The loss depends on the output, the output depends on the last layer's weights and on the previous activations, which depend on earlier weights, and so on back to the input. Calculus has exactly one tool for differentiating a composition, and it is the chain rule:

Read it as a chain of "how much does this affect the thing in front of it". A change in $w$ changes $z$ by some rate; that change in $z$ changes $a$ by some rate; that changes the loss by some rate. Multiply the rates and you get the end-to-end rate.

The insight that makes it fast: the middle terms are shared. Every weight in layer 2 needs "how much does $L$ change when layer 2's output changes". Compute that quantity once, store it, and reuse it for all of layer 2's weights — and again, propagated one step further back, for all of layer 1's weights. That stored quantity is conventionally called $\delta$, and computing it layer by layer from the output backwards is the entire algorithm.

A complete numerical example

Take the smallest network that still has everything: one input, one hidden unit with ReLU, one linear output, squared-error loss.

x ---[w1, b1]---> z1 ---ReLU---> a1 ---[w2, b2]---> yhat ---> L = (yhat - y)^2x  = 1.0     w1 = 0.5    b1 = 0.1y  = 1.0     w2 = -0.3   b2 = 0.2

Forward

Backward

Start at the loss and walk backwards. Each row uses the row above it.

Now take one gradient-descent step with learning rate 0.1, moving each parameter against its gradient:

Run the forward pass again with the new values: $z_1 = 0.4824$, $a_1 = 0.4824$, $\hat{y} = 0.308$, and $L = 0.4788$. The loss fell from 0.9604 to 0.4788 in one step. That is the entire training algorithm, and everything else in deep learning is this loop with more parameters and better bookkeeping.

Two details in that table are worth pausing on. First, $\partial L/\partial w_2 = -1.96 \times 0.6$ — the gradient of a weight is the incoming activation times the error signal arriving from above. If the incoming activation is zero, the gradient is zero and that weight does not move at all. Second, $\text{ReLU}'(z_1)$ acted as a gate: it was 1 because $z_1$ was positive. Had $z_1$ been negative, the derivative would be 0 and the entire error signal would have been blocked from reaching $w_1$ and $b_1$.

The general algorithm, in matrix form

Scaled up to layers of many units, with $\odot$ meaning element-by-element multiplication:

Written as NumPy, for the three-layer network from earlier:

1def backward(y_true, logits, cache, params):2    batch = y_true.shape[0]3    grads = {}45    # Softmax + cross-entropy have a famously clean combined gradient:6    # d(loss)/d(logits) = (predicted_probs - one_hot_true) / batch7    probs = np.exp(logits - logits.max(axis=1, keepdims=True))8    probs /= probs.sum(axis=1, keepdims=True)9    delta3 = (probs - y_true) / batch                       # (batch, 10)1011    grads["W3"] = cache["a2"].T @ delta312    grads["b3"] = delta3.sum(axis=0)1314    delta2 = (delta3 @ params["W3"].T) * (cache["z2"] > 0)  # ReLU' is a 0/1 mask15    grads["W2"] = cache["a1"].T @ delta216    grads["b2"] = delta2.sum(axis=0)1718    delta1 = (delta2 @ params["W2"].T) * (cache["z1"] > 0)19    grads["W1"] = cache["a0"].T @ delta120    grads["b1"] = delta1.sum(axis=0)21    return grads

Count the matrix multiplications: three going forward, five coming back (a weight gradient for each of the three layers, plus propagating $\delta$ down through the top two — there is no need to send it into the input). Backprop costs roughly twice a forward pass — which is where the earlier claim comes from.

Why the division by batch size

Dividing $\delta^{(3)}$ by the batch size makes the gradient an average over the batch rather than a sum. Skip it and your effective learning rate scales with batch size: switch from batches of 32 to batches of 256 and your updates become eight times larger, which usually shows up as a loss that explodes to NaN within a few steps. It is a one-character bug with a dramatic symptom.

What autograd does instead

You will almost never write the code above. Both major frameworks record the operations performed during the forward pass and differentiate that recording automatically.

1import torch23x = torch.tensor([1.0])4w1 = torch.tensor([0.5],  requires_grad=True)5b1 = torch.tensor([0.1],  requires_grad=True)6w2 = torch.tensor([-0.3], requires_grad=True)7b2 = torch.tensor([0.2],  requires_grad=True)8y  = torch.tensor([1.0])910z1   = w1 * x + b111a1   = torch.relu(z1)12yhat = w2 * a1 + b213loss = (yhat - y) ** 21415loss.backward()          # walks the recorded graph backwards16print(w1.grad, b1.grad)  # tensor([0.5880]) tensor([0.5880])17print(w2.grad, b2.grad)  # tensor([-1.1760]) tensor([-1.9600])

Identical to the hand computation. TensorFlow does the same thing with an explicit recorder:

1import tensorflow as tf23w1 = tf.Variable([0.5]); b1 = tf.Variable([0.1])4w2 = tf.Variable([-0.3]); b2 = tf.Variable([0.2])5x = tf.constant([1.0]); y = tf.constant([1.0])67with tf.GradientTape() as tape:8    a1   = tf.nn.relu(w1 * x + b1)9    yhat = w2 * a1 + b210    loss = tf.square(yhat - y)1112grads = tape.gradient(loss, [w1, b1, w2, b2])13print([g.numpy() for g in grads])   # [0.588, 0.588, -1.176, -1.96]

The mistake everyone makes once. PyTorch accumulates gradients into .grad rather than replacing them. Forget optimizer.zero_grad() at the top of your training loop and step $n$ uses the sum of the gradients from steps 1 through $n$. The loss will not error; it will simply behave as though the learning rate is growing without bound, and diverge.

Checking that a hand-written gradient is right

If you ever do implement a custom backward pass, verify it numerically before trusting it. The finite-difference method that was hopeless as a training algorithm is perfectly good as a test, because you only run it on a handful of parameters once.

1def gradient_check(f, params, key, index, eps=1e-5):2    """Compare an analytic gradient against a central-difference estimate."""3    original = params[key].flat[index]45    params[key].flat[index] = original + eps6    loss_plus = f(params)7    params[key].flat[index] = original - eps8    loss_minus = f(params)9    params[key].flat[index] = original1011    numeric = (loss_plus - loss_minus) / (2 * eps)12    return numeric1314# Relative error below 1e-7 means correct; above 1e-4 means a real bug.15rel_err = abs(numeric - analytic) / max(1e-12, abs(numeric) + abs(analytic))

Use the central difference $(f(x+\epsilon) - f(x-\epsilon)) / 2\epsilon$, not the one-sided version — its error shrinks as $\epsilon^2$ rather than $\epsilon$, which matters when you are fighting floating-point noise. And do not check ReLU networks at points near zero, where the derivative genuinely does not exist and the numeric estimate will disagree for legitimate reasons.

Why gradients die on the way back

Look again at the propagation rule: $\delta^{(l-1)} = (\delta^{(l)} W^{\top}) \odot f'(z^{(l-1)})$. Going back one layer multiplies the error signal by weights and by the activation's derivative. Going back twenty layers multiplies by twenty such factors.

If each factor is typically 0.25 — which is exactly the maximum derivative of the sigmoid function — then after 20 layers the signal has been multiplied by $0.25^{20} \approx 10^{-12}$. The early layers receive a gradient indistinguishable from zero and stop learning entirely. This is the vanishing gradient problem, and it is the direct reason deep sigmoid networks failed to train through the 1990s.

If each factor is typically 1.5, then $1.5^{20} \approx 3300$ and the early layers receive a gradient thousands of times too large, producing weight updates that overshoot into numerical overflow. That is the exploding gradient problem, and it announces itself as a loss of nan.

ReLU's derivative is exactly 1 for positive inputs. That single fact — no shrinking factor on the way back — is why replacing sigmoid with ReLU made deep networks trainable.

Using this when you are debugging at 2am

The reason to understand backprop, given that autograd writes it for you, is that gradients are the thing you inspect when training misbehaves. Print them.

1loss.backward()2for name, p in model.named_parameters():3    if p.grad is not None:4        print(f"{name:20s} grad_norm={p.grad.norm().item():.3e}")

Reading that output tells you almost immediately what class of problem you have. Norms around $10^{-8}$ in the early layers and $10^{-2}$ in the late layers means vanishing gradients — change activations or initialisation. Norms of $10^{3}$ or higher means exploding gradients — clip them. A gradient norm of exactly zero for a whole layer means that layer is either frozen, disconnected from the loss, or has been killed by dead ReLUs. And gradients that are None mean the tensor never entered the computation graph at all, usually because someone converted it to NumPy and back in the middle of the forward pass, silently severing the chain.

None of those diagnoses are available to you if backprop is a black box. They are all obvious the moment you picture the error signal being multiplied, layer by layer, on its way from the loss back to the input.