47 - Backpropagation derived in full, step by step

An assembly line produces a defective object. Who is responsible? It is not enough to say "the last station": perhaps it did its job well on a piece already ruined by the previous one. To fix the line you have to work back, station after station, computing for each one how much it contributed to the final defect. If you know that station 3 contributed to a certain degree, and that station 3 depends on station 2, you can obtain the contribution of station 2 by combining the contribution of station 3 with how much station 3 depends on station 2. It is an accounting of responsibilities that proceeds backward. Backpropagation (backward propagation of the error) is exactly this accounting, and its mathematical language is only one: the chain rule of differential calculus.

This is the pivotal lesson of the course. A neural network is trained with gradient descent, and gradient descent needs the gradient of the cost function with respect to every weight of the network, which can number in the millions. Backpropagation computes them all at a cost of the same order as a single forward pass. The algorithm was made famous by D. Rumelhart, G. Hinton and R. Williams in 1986 (the backward form of computation, generally known as the reverse mode of automatic differentiation, had been described earlier, for example by S. Linnainmaa in 1970). There is nothing mysterious about it: it is the chain rule applied in an orderly way. In this lesson you derive it in full, check it with a numerical example worked by hand, implement it, and verify the result in two independent ways.

Before we start: the chain rule

You have already used the chain rule in the previous lessons. We review it precisely, because everything rests on it. If y depends on u, and u depends on x, that is y = f(u) and u = g(x), the derivative of y with respect to x is the product of the derivatives along the chain:

dy/dx = (dy/du) · (du/dx)

The intuition: if u changes 3 times faster than x, and y changes 2 times faster than u, then y changes 6 times faster than x. The "rates of change" multiply. Minimal numerical check: y = (3x)² with u = 3x, y = u². Then dy/du = 2u = 6x and du/dx = 3, so dy/dx = 18x. Differentiating y = 9x² directly gives 18x: it matches.

With more variables a novelty enters that you will use constantly. If y depends on two intermediate variables u₁ and u₂, which both depend on x, the total derivative is the sum of the contributions of all the paths:

dy/dx = (∂y/∂u₁)(du₁/dx) + (∂y/∂u₂)(du₂/dx)

Every neuron of a layer influences the loss through all the neurons of the next layer it is connected to: this is why sums will appear in the formulas that follow (which will become matrix products).

The computational graph: a single neuron

A computation can be drawn as a computational graph: every node is an elementary operation, every arrow carries a value to the next operation. Let us take the smallest case that contains all the ingredients, a neuron with a sigmoid activation and cross-entropy loss, with a single input:

z = w · x + b              (weighted sum)
a = σ(z)                   (sigmoid activation)
L = −log(a)  if y = 1      (cross-entropy for the label y = 1)

Let us choose concrete numbers: x = 1.5, w = 0.8, b = −0.5, y = 1. The forward pass gives z = 0.8·1.5 − 0.5 = 0.7; a = σ(0.7) ≈ 0.6682; L = −log(0.6682) ≈ 0.4032. This is the only calculation that needs to be done "forward"; every intermediate value must be remembered, because it will be needed on the way back.

The drawing shows in blue above each node the value of the forward pass and in red below each node the derivative of the loss with respect to that value, computed backward. Let us retrace the red calculations, from right to left, applying the chain rule one step at a time. We start from the end: the derivative of the loss with respect to itself is 1.

The first step, from the loss to the value a, uses L = −log a, so ∂L/∂a = −1/a = −1/0.6682 ≈ −1.497: if a increased by a small amount, the loss would decrease by 1.497 times that amount. The second step crosses the sigmoid: ∂a/∂z = a(1 − a) = 0.6682·0.3318 ≈ 0.2217 (the derivative of the sigmoid, proved in the lesson on logistic regression), and by the chain rule ∂L/∂z = (∂L/∂a)·(∂a/∂z) = −1.497 · 0.2217 ≈ −0.3318. This coincides with a − y = 0.6682 − 1 = −0.3318, the "output minus target" simplification that you had derived for cross-entropy, and it is the reason why it is convenient to stop at z instead of going further.

From z the gradient is distributed to the ingredients of the weighted sum, and here the order is no longer relevant because the three calculations are independent. For the bias, z = wx + b gives ∂z/∂b = 1, so ∂L/∂b = (∂L/∂z)·1 = −0.3318: the sum "lets through" the gradient unchanged. For the weight, ∂z/∂w = x = 1.5, so ∂L/∂w = −0.3318·1.5 ≈ −0.4977: the product sends the gradient to the weight scaled by the other factor, the input. For the input, ∂z/∂x = w = 0.8, so ∂L/∂x ≈ −0.2654; in the single neuron the input is data and is not needed, but in a network x is the output of the previous layer, and this is the gradient that propagates backward.

Two practical rules emerge, and they are all that is needed for any graph: a sum distributes the gradient identically to all the addends; a product distributes it by multiplying it by the other factor. Let us verify the numbers with code, and check them independently with the finite differences seen in the lesson on gradient descent:

import numpy as np

sigmoid = lambda z: 1 / (1 + np.exp(-z))
x, w, b, y = 1.5, 0.8, -0.5, 1.0

def loss(w, b, x):                                 # cross-entropy with label y = 1
    return -np.log(sigmoid(w * x + b))

# forward pass, remembering the intermediate values
z = w * x + b
a = sigmoid(z)
L = loss(w, b, x)

# backward pass, node by node
dL_da = -1 / a
da_dz = a * (1 - a)
dL_dz = dL_da * da_dz                              # chain rule
dL_dw, dL_db, dL_dx = dL_dz * x, dL_dz * 1.0, dL_dz * w

print("forward:   z = %.4f   a = %.4f   L = %.4f" % (z, a, L))
print("backward:  dL/da = %.4f   da/dz = %.4f   dL/dz = %.4f  (a - y = %.4f)" % (dL_da, da_dz, dL_dz, a - y))
print("gradients: dL/dw = %.5f   dL/db = %.5f   dL/dx = %.5f" % (dL_dw, dL_db, dL_dx))

eps = 1e-6
num = ((loss(w + eps, b, x) - loss(w - eps, b, x)) / (2 * eps),
       (loss(w, b + eps, x) - loss(w, b - eps, x)) / (2 * eps),
       (loss(w, b, x + eps) - loss(w, b, x - eps)) / (2 * eps))
print("numeric:   dL/dw = %.5f   dL/db = %.5f   dL/dx = %.5f" % num)

forward:   z = 0.7000   a = 0.6682   L = 0.4032
backward:  dL/da = -1.4966   da/dz = 0.2217   dL/dz = -0.3318  (a - y = -0.3318)
gradients: dL/dw = -0.49772   dL/db = -0.33181   dL/dx = -0.26545
numeric:   dL/dw = -0.49772   dL/db = -0.33181   dL/dx = -0.26545

What the code does. First it computes the forward pass, saving z and a (the two intermediate quantities needed on the way back). Then it walks the graph backward: each line is a multiplication of a derivative "arriving" from the right by a "local" derivative of the node. The three final gradients coincide with those obtained by perturbing the parameters by a tiny amount (central finite differences), down to the last printed digit (five decimals). It is the first confirmation: the chain rule works as promised.

A layered network: the general derivation

Now the same idea on a network with L layers and a batch of m examples. We repeat the notation of the previous lesson, with the first dimension always the batch:

A⁰ = X                                  (m, n₀)
Z^l = A^{l−1} · W^l + b^l               (m, n_l)       pre-activation of layer l
A^l = f_l(Z^l)                          (m, n_l)       activation (applied element by element)
J   = (1/m) Σ_i ℓ(A^L_i, y_i)                          mean loss over the batch

We want ∂J/∂W^l and ∂J/∂b^l for every layer l. The key is to introduce an intermediate quantity, the error term (or "delta") of layer l: the derivative of the loss with respect to the pre-activation.

δ^l = ∂J/∂Z^l          shape (m, n_l)   (one row per example, one column per neuron)

Choosing the pre-activation Z (and not the activation A) as the reference variable is convenient because Z is the point where the weights "enter": once we know δ^l, the gradients of the weights of layer l come out immediately. We proceed in three steps, examining a single example and a single neuron, and then rewriting in matrix form.

Step 1: the gradients of the weights, given δ

The pre-activation of neuron k of layer l is Z^l_k = Σ_j A^{l−1}_j W^l_{jk} + b^l_k. It depends on the weight W^l_{jk} only through the term A^{l−1}_j · W^l_{jk}, so ∂Z^l_k/∂W^l_{jk} = A^{l−1}_j and ∂Z^l_k/∂b^l_k = 1. By the chain rule (and since W^l_{jk} influences the loss only through Z^l_k):

∂J/∂W^l_{jk} = δ^l_k · A^{l−1}_j          and          ∂J/∂b^l_k = δ^l_k

In words: the gradient of a weight is the product of the error of the neuron the weight arrives at by the activation of the neuron the weight starts from. Over the whole batch the contributions of all the examples are summed (and averaged). In matrix form:

∂J/∂W^l = (1/m) · (A^{l−1})ᵀ · δ^l              shape:  (n_{l−1}, m) · (m, n_l) = (n_{l−1}, n_l)   ✓ like W^l
∂J/∂b^l = (1/m) · Σ_rows δ^l                     shape:  (n_l,)                                      ✓ like b^l

The transpose and the sum over the rows come from the sum over all the examples of the batch. Checking the shapes is already a verification: the gradient of a parameter must always have the same shape as the parameter.

Step 2: the error of the final layer

In the last layer, δ^L = ∂J/∂Z^L is computed directly. For the three standard activation-loss combinations you have already derived the result in the previous lessons, and it is always the same:

linear output + squared error (½(ŷ − y)²):            δ^L = ŷ − y
sigmoid output + binary cross-entropy:                 δ^L = a − y
softmax output + categorical cross-entropy:            δ^L = P − Y     (Y = labels in one-hot encoding)

(The factor 1/m of the batch mean is applied only once, in the formula for the gradient of the weights.) This is the starting point of the path backward: output minus target.

Step 3: the backward recurrence

We have δ^L. How do we obtain δ^l starting from δ^{l+1}? Neuron j of layer l influences the loss through all the neurons k of the next layer (sum over the paths), because Z^{l+1}_k = Σ_j A^l_j W^{l+1}_{jk} + b^{l+1}_k depends on A^l_j with coefficient W^{l+1}_{jk}. We apply the chain rule in two stretches, first up to A^l_j and then up to Z^l_j (where A^l_j = f_l(Z^l_j)):

∂J/∂A^l_j = Σ_k (∂J/∂Z^{l+1}_k) · (∂Z^{l+1}_k/∂A^l_j) = Σ_k δ^{l+1}_k · W^{l+1}_{jk}

δ^l_j = ∂J/∂Z^l_j = (∂J/∂A^l_j) · f_l'(Z^l_j) = ( Σ_k δ^{l+1}_k · W^{l+1}_{jk} ) · f_l'(Z^l_j)

In matrix form, the sum over k is a product of δ^{l+1} by the transpose of W^{l+1}, and the product by f' is element by element (denoted by ⊙):

δ^l = ( δ^{l+1} · (W^{l+1})ᵀ ) ⊙ f_l'(Z^l)          shape: (m, n_{l+1}) · (n_{l+1}, n_l) = (m, n_l)  ✓

This is the backpropagation equation. Read it like a sentence: "the error of a layer is the error of the next layer, sent back through the weights (the transpose), and filtered by the derivative of the non-linearity". A neuron with a zero derivative (for example a switched-off ReLU) completely blocks the passage of the error: hence the phenomenon of "dead" ReLUs. And the weights act as return channels: the same weight that carried the signal forward carries the responsibility backward.

The four equations, in summary

(1) δ^L     = A^L − Y                                       (output minus target, for the three standard pairs)
(2) δ^l     = ( δ^{l+1} · (W^{l+1})ᵀ ) ⊙ f_l'(Z^l)         (for l = L−1, ..., 1)
(3) ∂J/∂W^l = (1/m) · (A^{l−1})ᵀ · δ^l
(4) ∂J/∂b^l = (1/m) · sum over the rows of δ^l

That is all. The algorithm is: (a) forward pass saving all the Z^l and A^l; (b) compute δ^L with (1); (c) for each layer, from the last to the first, compute the gradients with (3) and (4) and then, if l > 1, propagate δ backward with (2). The cost of the backward pass is of the same order as the forward one, because each layer requires matrix products of the same dimensions.

A worked example by hand: the 2 → 2 → 1 network

Let us take the network of the previous lesson: 2 inputs, 2 hidden neurons with ReLU, 1 linear output, with x = (1, 2), W¹ = [[0.5; −1.0], [0.25; 0.5]], b¹ = (0, 1), W² = [[2.0], [−1.0]], b² = 0.5. We had computed Z¹ = (1, 1), A¹ = (1, 1) and output ŷ = 1.5. We add the target y = 1 and the squared loss J = ½(ŷ − y)² = ½·0.25 = 0.125.

Let us apply the four equations, one at a time (the values are those written in red in the drawing):

The order of the steps matters, because each one uses the result of the previous one. The output error (equation 1) is δ² = ŷ − y = 1.5 − 1 = 0.5, the only quantity that starts from the loss. With equations (3-4) the gradients of layer 2 are ∂J/∂W² = (A¹)ᵀ δ² = (1, 1)ᵀ · 0.5 = (0.5; 0.5)ᵀ, of shape (2, 1) like W², and ∂J/∂b² = 0.5.

To propagate the error (equation 2) we compute δ² · (W²)ᵀ = 0.5 · (2.0; −1.0) = (1.0; −0.5): hidden neuron 1 receives 100% of the error signal multiplied by the weight +2, neuron 2 receives −50% because of the weight −1. Multiplying by the derivative of ReLU at Z¹ = (1, 1), which is 1 in both cases (positive Z), we get δ¹ = (1.0; −0.5). With (3-4) on layer 1: ∂J/∂W¹ = xᵀ δ¹ = (1; 2)ᵀ · (1.0; −0.5) = [[1.0; −0.5], [2.0; −1.0]] and ∂J/∂b¹ = (1.0; −0.5).

A readable interpretation: the gradient of the weight connecting input 2 to hidden neuron 1 is 2.0, twice that of input 1 toward the same neuron (1.0), because input 2 is twice as large. The negative sign of the gradients toward hidden neuron 2 says that to reduce the loss that neuron should increase its activation (the network is underestimating in that branch). Let us verify everything with code, once again with finite differences as an independent referee:

x = np.array([[1.0, 2.0]]); y = np.array([[1.0]])
W1 = np.array([[0.5, -1.0], [0.25, 0.5]]); b1 = np.array([0.0, 1.0])
W2 = np.array([[2.0], [-1.0]]);            b2 = np.array([0.5])

def forward_pass(W1, b1, W2, b2):
    Z1 = x @ W1 + b1; A1 = np.maximum(0, Z1); Z2 = A1 @ W2 + b2
    return Z1, A1, Z2, 0.5 * ((Z2 - y) ** 2).mean()

Z1, A1, Z2, J = forward_pass(W1, b1, W2, b2)

# backward pass, with equations (1)-(4)
d2 = Z2 - y                                        # (1) output minus target:                 (1, 1)
dW2 = A1.T @ d2                                    # (3) (A1)^T d2:                          (2, 1)
db2 = d2.sum(axis=0)                               # (4)
d1 = (d2 @ W2.T) * (Z1 > 0)                        # (2) d2 W2^T, then * ReLU derivative:    (1, 2)
dW1 = x.T @ d1                                     # (3) x^T d1:                             (2, 2)
db1 = d1.sum(axis=0)                               # (4)

print("J = %.3f   delta2 = %s   delta1 = %s" % (J, d2.ravel(), d1.ravel()))
print("dW2 =", dW2.ravel(), "  db2 =", db2, "\ndW1 =", dW1.tolist(), "  db1 =", db1)

def numeric_gradient(index, eps=1e-6):             # perturbs one parameter at a time
    params = [W1.copy(), b1.copy(), W2.copy(), b2.copy()]
    out = np.zeros_like(params[index])
    for i in np.ndindex(out.shape):
        plus = [p.copy() for p in params]; minus = [p.copy() for p in params]
        plus[index][i] += eps; minus[index][i] -= eps
        out[i] = (forward_pass(*plus)[3] - forward_pass(*minus)[3]) / (2 * eps)
    return out

for name, an, k in (("W1", dW1, 0), ("b1", db1, 1), ("W2", dW2, 2), ("b2", db2, 3)):
    print("%s: maximum difference between analytic and numeric = %.2e" % (name, np.abs(an - numeric_gradient(k)).max()))

J = 0.125   delta2 = [0.5]   delta1 = [ 1.  -0.5]
dW2 = [0.5 0.5]   db2 = [0.5] 
dW1 = [[1.0, -0.5], [2.0, -1.0]]   db1 = [ 1.  -0.5]
W1: maximum difference between analytic and numeric = 6.04e-11
b1: maximum difference between analytic and numeric = 1.98e-11
W2: maximum difference between analytic and numeric = 7.44e-12
b2: maximum difference between analytic and numeric = 7.44e-12

What the code does. The lines with the comments (1)-(4) are the four equations, transcribed almost literally. Note (Z1 > 0): it is the derivative of ReLU (true where Z is positive, false elsewhere), which NumPy converts to 1 and 0 when multiplying it. The function numeric_gradient perturbs one parameter at a time forward and backward, calling the forward pass again: it is slow (two forward passes for every parameter) but independent of the theory, and for this reason it is the ideal referee. The maximum differences are on the order of 10⁻¹¹ or less: the gradients computed by hand coincide with the numerical ones.

General implementation and gradient check on a 4 → 8 → 3 network

Let us put the four equations in a class with any number of layers. For the last layer we use softmax and cross-entropy (equation 1: δ = P − Y), for the hidden layers ReLU or sigmoid. After the forward pass, the class keeps all the Z and A: they are indispensable on the way back, and this is the reason why training a network consumes much more memory than using it for inference.

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

Xi, yi = load_iris(return_X_y=True)
Xi = StandardScaler().fit_transform(Xi)

def softmax(Z):
    Z = Z - Z.max(axis=1, keepdims=True)
    E = np.exp(Z)
    return E / E.sum(axis=1, keepdims=True)

class Network:
    def __init__(self, sizes, activation="relu", seed=0):
        rng = np.random.default_rng(seed)
        self.activation = activation
        self.W = [rng.normal(size=(i, o)) * np.sqrt(2.0 / i) for i, o in zip(sizes[:-1], sizes[1:])]
        self.b = [np.zeros(o) for o in sizes[1:]]

    def f(self, Z):
        return np.maximum(0, Z) if self.activation == "relu" else 1 / (1 + np.exp(-Z))

    def df(self, Z, A):                                    # derivative of f, in terms of Z (ReLU) or of A (sigmoid: A(1-A))
        return (Z > 0).astype(float) if self.activation == "relu" else A * (1 - A)

    def forward(self, X):
        self.Z, self.A = [], [X]                           # all pre-activations and activations are kept
        A = X
        for l, (W, b) in enumerate(zip(self.W, self.b)):
            Z = A @ W + b
            A = softmax(Z) if l == len(self.W) - 1 else self.f(Z)
            self.Z.append(Z); self.A.append(A)
        return A

    def loss(self, X, y):
        P = self.forward(X)
        return -np.log(P[np.arange(len(y)), y]).mean()

    def backward(self, y):
        m = len(y)
        delta = self.A[-1] - np.eye(self.A[-1].shape[1])[y]            # (1): output minus target (softmax + cross-entropy)
        gW, gb = [None] * len(self.W), [None] * len(self.W)
        for l in range(len(self.W) - 1, -1, -1):
            gW[l] = self.A[l].T @ delta / m                            # (3)
            gb[l] = delta.sum(axis=0) / m                              # (4)
            if l > 0:                                                  # (2): propagate backward, except beyond the first layer
                delta = (delta @ self.W[l].T) * self.df(self.Z[l - 1], self.A[l])
        return gW, gb

network = Network([4, 8, 3])
network.loss(Xi, yi)
gW, gb = network.backward(yi)

def numeric_gradient_network(network, X, y, eps=1e-5):
    result = []
    for P in network.W + network.b:
        g = np.zeros_like(P)
        for idx in np.ndindex(P.shape):
            old = P[idx]
            P[idx] = old + eps; plus = network.loss(X, y)
            P[idx] = old - eps; minus = network.loss(X, y)
            P[idx] = old; g[idx] = (plus - minus) / (2 * eps)
        result.append(g)
    return result

num = numeric_gradient_network(network, Xi, yi)
print("relative error analytic/numeric (norm of the difference / norm of the sum):")
for name, a, n in zip(["W1 (4x8)", "W2 (8x3)", "b1 (8,)", "b2 (3,)"], gW + gb, num):
    print("  %-9s %.2e" % (name, np.linalg.norm(a - n) / np.linalg.norm(a + n)))
print("parameters: %d  ->  the numeric gradient needed %d forward passes, the analytic one a forward pass plus a backward one" % (
    sum(W.size + b.size for W, b in zip(network.W, network.b)), 2 * sum(W.size + b.size for W, b in zip(network.W, network.b))))

relative error analytic/numeric (norm of the difference / norm of the sum):
  W1 (4x8)  4.80e-11
  W2 (8x3)  2.65e-11
  b1 (8,)   3.73e-11
  b2 (3,)   4.10e-11
parameters: 67  ->  the numeric gradient needed 134 forward passes, the analytic one a forward pass plus a backward one

What the code does. forward keeps all the intermediate values in self.Z and self.A. backward is a loop that goes from the last layer to the first and for each one applies (3) and (4); then, if we are not at the first layer, it propagates δ backward with (2): delta @ self.W[l].T is the product with the transpose, and self.df(self.Z[l - 1], self.A[l]) is the derivative of the activation of the previous layer. The index l - 1 for the Z (and l for the A, which have one more element, the input) is the part that is usually gotten wrong: the alignment of the indices deserves a check, and that is exactly what the gradient check does.

The relative error is on the order of 10⁻¹⁰ for all four groups of parameters: the implemented backpropagation is correct. The last line contrasts the costs: for the 67 parameters of the network the numerical gradient needed 134 forward passes; the analytic one, one forward and one backward. With a million parameters it would be two million against two: this is the reason why backpropagation made training networks practicable.

A second referee: PyTorch's automatic differentiation

The gradient check with finite differences is a referee independent of the theory. A second referee is the automatic differentiation of a library, which builds the graph while computing forward and applies the chain rule backward. The listing copies the same weights into float64 tensors with requires_grad=True (float64 because our network works in double precision; the library would use float32 by itself), redoes the forward pass with the library's cross-entropy loss (which contains the softmax) and calls backward(). It reuses network, Xi, yi, gW and gb from the previous block. It was run with PyTorch 2.12.0 on CPU (CUDA_VISIBLE_DEVICES="", 2 threads), twice, with identical output; the differences, being rounding errors, depend on the machine and the version, and have no value worth remembering.

import torch

tW = [torch.tensor(W, requires_grad=True) for W in network.W]    # same weights, in float64
tb = [torch.tensor(b, requires_grad=True) for b in network.b]
Xt, yt = torch.tensor(Xi), torch.tensor(yi)

Z1 = torch.relu(Xt @ tW[0] + tb[0])
Z2 = Z1 @ tW[1] + tb[1]
J = torch.nn.functional.cross_entropy(Z2, yt)                    # applies softmax + cross-entropy internally
J.backward()                                                     # automatic backpropagation

print("PyTorch loss %.6f   hand-written loss %.6f" % (J.item(), network.loss(Xi, yi)))
for name, mine, p in zip(["W1", "W2", "b1", "b2"], gW + gb, tW + tb):
    print("  %s  maximum difference between my gradients and PyTorch's: %.2e" % (name, np.abs(mine - p.grad.numpy()).max()))

PyTorch loss 1.572833   hand-written loss 1.572833
  W1  maximum difference between my gradients and PyTorch's: 9.71e-17
  W2  maximum difference between my gradients and PyTorch's: 1.67e-16
  b1  maximum difference between my gradients and PyTorch's: 5.55e-17
  b2  maximum difference between my gradients and PyTorch's: 1.04e-16

What the code does. It copies the same weights into PyTorch tensors, redoes the forward pass with PyTorch operations and calls J.backward(), which walks the computational graph built during the forward computation, applying the chain rule node by node. The losses coincide (1.572833) and the gradients coincide with ours up to rounding, with maximum differences on the order of 10⁻¹⁶. They are two different measures (maximum difference in absolute value here, relative error on the norms in the gradient check), but they say the same thing: the backpropagation written by hand is correct.

What backward() does is no different from what you have written: it is the same algorithm applied automatically to any graph of differentiable operations. The advantage is that you do not have to derive every new architecture by hand; the risk is not understanding what happens when something goes wrong.

Why the gradient vanishes in deep networks: you see it in the recurrence

Equation (2) shows the mechanism of the vanishing gradient (and of the exploding one). To get from the last layer to the first, δ is multiplied by a matrix (W)ᵀ and by a derivative f' at every layer: essentially a product of L−1 factors. If the typical factors are smaller than 1, the product decreases exponentially; if they are larger than 1, it grows. The derivative of the sigmoid never exceeds 0.25 (previous lesson), so a typical factor per layer is smaller than 1 and the gradients of the first layers are tiny. Let us measure it in a network with 10 layers, comparing sigmoid and ReLU (with He):

rng = np.random.default_rng(2)
Xr = rng.normal(size=(200, 20))
yr = rng.integers(0, 3, size=200)

for activation in ("sigmoid", "relu"):
    deep = Network([20] + [30] * 9 + [3], activation=activation, seed=1)
    if activation == "sigmoid":                                 # for the sigmoid I use the Xavier scale, as in the literature
        deep.W = [np.random.default_rng(i).normal(size=W.shape) * np.sqrt(1.0 / W.shape[0]) for i, W in enumerate(deep.W)]
    deep.loss(Xr, yr)
    gW_deep, _ = deep.backward(yr)
    print("%-8s norm of the weight gradient, from the first to the last layer:" % activation)
    print("         ", ["%.1e" % np.linalg.norm(g) for g in gW_deep])

sigmoid  norm of the weight gradient, from the first to the last layer:
          ['4.0e-07', '1.4e-06', '6.4e-06', '2.5e-05', '1.1e-04', '5.3e-04', '1.9e-03', '1.0e-02', '4.3e-02', '2.5e-01']
relu     norm of the weight gradient, from the first to the last layer:
          ['2.5e-01', '4.5e-01', '5.0e-01', '6.3e-01', '8.4e-01', '1.0e+00', '1.8e+00', '2.1e+00', '1.4e+00', '1.6e+00']

What the code does. It creates a network with 10 layers of weights (20 → 30 → ... → 30 → 3) and computes the gradient on 200 random examples, printing the norm of the weight gradient of every layer, from the first (near the input) to the last. With the sigmoid the norm grows by a factor between 3.5 and 6 at each layer toward the output: the first layer has a gradient of 4·10⁻⁷ against 0.25 for the last, a ratio of about six hundred thousand to one.

The initial layers practically do not learn. With ReLU and He initialization all the norms are between 0.25 and 2.1: the error signal arrives everywhere. What you have measured here on one example is, in miniature, the reason why ReLU and correct initialization made deep networks trainable.

Trade-offs and real-world traps

• Not saving the activations: the backward pass needs them. Training consumes memory proportional to depth × batch × width; this is why a model that "fits" in memory at inference may not fit in training. Techniques such as gradient checkpointing trade computation for memory.
• Getting the alignment of the indices wrong (Z^l versus A^l versus A^{l−1}): it is the typical error of a hand implementation. The gradient check always catches it; do not write a backpropagation without having run it at least once.
• Using the derivative of the sigmoid as a function of Z instead of A: the formula A(1 − A) uses the activation, not the pre-activation. Applying σ(1 − σ) to a value that is already an activation gives a wrong result without reporting it.
• Forgetting the division by m (the batch mean): the gradients come out m times larger, the effective learning rate changes with the batch size, and the gradient check fails.
• Reducing the bias dimensions wrongly: the gradient of b is obtained by summing δ over the rows (axis 0); summing over the wrong axis gives wrong shapes or, worse, "right" shapes by chance when n_l = m.
• Gradient check with ReLU: near Z = 0 the derivative does not exist, and finite differences may deviate. If a check fails by a small amount on a single element, try another seed before looking for a bug.
• The gradient is not the update: backpropagation computes ∂J/∂θ. It is the optimizer (SGD, Adam) that decides how to use it to update the weights.
• Exploding gradients: the opposite case of the vanishing gradient. They are limited with gradient clipping (rescaling the gradient if its norm exceeds a threshold) and with correct initialization.
• Libraries accumulate gradients: in PyTorch, without optimizer.zero_grad() at every step, backward() adds them to the previous ones (this is intended, to sum over several batches, but it is a source of bugs when forgotten).

The nine entries fall into three families. Some are implementation errors that the gradient check always intercepts (indices, derivative of the sigmoid, division by m, bias axis). Others are costs and limits of the algorithm (memory for the activations, vanishing or exploding gradients). The last are conceptual confusions: the gradient is not the update, and the gradient check with ReLU can fail by a small amount near zero without there being a bug.

Frequently asked questions

Is backpropagation a learning algorithm? No, it is an algorithm for computing gradients. Learning is gradient descent (or a variant of it) that uses those gradients. They are often confused because they are presented together.

Why do we go backward and not forward? There would also be a "forward" way of applying the chain rule (automatic differentiation in forward mode), but it requires one pass per parameter. In a network there is a single scalar loss and very many parameters: the reverse mode computes the derivative of one output with respect to all the inputs in a single pass. It is a case of choosing the order of multiplication that makes the computation cheap.

Does the brain use backpropagation? It is an open research question, and it is not necessary to solve it to use the algorithm. There is no simple evidence that the brain propagates exact gradients backward through the synapses in the way described here; backpropagation should be considered an engineering method and not a model of the brain.

Do I have to implement it by hand every time? No. In daily work you use the automatic differentiation of PyTorch or TensorFlow. Knowing how to write it helps you understand what the library does, diagnose problems (vanishing gradients, dead ReLUs, wrong shapes) and read other people's code with a critical eye.

What you can do now

You can derive backpropagation with the chain rule, apply it by hand to a neuron and to a 2 → 2 → 1 network, and implement it for a network with any number of layers, checking the shape of every gradient. You can verify the result with finite differences and recognize from the recurrence why the gradient vanishes with the sigmoid: in the ten-layer experiment the norm of the first layer was about six hundred thousand times smaller than that of the last. The reference lessons are layered networks and the forward pass and gradient descent; in the next one, a neural network from scratch in NumPy, you put it all together on real data.

Other lessons in the series that link to this one: Calculus 1: Limits and Derivatives, Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch, The perceptron and activation functions: what non-linearity is for, Convolution: filters, kernels and why look locally, Recurrent Neural Networks (RNN): Giving the Network a Memory, Backpropagation Through Time: Vanishing and Exploding Gradients, Attention, Derived in Full: Query, Key, Value and the Transformer Block Written From Scratch.