Backpropagation

Problem

Training needs ∂L/∂θ for every parameter θ, at every step. A modern network has millions or billions of parameters. Estimating each derivative by nudging that parameter and re-running the network would cost one forward pass per parameter: a billion forward passes for one step.

Backpropagation gets all of them for the price of about one extra forward pass. The chain rule page showed why: reverse mode differentiation sends one gradient backwards from the scalar loss, and each operation turns the gradient of its output into the gradients of its inputs. This page is about doing that for a real network: with a batch, a softmax output and a parameter update; with dozens of layers, where the backward product can shrink to nothing or blow up; and with a way to know your gradients are right.

Intuition

The forward pass computes the network's prediction and the loss, layer by layer, and keeps the intermediate values. The backward pass walks the same layers in reverse. At each step it holds one quantity, the gradient of the loss with respect to that layer's output, which says "if this value went up a little, the loss would change by this much". Each layer uses its saved inputs to turn that into two things: the gradient for its own parameters (which the optimizer uses) and the gradient for its input (which it hands to the layer below).

Credit assignment is the point. The output layer knows directly how wrong it was. A hidden unit deep in the network only learns how much it contributed through the weights that connect it to the output, and a ReLU unit that was switched off contributed nothing, so it gets no gradient.

Forward values in green, the gradient of the loss with respect to each value in amber. The dashed hidden unit had a negative input, so ReLU blocked it in both directions: its output was 0 and its gradient is 0, and the faint weights leaving it get no update.

Mechanics

The algorithm

Write the network as a chain of operations v₁ = f₁(x), v₂ = f₂(v₁), …, L = f_K(v_K−1), each with its own parameters. Write v̄ for ∂L/∂v, a tensor with the same shape as v.

1. Forward. Compute each v_k in order, saving whatever each operation's backward rule will need (usually its inputs).
2. Seed. Set L̄ = 1.
3. Backward. For k = K down to 1, apply f_k's local rule: from v̄_k and the saved values, compute v̄_k−1 and the gradient of f_k's parameters.
4. Accumulate. If a value feeds several operations (a residual connection, a shared weight, an embedding used twice), its gradient is the sum of the gradients arriving from each use.

On a general graph, "in reverse" means reverse topological order: a node's backward rule runs only after every operation that used its output has sent its gradient back. That is the whole algorithm. Everything else is the catalogue of local rules and engineering.

With a batch

Real training uses a batch of N examples at once. For a linear layer Y = XWᵀ + b with X of shape (N, d_in) and W of shape (d_out, d_in), the rules become matrix products:

X̄ = Ȳ W,    W̄ = Ȳᵀ X,    b̄ = Σ_rows Ȳ

Check the shapes and you can rederive them: W̄ must be (d_out, d_in), and Ȳᵀ X is the only product of the available tensors with that shape. The product Ȳᵀ X sums over the batch dimension, so a parameter's gradient adds up every example's contribution, and the bias, broadcast across rows in the forward pass, gets a sum over rows in the backward pass. If the loss is a mean over the batch (PyTorch's default reduction="mean"), Ȳ already carries the 1/N.

Local gradients per operation

Every operation needs a backward rule, a vector-Jacobian product that never builds the Jacobian. These cover most of a standard network (⊙ is element-wise multiplication):

Two patterns are worth remembering. Add copies its gradient to both inputs unchanged, which is why residual connections are such good gradient highways. And the fused softmax-plus-cross-entropy rule, "probabilities minus target", is both simpler and more numerically stable than chaining the two separate rules, which is why PyTorch's cross_entropy takes logits (see loss functions).

Time and memory

The backward pass of a linear layer does two matrix multiplies (X̄ and W̄) for the forward pass's one, so a training step costs about three forward passes: the standard estimate is 6 floating-point operations per parameter per training token, against 2 for inference alone.

Memory is the harder constraint. Every value a backward rule needs must stay alive from the moment the forward pass creates it until the backward pass consumes it, so activation memory grows with depth × batch size × width. That is why inference (inside torch.no_grad(), which records nothing) fits much larger batches than training, and why large models use activation checkpointing: save only some layers' inputs and recompute the rest during the backward pass, trading about one extra forward pass for a large cut in memory.

Vanishing and exploding gradients

Going down through layer ℓ multiplies the gradient by that layer's transposed Jacobian: h̄_ℓ−1 = W_ℓᵀ (h̄_ℓ ⊙ φ′(z_ℓ)). After L layers the gradient at the input is a product of L such factors. If each factor shrinks vectors by a typical ratio r, the gradient scales like r^L: with r = 0.7 and 30 layers that is 2 × 10⁻⁵; with r = 1.4 it is 2 × 10⁴. Early layers then learn nothing (vanishing), or updates are so large that training diverges and the loss turns into nan (exploding).

Two things set r. The weight scale: a matrix with entries of variance σ²/n scales a vector's length by about σ. And the activation's slope: ReLU passes about half the units with slope 1, so it halves the squared length, which the factor √2 in He initialization exactly compensates. The sigmoid's slope is at most ¼, so a sigmoid network shrinks the gradient by at least 4 per layer whatever the weights; this is why deep networks moved from sigmoid and tanh to ReLU.

Same weights, four scalings. Only He initialization with ReLU keeps the gradient near 1 all the way down (1.4e+00 at the input). Scale the weights by 0.7 or 1.4 and the gradient shrinks to 3.1e-05 or grows to 3.3e+04; the sigmoid network is at 3.3e-20.

The fixes are now standard, and most of them get their own page: careful initialization (Xavier/Glorot for tanh, He/Kaiming for ReLU) so that r ≈ 1 at the start; ReLU-family activations, which don't saturate for positive inputs; normalization layers, which reset each layer's scale; residual connections, which make each layer compute x + f(x) so its Jacobian is I + ∂f/∂x and the gradient always has a path that isn't multiplied down; and gradient clipping (clip_grad_norm_), which rescales the whole gradient when its norm exceeds a limit, as a guard against the occasional explosive step. Recurrent networks, which apply the same matrix at every time step, suffer the worst version of this, which is what the LSTM was invented to fix.

Gradient checking

A wrong backward rule doesn't crash: the network still trains, just badly, and the bug can hide for weeks. The test is to compare each analytic gradient with a finite difference. The forward difference (L(θ + h) − L(θ))/h has error proportional to h; the central difference

∂L/∂θ ≈ (L(θ + h) − L(θ − h)) / 2h

has error proportional to h², because the second-order terms cancel. But h can't go to zero: each loss is only accurate to about 10⁻¹⁶ of its size in float64, so the rounding error of the difference grows like 10⁻¹⁶/h. The best h balances the two, and the result is a V shape. Compare with the relative error |a − n| / max(|a|, |n|), not the absolute one, so that large and small gradients are judged alike; in float64, a relative error below about 10⁻⁷ means the backward pass is right, and one near 10⁻² means a bug.

Too large an h measures curvature, too small an h measures rounding. Central differences reach about 10⁻¹¹ here at h around 10⁻⁵; in float32, rounding is 10⁹ times worse, which is why gradient checks run in float64.

Two cautions. Kinks: if a ReLU input sits within h of zero, the two evaluations straddle the kink and the check fails even for a correct backward pass, so check at points away from kinks or ignore isolated failures. And cost: one check needs two forward passes per parameter, so check a tiny network, or a random sample of parameters, never a full model. PyTorch's torch.autograd.gradcheck does exactly this, perturbing each input entry by eps=1e-6 with central differences, requiring float64 inputs, and comparing the numerical and analytic Jacobians with atol=1e-5 and rtol=1e-3.

In PyTorch

PyTorch builds the graph while the forward pass runs: every operation on a tensor that requires grad records a node holding its backward rule and saved tensors (the chain rule page shows an autograd.Function written by hand). loss.backward() calls torch.autograd.backward, which hands the graph to the C++ engine. The engine runs the nodes in dependency order and, at each leaf tensor (a parameter), adds the result into .grad. Adding rather than overwriting is what makes shared weights and gradient accumulation work, and it's why every training loop calls optimizer.zero_grad() (or model.zero_grad()) before the next backward pass. After backward, the graph's saved tensors are freed; calling backward through it a second time needs retain_graph=True.

Worked example

Implementation

import torch
from torch import nn
import torch.nn.functional as F

model = nn.Sequential(nn.Linear(2, 3), nn.ReLU(), nn.Linear(3, 2))
with torch.no_grad():
    model[0].weight.copy_(torch.tensor([[1., 0.], [-1., 1.], [0.5, -1.]]))
    model[0].bias.copy_(torch.tensor([0., 0., 0.5]))
    model[2].weight.copy_(torch.tensor([[1., -1., 2.], [0., 1., 1.]]))
    model[2].bias.zero_()

x, t = torch.tensor([[1., 2.]]), torch.tensor([0])
loss = F.cross_entropy(model(x), t)       # 1.3133
loss.backward()                           # fills .grad on all four parameters
model[0].weight.grad                      # [[-0.731, -1.462], [1.462, 2.924], [0, 0]]

with torch.no_grad():                     # one SGD step, lr = 0.1
    for p in model.parameters():
        p -= 0.1 * p.grad
model.zero_grad()                         # .grad accumulates; clear it
F.cross_entropy(model(x), t)              # 0.1894

# gradient checking a small function in float64
def f(w):
    return F.cross_entropy(torch.relu(x.double() @ w.T) @ model[2].weight.double().T, t)
w = model[0].weight.detach().double().requires_grad_()
torch.autograd.gradcheck(f, (w,))         # True

The backward pass by hand, for a batch, in plain Python:

import math

def matmul(a, b):
    return [[sum(x * y for x, y in zip(row, col)) for col in zip(*b)] for row in a]

def transpose(m):
    return [list(col) for col in zip(*m)]

def forward_backward(X, T, W1, b1, W2, b2):
    N = len(X)
    Z = [[z + b for z, b in zip(row, b1)] for row in matmul(X, transpose(W1))]
    H = [[max(0.0, z) for z in row] for row in Z]
    S = [[s + b for s, b in zip(row, b2)] for row in matmul(H, transpose(W2))]
    P = [[math.exp(s) / sum(math.exp(v) for v in row) for s in row] for row in S]
    loss = -sum(math.log(P[i][T[i]]) for i in range(N)) / N
    dS = [[(p - (k == T[i])) / N for k, p in enumerate(P[i])] for i in range(N)]
    dW2, db2 = matmul(transpose(dS), H), [sum(c) for c in zip(*dS)]
    dH = matmul(dS, W2)
    dZ = [[g if z > 0 else 0.0 for g, z in zip(gr, zr)] for gr, zr in zip(dH, Z)]
    dW1, db1 = matmul(transpose(dZ), X), [sum(c) for c in zip(*dZ)]
    return loss, (dW1, db1, dW2, db2)

forward_backward([[1, 2]], [0], [[1, 0], [-1, 1], [0.5, -1]], [0, 0, 0.5],
                 [[1, -1, 2], [0, 1, 1]], [0, 0])[0]       # 1.3133

Tradeoffs

• Reverse mode versus finite differences. Backprop gives every gradient for about three forward passes, exactly up to rounding; finite differences cost two forward passes per parameter and are only approximate. Finite differences are a test, never a training method.
• Compute versus memory. Saving every activation is fastest; activation checkpointing recomputes them to save memory. Large-model training almost always checkpoints, paying roughly a third more compute.
• Fused versus composed rules. Writing softmax and cross-entropy as one operation gives a simpler, more stable gradient than chaining their rules; the same is true of log-sum-exp, attention and normalization layers, which is why frameworks ship fused kernels.
• Depth costs gradient health. Each layer multiplies the gradient by another Jacobian. Depth is only usable with the right initialization, activations, normalization and residual connections, and gradient clipping as a last guard.
• Biological plausibility. Backprop needs the exact transposed weights on the way back and a global error signal, which real neurons don't obviously have. That matters for neuroscience, not for engineering: no alternative training rule matches it at scale.

Check yourself

Connections

References

• Linnainmaa. The Representation of the Cumulative Rounding Error of an Algorithm as a Taylor Expansion of the Local Rounding Errors. Master's thesis, University of Helsinki, 1970 (reverse-mode differentiation).
• Werbos. Beyond Regression: New Tools for Prediction and Analysis in the Behavioral Sciences. PhD thesis, Harvard, 1974.
• Rumelhart, Hinton and Williams. Learning Representations by Back-Propagating Errors. Nature, 1986.
• Hochreiter. Untersuchungen zu dynamischen neuronalen Netzen. Diploma thesis, TU Munich, 1991; Bengio, Simard and Frasconi. Learning Long-Term Dependencies with Gradient Descent is Difficult. IEEE Transactions on Neural Networks, 1994 (vanishing gradients).
• Glorot and Bengio. Understanding the Difficulty of Training Deep Feedforward Neural Networks. AISTATS, 2010.
• He, Zhang, Ren and Sun. Delving Deep into Rectifiers. ICCV, 2015 (arXiv 1502.01852).
• Pascanu, Mikolov and Bengio. On the Difficulty of Training Recurrent Neural Networks. ICML, 2013 (arXiv 1211.5063; gradient clipping).
• Chen, Xu, Zhang and Guestrin. Training Deep Nets with Sublinear Memory Cost. 2016 (arXiv 1604.06174; activation checkpointing).
• Baydin, Pearlmutter, Radul and Siskind. Automatic Differentiation in Machine Learning: a Survey. JMLR, 2018 (arXiv 1502.05767).
• PyTorch v2.15.0: torch/autograd/__init__.py, torch/autograd/function.py and torch/autograd/gradcheck.py.