5  Backpropagation

Recall the update rule that ended Chapter 1: new parameters come from old parameters by stepping against the gradient, \(\vect{w}^{(t+1)} = \vect{w}^{(t)} - \eta\,\nabla_{\vect{w}}\loss\). For linear regression we could write that gradient down by hand. But a modern network has millions, sometimes billions, of knobs, tangled through layer after layer of composed functions \(\rightarrow\) computing each partial derivative naively, one at a time, is hopeless.

Backpropagation is the algorithm that computes all of them, exactly, in one backward sweep whose cost is a small constant multiple of the graph’s elementary operations. It is not a new kind of calculus. It is the chain rule, organized. By the end of this chapter you will have derived it, implemented it by hand, rebuilt it as a tiny scalar autograd engine, and then checked that torch.autograd agrees with you to machine precision.

This chapter owns the mechanics. For reverse accumulation as a memory and numerical diagnostic—including fan-out accumulation, vector–Jacobian products, retained state, and checkpoint replay—continue in The Gradient Has a Memory.

5.1 One neuron, one chain

Start with the smallest possible case: one neuron on one training example. The previous activation \(a^{(l-1)}\) is multiplied by a weight \(w\), giving the pre-activation \(z = w\,a^{(l-1)} + b\); the activation function produces \(a = \sigma(z)\); and \(a\) feeds a loss \(L\).

Now turn the knob \(w\) a little. How much does \(L\) change? Follow the chain: turning \(w\) changes \(z\); the change in \(z\) changes \(a\); the change in \(a\) changes \(L\). The chain rule multiplies these local sensitivities:

\[ \frac{\partial L}{\partial w} = \frac{\partial L}{\partial a}\cdot \frac{\partial a}{\partial z}\cdot \frac{\partial z}{\partial w} . \tag{5.1}\]

Figure 5.1: Computation as a graph: the forward pass (top, blue) computes and caches values; the backward pass (bottom, orange) multiplies local derivatives along the same edges, in reverse.

Two things about this picture generalize to any network, however deep:

1. The forward pass caches. Computing \(L\) produces the intermediate values (\(z\), \(a\)) along the way; we keep them, because the backward pass will need them.
2. The backward pass reuses. Each edge contributes one local derivative, and the derivative of \(L\) with respect to anything is the product of local derivatives along the path back to it. Nothing is recomputed \(\rightarrow\) the whole sweep costs about as much as one forward pass.

5.2 The game of blame

For a full layer of neurons the bookkeeping threatens to become a tangled mess of sums. The cure is to pick the right intermediate quantity and track only that. At the output of the network the error is obvious: it is just the gap to the target. What is not obvious is how much of that error is the fault of a weight buried ten layers deep. Backpropagation sends this blame backward, from where it is obvious to where it is not.

NoteThe blame signal

Define, for each layer \(l\), the blame signal

\[ \delta^{(l)} \;:=\; \frac{\partial L}{\partial \vect{z}^{(l)}} , \]

the sensitivity of the loss to that layer’s pre-activations. It answers: how much did each neuron’s \(z\) contribute to the final loss?

Everything now unfolds in four short equations. With \(\odot\) denoting elementwise (Hadamard) product:

1. Start where blame is obvious (output layer \(L\)):

\[ \delta^{(L)} = \frac{\partial L}{\partial \vect{a}^{(L)}} \odot \sigma'(\vect{z}^{(L)}) . \tag{5.2}\]

For the half-squared-error convention \(L=\tfrac12\|\vect{a}^{(L)}-\vect{y}\|^2\), \(\partial L / \partial \vect{a}^{(L)}\) is exactly the prediction error \(\vect{a}^{(L)} - \vect{y}\). Other conventions contribute a known scale: mean squared error, for example, contributes \(2/n\) for \(n\) scalar outputs. In every case the output error is then gated by how sensitive the activation was to its input.

2. Propagate blame one layer back:

\[ \delta^{(l)} = \Bigl( (\matr{W}^{(l+1)})^\top \delta^{(l+1)} \Bigr) \odot \sigma'(\vect{z}^{(l)}) . \tag{5.3}\]

Blame flows backward through the same weights the signal flowed forward through, hence the transpose; then the local gate \(\sigma'(\vect{z}^{(l)})\) scales it. This recursion runs from the last layer to the first.

1. Prepare the inputs and fixed settings for the example.
2. Compute the backward blame signal through every layer.

import torch
import matplotlib.pyplot as plt

# [1]
torch.manual_seed(6050)
sizes = [3, 4, 4, 1]
Ws = [torch.randn(sizes[i + 1], sizes[i]) for i in range(3)]

x = torch.randn(3)                       # forward pass, cached
zs, a = [], x
# [2]
for W in Ws[:-1]:                        # hidden layers: sigmoid
    zs.append(W @ a)
    a = torch.sigmoid(zs[-1])
zs.append(Ws[-1] @ a)                    # identity output, as in our 2-layer build

delta = [None] * 3                       # backward pass: eqs. 1 and 2
sp = lambda z: torch.sigmoid(z) * (1 - torch.sigmoid(z))
delta[2] = zs[2] - 1.0                               # eq. 1 (identity output, y = 1)
delta[1] = (Ws[2].T @ delta[2]) * sp(zs[1])          # eq. 2
delta[0] = (Ws[1].T @ delta[1]) * sp(zs[0])          # eq. 2 again

Figure 5.2: Equations 1 and 2 in action on a real network (sizes 3–4–4–1, random weights, one example). Node size and color show each neuron’s blame \(|\delta_j|\): it starts at the output, where the error is obvious, and flows backward through \(W^\top\), shrinking wherever a sigmoid gate is nearly closed.

3. Convert blame into weight gradients:

\[ \frac{\partial L}{\partial \matr{W}^{(l)}} = \delta^{(l)} \, (\vect{a}^{(l-1)})^\top . \tag{5.4}\]

4. And into bias gradients:

\[ \frac{\partial L}{\partial \vect{b}^{(l)}} = \delta^{(l)} . \tag{5.5}\]

Equation 5.4 deserves a pause, because its shapes tell the story. With \(m\) neurons in layer \(l\) and \(n\) in layer \(l-1\): \(\delta^{(l)}\) is \((m \times 1)\), and \((\vect{a}^{(l-1)})^\top\) is \((1 \times n)\). Their outer product is an \(m \times n\) matrix, exactly the shape of \(\matr{W}^{(l)}\), and entry \((i,j)\) is the blame of neuron \(i\) times the activation of neuron \(j\) that fed it. The update for a weight is blame out times signal in:

1. Form the weight-gradient outer product.

# [1]
d2 = delta[1]                              # layer-2 blame        (4,)
a1 = torch.sigmoid(zs[0])                  # layer-1 activations  (4,)
G = torch.outer(d2, a1)                    # eq. 3

Figure 5.3: Equation 3, pictured with the blame and activations from the network above: a \((4 \times 1)\) blame column times a \((1 \times 4)\) activation row yields the full \((4 \times 4)\) gradient matrix in one stroke — entry \((i,j)\) is blame out \(\times\) signal in.

No messy summations: the matrix form performs them all at once, and it maps directly onto the parallel matrix multiplies GPUs are built for. The same four equations govern a toy network and a billion-parameter model; only the matrix sizes change.

5.3 Build it, then check it against autograd

Let us implement the four equations on a two-layer network, with no autograd anywhere, and then ask torch.autograd to differentiate the same network. If our blame algebra is right, the two answers must agree to machine precision.

1. Load the chapter dependencies and establish reproducible state.

import torch

# [1]
_ = torch.manual_seed(6050)

1. Prepare the inputs and fixed settings for the example.
2. Cache the forward pass, propagate blame backward, and form both gradients.
3. Report or visualize the measured result.

# [1]
n, d, h = 32, 5, 8                            # batch, input dim, hidden dim
X = torch.randn(n, d)
y = torch.randn(n)

W1, b1 = torch.randn(h, d) * 0.3, torch.zeros(h)
W2, b2 = torch.randn(1, h) * 0.3, torch.zeros(1)

# forward pass -- cache everything the backward pass will need
Z1 = X @ W1.T + b1                            # (n, h)
A1 = torch.sigmoid(Z1)                        # (n, h)
Z2 = A1 @ W2.T + b2                           # (n, 1)
L = ((Z2.squeeze() - y) ** 2).mean()

# backward pass -- the four equations (batch-averaged)
sig_prime = A1 * (1 - A1)                     # sigma'(z) = sigma(z)(1 - sigma(z))
delta2 = 2 * (Z2.squeeze() - y)[:, None] / n  # eq. 1 (identity output: no gate)
delta1 = (delta2 @ W2) * sig_prime            # eq. 2: blame through W^T, gated
grads = {                                     # eqs. 3 & 4: blame out x signal in
    "W2": delta2.T @ A1, "b2": delta2.sum(0),
    "W1": delta1.T @ X,  "b1": delta1.sum(0),
}

# same network, autograd's turn
params = {k: v.clone().requires_grad_(True) for k, v in
          {"W1": W1, "b1": b1, "W2": W2, "b2": b2}.items()}
A1_ = torch.sigmoid(X @ params["W1"].T + params["b1"])
L_ = (((A1_ @ params["W2"].T + params["b2"]).squeeze() - y) ** 2).mean()
# [2]
L_.backward()

# [3]
for k in grads:
    gap = (grads[k] - params[k].grad).abs().max()
    print(f"{k}: max |ours - autograd| = {gap:.2e}")

W2: max |ours - autograd| = 0.00e+00
b2: max |ours - autograd| = 0.00e+00
W1: max |ours - autograd| = 3.73e-09
b1: max |ours - autograd| = 1.86e-09

The largest discrepancy is \(3.73\times10^{-9}\), consistent with float32 rounding. The blame equations are what autograd computes; the framework simply builds the computation graph for us as we go and then walks it backward, caching and reusing exactly as we did.

5.4 A tiny scalar autograd engine

To remove the last trace of magic, let us build the graph ourselves. Each Value node remembers how it was made (its parents and the local derivative rule); calling backward() walks the graph in reverse topological order, accumulating blame into every node’s grad.

This teaching implementation is adapted from Andrej Karpathy’s micrograd, released under the MIT License. We use only the scalar reverse-mode pattern needed for this derivation; the upstream copyright and license notice are preserved in the repository’s third-party notices.

1. Implement scalar reverse-mode autodiff with a topological backward traversal.

# [1]
class Value:
    """A scalar that remembers its history -- a node in the computation graph."""

    def __init__(self, data: float, parents=(), backward_rule=lambda: None):
        self.data = data
        self.grad = 0.0
        self._parents = parents
        self._backward_rule = backward_rule

    def __add__(self, other: "Value | float") -> "Value":
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data + other.data, (self, other))
        def rule():                      # d(out)/d(self) = d(out)/d(other) = 1
            self.grad += out.grad
            other.grad += out.grad
        out._backward_rule = rule
        return out

    def __mul__(self, other: "Value | float") -> "Value":
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data * other.data, (self, other))
        def rule():                      # product rule, one side at a time
            self.grad += other.data * out.grad
            other.grad += self.data * out.grad
        out._backward_rule = rule
        return out

    def sigmoid(self) -> "Value":
        s = 1 / (1 + 2.718281828459045 ** (-self.data))
        out = Value(s, (self,))
        def rule():                      # sigma' = s (1 - s), the gate
            self.grad += s * (1 - s) * out.grad
        out._backward_rule = rule
        return out

    def backward(self) -> None:
        order, seen = [], set()
        def topo(v):                     # creators land in order first
            if v not in seen:
                seen.add(v)
                for p in v._parents:
                    topo(p)
                order.append(v)
        topo(self)
        self.grad = 1.0                  # dL/dL = 1: blame starts here
        for v in reversed(order):
            v._backward_rule()

1. Prepare the inputs and fixed settings for the example.
2. Verify micro autograd.
3. Report or visualize the measured result.

# [1]
w, x, b = Value(0.7), Value(2.0), Value(-0.5)
a = ((w * x) + b).sigmoid()             # one neuron, forward
loss = (a + (-0.3)) * (a + (-0.3))      # squared error vs target 0.3
# [2]
loss.backward()

# analytic check: dL/dw = 2(a - y) * sigma'(z) * x
z = 0.7 * 2.0 - 0.5
s = 1 / (1 + 2.718281828459045 ** (-z))
# [3]
print(f"micro-autograd: dL/dw = {w.grad:.6f}")
print(f"by hand:        dL/dw = {2 * (s - 0.3) * s * (1 - s) * 2.0:.6f}")

micro-autograd: dL/dw = 0.337801
by hand:        dL/dw = 0.337801

In a few dozen lines, it is the real thing: PyTorch’s autograd is this idea with tensors instead of scalars, a compiled graph walker, and thousands of derivative rules.

5.5 torch.autograd in practice: five rules

Autograd is doing exactly what we just built by hand, and knowing that mechanics tells you why each of its rules exists. These five cover nearly every autograd bug you will meet; the worked cases below isolate each path before the full network combines them.

Rule 1 — by default, only leaf tensors retain gradients. A leaf is a tensor you created directly (your parameters); a non-leaf is the result of an operation (activations, losses). After backward(), leaves have .grad populated while intermediate blame is used and released. Call .retain_grad() on a non-leaf before backward() only when you explicitly need its gradient for inspection.

Rule 2 — gradients accumulate. backward() adds into .grad. Zero them before each step (optimizer.zero_grad()), or every update uses the sum of all past gradients.

1. Prepare the inputs and fixed settings for the example.
2. Backpropagate twice without clearing the leaf gradient.
3. Report or visualize the measured result.

# [1]
w = torch.randn(3, requires_grad=True)      # leaf: our parameter
y_hat = (w * 2).sum()                       # non-leaf: result of an operation

# [2]
y_hat.backward()
# [3]
print(f"leaf grad:      {w.grad.tolist()}")            # populated
print(f"non-leaf grad:  {y_hat.grad}")                 # None -- rule 1

(w * 2).sum().backward()                    # backward again, WITHOUT zeroing
print(f"after 2nd pass: {w.grad.tolist()}")            # doubled -- rule 2

leaf grad:      [2.0, 2.0, 2.0]
non-leaf grad:  None
after 2nd pass: [4.0, 4.0, 4.0]

Rule 3 — parameter updates happen outside the graph. For a raw leaf tensor, the update w = w - lr * w.grad inside autograd’s view records the update and rebinds w to a non-leaf, so its .grad will not be retained by default on the next pass. Assigning a plain tensor over a registered nn.Parameter may instead fail loudly. The safe idiom below covers both cases:

1. Save the leaf parameter and choose a learning rate.
2. Pause graph recording and update the parameter in place.
3. Verify that the update changed the value without changing leaf status.

# [1]
before = w.detach().clone()
lr = 0.1

# [2]
with torch.no_grad():
    w -= lr * w.grad

# [3]
assert not torch.equal(w, before)
assert w.is_leaf

Rule 4 — in-place operations can corrupt the tape. Methods ending in _ (add_, relu_) overwrite values that the backward pass may still need from the forward cache. Inside no_grad() parameter updates they are safe; anywhere else, prefer creating a new tensor.

Rule 5 — the graph lives until you let it go. Every forward pass builds a fresh graph, freed after backward(). But any tensor you keep (say, appending loss to a list for plotting) keeps its whole graph alive with it. Store loss.item() (a plain float) or loss.detach() instead: detach() returns the same value with the history cut.

WarningKeep updates outside the graph

Rules 1 and 3 combine into a classic raw-tensor bug: rebinding w = w - lr * w.grad creates a new non-leaf. Its .grad then comes back None by default on the next pass. A registered nn.Parameter assignment can raise an error instead. In either case, update the existing parameter under torch.no_grad().

5.6 The gate, and the problem with deep chains

Look again at Equation 5.3. Every backward step multiplies by \(\sigma'(\vect{z}^{(l)})\). That derivative acts as a gate on the blame: wide open, and the signal passes; nearly closed, and it dies.

Now examine the sigmoid’s gate. From \(\sigma'(z) = \sigma(z)(1-\sigma(z))\), its maximum value, at \(z = 0\), is exactly \(0.25\); away from zero the neuron saturates and the gate approaches \(0\). The sigmoid gate is at best a quarter open, and often nearly closed.

1. Evaluate the sigmoid and ReLU derivative gates.

import matplotlib.pyplot as plt

# [1]
z = torch.linspace(-6, 6, 300)
sig = torch.sigmoid(z)
sigmoid_gate = sig * (1 - sig)
relu_gate = (z > 0).float()

Figure 5.4: Two gates. The sigmoid’s derivative never exceeds 0.25 and vanishes under saturation; ReLU’s is fully open (exactly 1) for every active neuron.

Here is the consequence, and it is one of deep learning’s most important failure modes. The chain rule multiplies local Jacobians across layers. The activation-derivative part of a sigmoid path contributes at most \(0.25\) per layer, so that part is bounded by \(0.25^k\) after \(k\) gates. Ten gates in, its ceiling is below one in a million. The weight matrices and branching paths also scale the full network Jacobian, so \(0.25^k\) is not a bound on the complete gradient. It isolates why repeated sigmoid gates strongly favor vanishing signals.

The fix is almost embarrassingly simple. \(\mathrm{ReLU}(z) = \max(0, z)\) has derivative \(1\) for \(z > 0\) and \(0\) for \(z < 0\). At the kink \(z=0\), PyTorch uses derivative \(0\); exact zeros can occur with zero biases or sparse activations. The gate is either fully open or fully closed. For every neuron that was active in the forward pass, the activation derivative contributes a factor of exactly one. The surrounding weight matrices can still scale the full signal, but the active ReLU does not shrink it: a gradient superhighway through the network. This, more than anything, is why ReLU and its variants are the default in deep networks.

Do not take the argument on faith; measure it. We build the same 30-layer network twice (sigmoid versus ReLU), give both copies the same He-scaled weights and input batch, and record how much gradient reaches each layer. We compute the diagnostic in float64 so the norm itself does not underflow; then we ask the practical float32 question: would a learning-rate-\(0.1\) update actually change the first layer’s stored weights?

1. Define the reusable gradient_diagnostic helper.
2. Prepare the inputs and fixed settings for the example.
3. Implement the depth experiment.
4. Report or visualize the measured result.

from torch import nn

# [1]
def gradient_diagnostic(
    activation: type[nn.Module], depth: int = 30, width: int = 64
) -> dict[str, object]:
    torch.manual_seed(6050)
    layers = []
    for _ in range(depth):
        layers += [nn.Linear(width, width), activation()]
    net = nn.Sequential(*layers, nn.Linear(width, 1)).double()
    with torch.no_grad():
        for layer in net:
            if isinstance(layer, nn.Linear):
                nn.init.normal_(layer.weight, std=(2 / layer.in_features) ** 0.5)
                nn.init.zeros_(layer.bias)

    out = net(torch.randn(128, width, dtype=torch.float64)).mean()
    out.backward()
    linears = [layer for layer in net if isinstance(layer, nn.Linear)][:-1]
    first = linears[0]
    grad32 = first.weight.grad.float()
    weight32 = first.weight.detach().float()
    changed = ((weight32 - 0.1 * grad32) != weight32).sum().item()
    return {
        "norms": [layer.weight.grad.norm().item() for layer in linears],
        "first_max": first.weight.grad.abs().max().item(),
        "first_nonzero": grad32.count_nonzero().item(),
        "first_changed": changed,
        "first_size": grad32.numel(),
    }

# [2]
diagnostics = {}
plt.figure(figsize=(6.2, 3.4))
# [3]
for act, color in [(nn.Sigmoid, "#232D4B"), (nn.ReLU, "#E57200")]:
    diagnostics[act.__name__] = gradient_diagnostic(act)
    plt.semilogy(diagnostics[act.__name__]["norms"], "o-", ms=3,
                 color=color, label=act.__name__)
plt.xlabel("layer (0 = closest to input)")
plt.ylabel(r"$\|\partial L / \partial W^{(l)}\|$")
# [4]
plt.legend(); plt.tight_layout(); plt.show()

for name, result in diagnostics.items():
    print(f"{name:7s}: max|g0|={result['first_max']:.2e}; "
          f"nonzero={result['first_nonzero']}/{result['first_size']}; "
          f"changed={result['first_changed']}/{result['first_size']}")

Figure 5.5: Gradient reaching each layer of two matched 30-layer MLPs (float64 norm, log scale). With the same He-scaled weights and inputs, repeated sigmoid gates leave a nonzero but unusably small first-layer signal; ReLU preserves enough signal for a float32 update to change every first-layer weight.

Sigmoid: max|g0|=9.47e-17; nonzero=4096/4096; changed=0/4096
ReLU   : max|g0|=2.54e-02; nonzero=4096/4096; changed=4096/4096

Read the slopes, then read the update audit. The sigmoid first-layer entries are not mathematically or representationally zero: all 4,096 survive the float32 cast. They are so small, however, that subtracting \(0.1g\) changes none of the stored float32 weights. The ReLU copy changes all 4,096. That distinction matters: a displayed float32 norm can underflow while individual entries remain nonzero, and nonzero gradients can still be too small to move a parameter at its current magnitude. ReLU is not magic, but its open gates make the signal usable in this matched experiment. Appendix C names these as separate representation failures: a nonzero value can survive near zero while the same update rounds away against a larger stored weight.

5.7 Initialization: setting the signal’s energy

One more lever lives in this chapter, because it falls out of the same variance bookkeeping. Think of the network as a pipeline and the variance of the signal as its energy. If each layer multiplies the variance by more than one, activations explode after enough layers; by less than one, they vanish, and by Equation 5.3 the gradients follow. A good initialization keeps the energy roughly constant in both directions.

The analysis is one line. For \(z = \sum_{i=1}^{n_{\text{in}}} w_i x_i\) with independent, zero-mean inputs and weights, \(\var(z) = n_{\text{in}} \var(w) \var(x)\). Keeping \(\var(z) = \var(x)\) demands

\[ \var(w) = \frac{1}{n_{\text{in}}} \quad \text{(forward)}, \qquad \var(w) = \frac{1}{n_{\text{out}}} \quad \text{(backward)}. \]

Xavier initialization splits the difference for symmetric activations like tanh and sigmoid, \(\var(w) = 2/(n_{\text{in}} + n_{\text{out}})\). He initialization accounts for ReLU discarding half its input distribution by doubling the forward recipe, \(\var(w) = 2/n_{\text{in}}\). PyTorch applies sensible defaults for you; the difference between a good and a careless choice is smooth training versus a dead or exploding network:

1. Define the reusable signal_std helper.
2. Prepare the inputs and fixed settings for the example.
3. Signal energy vs depth.
4. Report or visualize the measured result.

from collections.abc import Callable

# [1]
def signal_std(
    scale_fn: Callable[[torch.Tensor, int], torch.Tensor],
    depth: int = 40, width: int = 256
) -> list[float]:
    x = torch.randn(512, width)
    stds = []
    for _ in range(depth):
        W = scale_fn(torch.randn(width, width), width)
        x = torch.relu(x @ W)
        stds.append(x.std().item())
    return stds

# [2]
torch.manual_seed(6050)
schemes = {
    "naive  (std = 0.05)": lambda W, n: W * 0.05,
    "naive  (std = 0.12)": lambda W, n: W * 0.12,
    "He  (std = sqrt(2/n))": lambda W, n: W * (2 / n) ** 0.5,
}
plt.figure(figsize=(6.2, 3.4))
# [3]
for (name, fn), color in zip(schemes.items(), ["#5379AA", "#722F37", "#E57200"]):
    plt.semilogy(signal_std(fn), color=color, label=name)
plt.xlabel("layer"); plt.ylabel("std of activations")
# [4]
plt.legend(); plt.tight_layout(); plt.show()

Figure 5.6: Signal energy through 40 ReLU layers. Careless scales lose or amplify the signal exponentially; He initialization holds it steady.

Initialization gives a stable start. Keeping the signal stable during training, as the weights move, is the job of normalization layers; and for very deep networks even that is not enough, and gradients get an architectural bypass: skip connections that add a block’s input straight to its output. Both belong to later chapters (normalization and the residual idea in Chapter 9, and layer normalization where it becomes essential, inside the Transformer of Chapter 14). What you should carry from here is the invariant they all serve: keep the signal, forward and backward, at constant energy.

NoteCheck yourself

Close the book for one minute and run the chain rule backward.

• What does each layer cache on the forward pass?
• Why does the blame signal pass through a transposed weight matrix?
• How do activation derivatives and initialization together control gradient scale with depth?

5.8 Okay, so — the chain rule, organized

1. The graph. Any network is a computation graph; the forward pass computes and caches, the backward pass multiplies local derivatives in reverse.
2. The blame signal. \(\delta^{(l)} = \partial L / \partial \vect{z}^{(l)}\) turns the tangle into four equations: start at the output (Equation 5.2), propagate through \(\matr{W}^\top\) and the gate (Equation 5.3), then read off weight and bias gradients (Equation 5.4, Equation 5.5). The backward sweep costs a small constant multiple of the forward graph’s operations, for every gradient at once.
3. We verified it: our hand-built equations match torch.autograd to \(10^{-8}\), and a tiny scalar Value class reproduces the machinery end to end.
4. The gate rules the depth. The sigmoid-derivative product along a path is bounded above by \(0.25^k\) after \(k\) gates; ReLU’s open gate is a gradient superhighway. Initialization sets the signal’s energy so neither direction explodes or dies at the start.

Backpropagation is the engine under every later chapter, through Chapter 20. The next time a deep model fails to learn, your first questions now have names: is the gradient reaching the early layers? what is gating it? what is its energy?

Sources and further reading

• Rumelhart, Hinton, and Williams, “Learning representations by back-propagating errors” (1986): the landmark short account of learning internal representations by reverse error propagation.
• Baydin et al., “Automatic Differentiation in Machine Learning: a Survey” (2018): separates reverse-mode automatic differentiation from symbolic and numerical differentiation.
• Glorot and Bengio, “Understanding the difficulty of training deep feedforward neural networks” (2010): the variance argument behind Xavier initialization.
• He et al., “Delving Deep into Rectifiers” (2015): derives the rectifier-aware initialization used in the matched depth experiment.
• CS231n, “Gradient Checks”: the centered-difference and relative-error safeguards used in Exercise 6.
• Ioffe and Szegedy, “Batch Normalization” (2015): the coupled minibatch normalization whose backward pass is audited in Exercise 7.

Exercises

1. (Pencil.) Derive the four blame equations for a network whose output layer uses the identity activation and squared-error loss (the network of our verification cell). Confirm that your \(\delta^{(2)}\) matches the delta2 line in the code.
2. (Pencil.) In Equation 5.3, why does the blame flow through \((\matr{W}^{(l+1)})^\top\) rather than \(\matr{W}^{(l+1)}\)? Argue from shapes: what are the dimensions of \(\delta^{(l)}\), \(\delta^{(l+1)}\), and \(\matr{W}^{(l+1)}\)?
3. (Code.) Extend the Value class with tanh and use it to train the one-neuron model from the check cell for 50 steps (write the update loop yourself). Verify the loss decreases.
4. (Code.) In the depth experiment, give the sigmoid network Xavier initialization (nn.init.xavier_normal_ on each linear layer). How much of the vanishing does good initialization recover at depth 30? What does that tell you about why both fixes, initialization and ReLU, became standard?
5. (Code.) Comment out optimizer.zero_grad() in any training loop from Chapter 1 and describe what happens to the loss curve. Explain it with one sentence about .grad accumulation.
6. (Pencil.) Derive why a forward finite difference has truncation error \(O(\varepsilon)\) while a centered difference has error \(O(\varepsilon^2)\). (Audit.) Gradient-check this chapter’s two-layer manual backward pass with \(|g_a-g_n|/\max(10^{-12},|g_a|,|g_n|)\) for \(\varepsilon=10^{-2},\ldots,10^{-10}\) in float32 and float64. Plot relative error, explain its U-shape as truncation meets rounding, and require agreement across several randomly selected parameters under seed 6050. Named wrong answer: “absolute error below \(10^{-4}\) proves the gradient” — absolute error ignores the scale of the quantity being checked.
7. (Pencil.) Derive \(\partial\mathcal{L}/\partial x\) through BatchNorm’s batch mean and variance. Mark the direct normalized-input term and the two terms created by coupling every example through those statistics. (Code.) Verify all three terms against autograd on the experimentation interlude’s BatchNorm module. Then delete the two coupling terms deliberately, compare the gradient error, and rerun the interlude’s training curve under the same data, initialization, minibatch schedule, and seed panel.