Backpropagation

Backpropagation, short for “backward propagation of errors”, computes the gradient of a neural network’s loss with respect to every weight and bias. It runs the network forward on a batch of images, keeps the intermediate results, and applies the chain rule layer by layer from the loss back to the input. Gradient descent and its stochastic variants use this gradient to update the weights. In numerical-analysis terms, backpropagation is the reverse mode of automatic differentiation applied to the network’s computation graph.

Problem

A network maps an image $\mathbf{x}$ to class scores through parameters $\theta$, and a loss $L(\theta)$ measures how wrong the scores are on a batch of labeled images. Gradient descent needs $\nabla_\theta L$, one partial derivative per parameter, often millions of them. Finite differences would need an extra forward pass or two per parameter; backpropagation needs one forward and one backward pass in total.

Inputs and Outputs

• Input: a differentiable network (a composition of layers), its current parameters, a batch of labeled images, and a scalar loss.
• Output: the loss and its gradient with respect to every parameter, in arrays shaped like the parameters. The gradient with respect to the input pixels costs one more backward step.

Intuition

Write the computation as a computation graph: nodes are intermediate arrays (pixels, weighted sums, activations, scores, the loss), and edges record which arrays each node was computed from. The backward pass asks, for each node, how much the loss would change if that node changed slightly. That answer is assembled from the answers already known for the nodes it feeds, so starting at the loss, where the answer is 1, and working in reverse order gives every answer once. Each intermediate derivative is computed once and reused by everything upstream of it, which makes the method cheap.

Algorithm

1. Forward pass. Evaluate the network layer by layer on the batch and store every layer’s pre-activation $\mathbf{z}^{(\ell)}$ and activation $\mathbf{a}^{(\ell)}$, then the loss.
2. Initialize the error signal. Compute $\boldsymbol{\delta}^{(L)} = \partial L / \partial \mathbf{z}^{(L)}$ at the output layer.
3. Backward pass. For $\ell = L, L-1, \dots, 1$: compute the layer’s parameter gradients from $\boldsymbol{\delta}^{(\ell)}$ and the stored input $\mathbf{a}^{(\ell-1)}$, then send the error signal to the previous layer to obtain $\boldsymbol{\delta}^{(\ell-1)}$.
4. Return the gradients; an optimizer then updates the parameters.

Mathematical Formulation

Chain rule. If $L$ depends on $\mathbf{u}$ only through $\mathbf{v} = g(\mathbf{u})$, then $\partial L/\partial \mathbf{u} = J_g^\top\, \partial L/\partial \mathbf{v}$, where $J_g$ is the Jacobian of $g$. For a chain of layers the gradient is a product of Jacobians. Backpropagation evaluates that product from the loss end, so every step multiplies a vector by a transposed Jacobian and no Jacobian matrix is ever formed. Evaluating from the input end (forward mode) needs one pass per input variable, the reverse order one pass per output; for a scalar loss the reverse order wins (Baydin et al., 2018). Goodfellow, Bengio and Courville (2016, §6.5.2) describe backpropagation as the chain rule computed in a highly efficient order.

Fully connected layers. Writing $\mathbf{a}^{(\ell)}$ for the neural network article’s $\mathbf{h}^{(\ell)}$, layer $\ell$ computes $\mathbf{z}^{(\ell)} = \mathbf{W}^{(\ell)}\mathbf{a}^{(\ell-1)} + \mathbf{b}^{(\ell)}$ and $\mathbf{a}^{(\ell)} = \phi(\mathbf{z}^{(\ell)})$. Define the error signal $\boldsymbol{\delta}^{(\ell)} = \partial L/\partial \mathbf{z}^{(\ell)}$. Then, for one image,

$$\frac{\partial L}{\partial \mathbf{W}^{(\ell)}} = \boldsymbol{\delta}^{(\ell)} \big(\mathbf{a}^{(\ell-1)}\big)^\top, \qquad \frac{\partial L}{\partial \mathbf{b}^{(\ell)}} = \boldsymbol{\delta}^{(\ell)}, \qquad \boldsymbol{\delta}^{(\ell-1)} = \Big(\big(\mathbf{W}^{(\ell)}\big)^\top \boldsymbol{\delta}^{(\ell)}\Big) \odot \phi'\big(\mathbf{z}^{(\ell-1)}\big),$$

where $\odot$ is the elementwise product and $\phi'$ the derivative of the activation function. The weight gradient is an outer product: weight $W_{ij}$ connects input $j$ to unit $i$, so its gradient is the error at the unit times the activity at the input. For a batch, with images stored as rows, the summed outer products form the matrix product $\boldsymbol{\Delta}^\top \mathbf{A}$. The error travels backward through the transpose of the weight matrix that carried the activations forward. For a softmax output with cross-entropy loss and one-hot label $\mathbf{y}$, the output error is $\boldsymbol{\delta}^{(L)} = \mathbf{p} - \mathbf{y}$, predicted probabilities minus label. Rumelhart, Hinton and Williams (1986) derived the same recursion for sigmoid units and a squared-error loss.

Convolution layers. A convolutional layer computes a cross-correlation (see convolution). For one input channel, kernel $w$ and “valid” output,

$$z[i,j] = \sum_{u,v} w[u,v]\, x[i+u,\, j+v] + b .$$

Differentiating with error signal $\delta[i,j] = \partial L/\partial z[i,j]$:

$$\frac{\partial L}{\partial w[u,v]} = \sum_{i,j} \delta[i,j]\, x[i+u,\, j+v], \qquad \frac{\partial L}{\partial b} = \sum_{i,j} \delta[i,j], \qquad \frac{\partial L}{\partial x[m,n]} = \sum_{i,j} \delta[i,j]\, w[m-i,\, n-j].$$

In the convolution article’s notation, the kernel gradient is the cross-correlation $x \star \delta$ of the input image with the error map. The input gradient is the full convolution $\delta * w$, the same as cross-correlating the zero-padded error map with the kernel rotated by 180 degrees. Weight sharing appears as the sum over positions $(i,j)$: a kernel weight collects contributions from every place it was used. With several channels, the input gradient also sums over output channels.

Complexity

Time. Each forward operation spawns one backward operation, so the gradient of a scalar function costs at most a constant multiple $\omega$ of evaluating it, whatever the number of parameters. Griewank (2012) calls this the cheap gradient principle; $\omega$ depends on the cost measure and the computing platform, and is exactly 3 when counting multiplications for polynomial operations. Baydin et al. (2018), citing Griewank and Walther’s textbook, give a guaranteed bound below 6, typically 2 to 3. A fully connected layer’s forward pass is one matrix product and its backward pass two (for $\partial L/\partial \mathbf{W}$ and $\boldsymbol{\delta}^{(\ell-1)}$); Chen et al. (2016) likewise call the backward operation nearly twice as costly as the forward one.

Memory. The backward pass needs the intermediate values in reverse order, so the basic reverse mode stores them all, and its memory grows with the length of the forward computation (Griewank, 2012). For a vision network that means every layer’s feature maps for the batch. Chen, Xu, Zhang and Guestrin (2016) adapted checkpointing from the automatic differentiation literature: keep only some feature maps and recompute the rest during the backward pass. Splitting an $n$-layer chain into $\sqrt{n}$ segments reduces feature-map memory to $O(\sqrt{n})$ for the price of one extra forward pass. In their experiments this cost about 30% more running time and cut the memory needed to train a 1,000-layer residual network on 224×224 ImageNet images from 48 GB to 7 GB.

Implementation

The example implements a small convolutional network in NumPy and trains it on the 8×8 handwritten-digit images bundled with scikit-learn: eight 3×3 kernels (each output 6×6, followed by ReLU), a fully connected layer from the 288 features to 10 scores, and a softmax cross-entropy loss. backward follows the formulas above, including the input-image gradient with the rotated kernel. Before training, every gradient (all 2,970 parameters and all 1,280 pixels of 20 images) is compared with centered finite differences; the printed relative error is $\lVert \mathbf{g} - \mathbf{g}_{\text{num}} \rVert / (\lVert \mathbf{g} \rVert + \lVert \mathbf{g}_{\text{num}} \rVert)$ per array.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
from sklearn.datasets import load_digits

digits = load_digits()
X = digits.images / 16.0                    # (1797, 8, 8), pixel values in [0, 1]
y = digits.target
rng = np.random.default_rng(2)
perm = rng.permutation(len(y))
train, test = perm[:1347], perm[1347:]      # 1,347 training and 450 test images

K = 8                                       # number of 3x3 filters
W1 = rng.normal(0, 1 / 3, (K, 3, 3))        # conv kernels (fan-in 9)
b1 = np.full(K, 0.1)                        # small positive bias, see the gradient check
W2 = rng.normal(0, 1 / np.sqrt(K * 36), (10, K * 36))   # fully connected 288 -> 10
b2 = np.zeros(10)
params = {"W1": W1, "b1": b1, "W2": W2, "b2": b2}

def forward(x, p):
    """x: (N, 8, 8) images. Returns softmax probabilities and what backward needs."""
    patches = sliding_window_view(x, (3, 3), axis=(1, 2))         # (N, 6, 6, 3, 3)
    z1 = np.einsum("nijuv,kuv->nkij", patches, p["W1"]) + p["b1"][:, None, None]
    a1 = np.maximum(0, z1)                                         # (N, K, 6, 6)
    h = a1.reshape(len(x), -1)                                     # flatten to (N, 288)
    s = h @ p["W2"].T + p["b2"]                                    # class scores (N, 10)
    s = s - s.max(axis=1, keepdims=True)
    P = np.exp(s) / np.exp(s).sum(axis=1, keepdims=True)           # softmax
    return P, (x, patches, z1, h)

def loss_fn(x, t, p):
    P, _ = forward(x, p)
    return -np.mean(np.log(P[np.arange(len(t)), t]))               # cross-entropy

def backward(P, t, cache, p):
    x, patches, z1, h = cache
    N = len(t)
    d_s = P.copy()
    d_s[np.arange(N), t] -= 1
    d_s /= N                                   # dL/ds: (softmax - one-hot) / N
    g = {"W2": d_s.T @ h,                      # dW = delta a^T, summed over the batch
         "b2": d_s.sum(axis=0)}
    d_h = d_s @ p["W2"]                        # error sent back to the hidden layer
    d_z1 = d_h.reshape(z1.shape) * (z1 > 0)    # times ReLU'(z)
    g["W1"] = np.einsum("nkij,nijuv->kuv", d_z1, patches)    # correlate input with delta
    g["b1"] = d_z1.sum(axis=(0, 2, 3))
    # Gradient w.r.t. the input image: full convolution of delta with each kernel,
    # i.e. correlation of the zero-padded delta with the kernel rotated 180 degrees.
    d_pad = np.pad(d_z1, ((0, 0), (0, 0), (2, 2), (2, 2)))
    win = sliding_window_view(d_pad, (3, 3), axis=(2, 3))   # (N, K, 8, 8, 3, 3)
    g["x"] = np.einsum("nkijuv,kuv->nij", win, p["W1"][:, ::-1, ::-1])
    return g

# ---- gradient check on 20 training images, every parameter and every pixel ----
xb, tb = X[train[:20]], y[train[:20]]
P, cache = forward(xb, params)
g = backward(P, tb, cache, params)
eps = 1e-6
for name in ["W1", "b1", "W2", "b2", "x"]:
    arr = xb if name == "x" else params[name]
    num = np.zeros_like(arr)
    for idx in np.ndindex(arr.shape):
        old = arr[idx]
        arr[idx] = old + eps; lp = loss_fn(xb, tb, params)
        arr[idx] = old - eps; lm = loss_fn(xb, tb, params)
        arr[idx] = old
        num[idx] = (lp - lm) / (2 * eps)                  # centered difference
    rel = np.linalg.norm(g[name] - num) / (np.linalg.norm(g[name]) + np.linalg.norm(num))
    print(f"{name:2s} {str(arr.shape):12s} {arr.size:5d} entries  relative error {rel:.1e}")

# ---- train with mini-batch gradient descent ----
def accuracy(idx):
    P, _ = forward(X[idx], params)
    return (P.argmax(axis=1) == y[idx]).mean()

lr, batch = 0.02, 16
print(f"epoch 0: loss {loss_fn(X[train], y[train], params):.3f}, test accuracy {accuracy(test):.3f}")
for epoch in range(1, 9):
    order = rng.permutation(train)
    for i in range(0, len(order), batch):
        idx = order[i:i + batch]
        P, cache = forward(X[idx], params)
        g = backward(P, y[idx], cache, params)
        for name in params:
            params[name] -= lr * g[name]
    print(f"epoch {epoch}: loss {loss_fn(X[train], y[train], params):.3f}, "
          f"test accuracy {accuracy(test):.3f}")

Output:

W1 (8, 3, 3)       72 entries  relative error 3.9e-09
b1 (8,)             8 entries  relative error 2.2e-09
W2 (10, 288)     2880 entries  relative error 3.2e-09
b2 (10,)           10 entries  relative error 2.0e-09
x  (20, 8, 8)    1280 entries  relative error 1.6e-08
epoch 0: loss 2.317, test accuracy 0.111
epoch 1: loss 1.379, test accuracy 0.747
epoch 2: loss 0.722, test accuracy 0.860
epoch 3: loss 0.492, test accuracy 0.904
epoch 4: loss 0.414, test accuracy 0.911
epoch 5: loss 0.348, test accuracy 0.878
epoch 6: loss 0.303, test accuracy 0.902
epoch 7: loss 0.249, test accuracy 0.944
epoch 8: loss 0.264, test accuracy 0.927

Every analytic gradient, including the input-image gradient, agrees with finite differences to about eight digits. Training lowers the loss from 2.32 to 0.26 and raises test accuracy from chance to 92.7% in eight epochs; over ten seeds the median was 92.8% (range 87% to 96%). With a step size of 0.1 the training loss rose between epochs in every seed (largest rise 0.13 to 1.48, against at most 0.10 at 0.02). The gradients are equally exact, so the rises come from the update, not from backpropagation; our reading is that some steps overshoot where the mini-batch loss curves sharply, the step-size limit described in the gradient descent article.

The convolution biases start at 0.1 for a reason. With zero biases, every all-black 3×3 patch gives a pre-activation of exactly zero, on the ReLU’s kink, and the centered difference straddles the kink. In that variant the check reported relative errors of about $4 \times 10^{-2}$ for b1 and the pixels, although the backward pass was correct.

PyTorch’s autograd records the computation graph during the forward pass and runs reverse-mode differentiation on loss.backward(). Continuing from the code above (run with PyTorch 2.14.1), with the trained weights and 64 test images:

import torch
import torch.nn.functional as F

xb, tb = X[test[:64]], y[test[:64]]
P, cache = forward(xb, params)
g = backward(P, tb, cache, params)              # NumPy gradients

T = {k: torch.tensor(v, requires_grad=True) for k, v in params.items()}
x_t = torch.tensor(xb, requires_grad=True)
z1 = F.conv2d(x_t[:, None], T["W1"][:, None], T["b1"])   # a cross-correlation, like ours
s = F.relu(z1).flatten(1) @ T["W2"].T + T["b2"]
loss = F.cross_entropy(s, torch.tensor(tb))
loss.backward()                                  # reverse-mode autodiff

print(f"torch {torch.__version__}, loss {loss.item():.6f}")
for k in ["W1", "b1", "W2", "b2"]:
    print(f"{k}: max |autograd - NumPy| = {(T[k].grad.numpy() - g[k]).__abs__().max():.1e}")
print(f"x : max |autograd - NumPy| = {abs(x_t.grad.numpy() - g['x']).max():.1e}")

Output:

torch 2.14.1, loss 0.205382
W1: max |autograd - NumPy| = 9.0e-17
b1: max |autograd - NumPy| = 7.8e-18
W2: max |autograd - NumPy| = 5.6e-17
b2: max |autograd - NumPy| = 1.4e-17
x : max |autograd - NumPy| = 2.1e-17

The automatically generated backward pass agrees with ours to float64 rounding error.

Properties and Behavior

• Exact, not approximate. Backpropagation computes the gradient of the function the network evaluates, up to floating-point rounding, with no step size or truncation. Where an activation is not differentiable, such as the ReLU at $z = 0$, an implementation picks a fixed value (0 in the code above).
• Gradient checking. A hand-written backward pass is verified against finite differences, preferably the centered difference $(L(\theta + h) - L(\theta - h))/2h$, whose error shrinks as $h^2$. Goodfellow et al. (2016, §11.5) recommend this test when implementing a new operation, with random projections when a full Jacobian is too expensive to check. The step $h$ must be small but not so small that rounding dominates; float64 and $h = 10^{-6}$ worked here.
• Vanishing and exploding gradients. The recursion multiplies the error signal by $(\mathbf{W}^{(\ell)})^\top$ and by $\phi'$ once per layer. If these factors tend to shrink vectors, the error reaching the early layers decays exponentially with depth; if they stretch them, it grows exponentially. Schmidhuber (2015) credits Sepp Hochreiter’s 1991 diploma thesis with formally identifying the problem, and Bengio, Simard and Frasconi (1994) analyzed why gradient-based learning in recurrent networks becomes harder as the time span of the dependencies grows. Schmidhuber lists LSTM-like recurrent networks among the ways around the problem; in feedforward networks, the choice of activation, normalization layers and residual connections affect how well gradients pass through many layers.

Limitations

• Differentiability. Discrete steps such as arg-max, rounding, or sampling pass no useful gradient and need surrogates.
• Biological plausibility. Rumelhart et al. (1986) themselves noted that the procedure, in its current form, is not a plausible model of learning in brains.

Variants

• Reverse-mode automatic differentiation is the general form, applied to any program built from differentiable operations, including loops and branches; frameworks such as PyTorch provide it for models written as ordinary programs (Baydin et al., 2018).
• Backpropagation through time unrolls a recurrent network over the steps of a sequence and backpropagates through the unrolled graph, with weights shared across steps; Rumelhart et al. (1986) already described this mapping.

History

Griewank (2012) recounts that Seppo Linnainmaa described the reverse mode in his 1970 master’s thesis at the University of Helsinki, as a tool for estimating how rounding errors propagate through a computation, and published results from it in English in 1976. Efficient chain-rule computations based on dynamic programming appeared in the 1960s and 1970s, mostly for control applications (Goodfellow et al., 2016, §6.6). For neural networks, Schmidhuber (2015) points to a preliminary discussion in Paul Werbos’s 1974 Harvard PhD thesis and to Werbos’s 1981 paper as the first network-specific application; Goodfellow et al. also credit Werbos (1981), and list LeCun (1985), Parker (1985) and Rumelhart et al. (1986) as independent rediscoveries.

Rumelhart, Hinton and Williams (1986) showed in a short Nature paper that backpropagation lets hidden units learn useful internal representations, for example on symmetry detection and family-tree relationships, and acknowledged Parker’s and LeCun’s independent work. Schmidhuber credits this paper, and Goodfellow et al. the companion Parallel Distributed Processing chapter, with popularizing the method. LeCun et al. (1989) brought it to images: their network read 16×16 normalized images of handwritten zip-code digits from U.S. mail directly, rather than feature vectors, used feature maps whose units shared 5×5 kernels, and was trained end to end. The authors reported a 5.0% error rate on 2,007 test digits after 23 passes through the 7,291 training digits.