Backpropagation: Chain Rule, Gradient Flow, and Autograd

How does a neural network learn? In the previous chapters, we built multilayer perceptrons that transform inputs through layers of weights and activations, and we defined loss functions that measure prediction quality. But we haven't answered the fundamental question: how do we update the millions of parameters in a neural network to reduce the loss?

The answer is backpropagation, short for "backward propagation of errors." This algorithm computes the gradient of the loss function with respect to every weight in the network, telling us exactly how to adjust each parameter to improve predictions. Backpropagation is the most important algorithm in deep learning. Without it, training neural networks would be computationally intractable.

This chapter builds backpropagation from the ground up. We start with computational graphs that visualize how operations compose, then derive the chain rule that makes gradient computation tractable. You'll trace through forward and backward passes step by step, understand gradient flow through activation functions, examine Jacobian matrices for vector operations, and implement backpropagation from scratch. We'll also cover numerical gradient checking to verify implementations, common backprop bugs, and how PyTorch autograd works internally. By the end, you'll see exactly how modern deep learning frameworks compute gradients automatically.

Neural networks have millions of parameters. GPT-3 has 175 billion. How do we figure out which direction to nudge each one to reduce the loss?

Calculus gives us the answer: the gradient. For a function $f(\mathbf{w})$ where $\mathbf{w}$ is a vector of parameters, the gradient $\nabla_\mathbf{w} f$ is a vector that points in the direction of steepest increase of $f$. Each component tells us how sensitive the function is to changes in the corresponding parameter. To minimize the loss, we move opposite to the gradient:

$$\mathbf{w}_{\text{new}} = \mathbf{w}_{\text{old}} - \eta \nabla_\mathbf{w} \mathcal{L}$$

where:

• $\mathbf{w}_{\text{new}}$: the updated parameter vector after one step
• $\mathbf{w}_{\text{old}}$: the current parameter vector before the update
• $\eta$: the learning rate, a small positive scalar controlling step size
• $\nabla_\mathbf{w} \mathcal{L}$: the gradient of the loss with respect to parameters, where each entry $\frac{\partial \mathcal{L}}{\partial w_i}$ indicates how much the loss increases if we increase $w_i$ by a tiny amount
• $\mathcal{L}$: the loss function, a scalar value measuring prediction error

The subtraction makes intuitive sense: if $\frac{\partial \mathcal{L}}{\partial w_i} > 0$, increasing $w_i$ would increase the loss, so we decrease it instead. This is gradient descent, simple in concept, but the challenge is computing $\nabla_\mathbf{w} \mathcal{L}$ when the loss depends on parameters through many nested operations.

Consider a 3-layer network. The loss depends on the output, which depends on the third layer activations, which depend on the second layer activations, which depend on the first layer activations, which depend on the input and first layer weights. This chain of dependencies seems hopelessly complex to differentiate.

Backpropagation solves this elegantly. By applying the chain rule in reverse order, from output back to input, we compute all gradients in time proportional to a single forward pass. This efficiency is what makes training deep networks practical.

The story of backpropagation stretches back to the 1960s and 1970s, though its impact on neural networks only became clear in the 1980s. The basic idea of differentiating through a chain of composed functions was understood in the automatic differentiation community, but applying it to neural networks required several conceptual steps.

The paper that really crystallized the algorithm for the machine learning community was Rumelhart, Hinton, and Williams' 1986 paper "Learning representations by back-propagating errors" in Nature. This paper demonstrated that backpropagation could efficiently train multilayer networks to learn useful internal representations, which was not obvious at the time. Before this, the dominant view in AI was that useful internal representations needed to be hand-designed by experts.

The key insight of the 1986 paper was simple but powerful: hidden units in a network could learn to represent features of the input that were useful for the task, and backpropagation would discover what those features should be automatically. This was a striking departure from prior approaches that either used only one-layer models or required laborious feature engineering.

Despite this theoretical breakthrough, neural networks fell out of fashion in the 1990s for two reasons. First, training was slow on the hardware available then. Second, networks large enough to solve interesting problems ran into the vanishing gradient problem, making deep networks hard to train in practice. Other methods like support vector machines and random forests worked better for many problems. The deep learning renaissance beginning around 2006, with better hardware (GPUs), better initialization, and better architectures (like ReLU activations and skip connections), finally made the potential suggested by the 1986 paper practical.

Before diving into backpropagation, we need a way to visualize how computations flow through a network. A computational graph represents a mathematical expression as a directed graph where nodes are operations and edges show data flow.

Every neural network computation breaks down into elementary operations: addition, multiplication, matrix products, and activation functions. Consider the simple expression $f(x, y, z) = (x + y) \cdot z$. This decomposes into two operations:

1. $a = x + y$ (addition)
2. $f = a \cdot z$ (multiplication)

where:

• $x, y, z$: input variables (leaf nodes in the graph)
• $a$: an intermediate value storing the result of the first operation
• $f$: the final output (the root node)

The graph shows data flowing left to right through intermediate nodes. Each node computes something simple, and the structure of the graph encodes the data dependencies. Each node only needs to know about its immediate inputs and output, not the rest of the network. This locality is what makes the chain rule tractable.

A neural network layer performs the computation $\mathbf{a} = \sigma(\mathbf{W}\mathbf{x} + \mathbf{b})$, where:

• $\mathbf{x}$: the input vector
• $\mathbf{W}$: the weight matrix
• $\mathbf{b}$: the bias vector
• $\sigma$: the activation function
• $\mathbf{a}$: the output activations

This breaks into three atomic operations:

1. $\mathbf{z}_1 = \mathbf{W}\mathbf{x}$ (matrix multiplication)
2. $\mathbf{z}_2 = \mathbf{z}_1 + \mathbf{b}$ (bias addition)
3. $\mathbf{a} = \sigma(\mathbf{z}_2)$ (nonlinear activation)

Each node has a well-defined local gradient. The power of computational graphs is that they make the structure of computation explicit, which is exactly what backpropagation needs. The forward pass computes values flowing left to right; the backward pass computes gradients flowing right to left.

When we compose multiple layers, we are simply building a larger graph. The loss node at the end is connected to the output layer, which is connected to the layer before it, and so on back to the input. The entire network is one big directed acyclic graph of operations, with parameters as leaf nodes and the loss as the root.

Computational graphs provide two benefits for gradient computation:

• Modularity: Each node has a well-defined local gradient. We don't need to derive the gradient of the entire network at once, only the gradient of each primitive operation. This is the key insight that makes backprop practical: instead of deriving one giant formula for the gradient of the loss with respect to every parameter, we simply need to know the derivative of each individual operation, and the chain rule handles the composition.
• Efficiency: By storing intermediate values during the forward pass, we reuse them during the backward pass, avoiding redundant computation. Without this storage, we would have to recompute every intermediate value every time we needed it, which would make the backward pass much more expensive.

Deep learning frameworks like PyTorch and TensorFlow build exactly this kind of dynamic computational graph as code runs. Every arithmetic operation on a tensor implicitly creates a node in the graph, recording what inputs were used so the backward pass can reconstruct the gradient computation. You never see this graph directly, but it is being built silently behind every +, @, and activation function call you write.

The chain rule from calculus is the mathematical engine powering backpropagation. It tells us how to differentiate composed functions, and since a neural network is one giant composition of functions, the chain rule is precisely what we need.

Let $u = g(x)$ and $y = f(u) = f(g(x))$. The chain rule states:

$$\frac{dy}{dx} = \frac{dy}{du} \cdot \frac{du}{dx}$$

where:

• $\frac{dy}{dx}$: the total derivative of the final output with respect to the original input
• $\frac{dy}{du}$: the local gradient of the outer function $f$ at the intermediate value $u$
• $\frac{du}{dx}$: the local gradient of the inner function $g$ at $x$

The intuition is multiplicative sensitivity: if $x$ changes by a small amount $\delta x$, it causes $u$ to change by $\frac{du}{dx} \cdot \delta x$, which in turn causes $y$ to change by $\frac{dy}{du} \cdot \frac{du}{dx} \cdot \delta x$. The total sensitivity is the product of each stage's sensitivity.

Think of it as a chain of dominoes. A small push on the first domino ($\delta x$) gets amplified or attenuated at each stage. The chain rule tells us the total amplification from one end of the chain to the other, which is just the product of each stage's individual amplification. This multiplicative structure is central to understanding both why backpropagation works and why deep networks can suffer from vanishing or exploding gradients.

For a chain of three functions $y = f(g(h(x)))$, the chain rule extends naturally:

$$\frac{dy}{dx} = \frac{dy}{df} \cdot \frac{df}{dg} \cdot \frac{dg}{dx}$$

This pattern extends to any depth. For a network with $n$ layers, the gradient of the loss with respect to the first layer's weights involves multiplying $n$ local gradients together. This is why deep networks can suffer from vanishing or exploding gradients: if each local gradient is slightly less than 1, their product shrinks exponentially; if each is slightly greater than 1, the product explodes.

Concretely, suppose we have a 20-layer network and the average gradient magnitude at each layer is 0.9. The gradient at the first layer would be multiplied by $0.9^{20} \approx 0.12$, reducing it to about 12% of the output gradient. That might seem manageable, but with 50 layers and an average gradient of 0.7, we get $0.7^{50} \approx 1.3 \times 10^{-8}$: essentially zero. The first layers of the network would receive essentially no learning signal.

Neural networks work with vectors and matrices, so we need the multivariate form. If $\mathbf{z} = \mathbf{f}(\mathbf{y})$ and $\mathbf{y} = \mathbf{g}(\mathbf{x})$, the chain rule becomes:

$$\frac{\partial z_i}{\partial x_k} = \sum_j \frac{\partial z_i}{\partial y_j} \cdot \frac{\partial y_j}{\partial x_k}$$

where:

• $\frac{\partial z_i}{\partial x_k}$: the partial derivative of the $i$-th output component with respect to the $k$-th input component
• $\frac{\partial z_i}{\partial y_j}$: the partial derivative of $z_i$ with respect to the $j$-th intermediate variable
• $\frac{\partial y_j}{\partial x_k}$: the partial derivative of $y_j$ with respect to the $k$-th input

In matrix form, this is exactly a matrix product: if $\mathbf{J}_f$ and $\mathbf{J}_g$ are the Jacobian matrices of $f$ and $g$ respectively, then the Jacobian of the composition is:

$$\mathbf{J}_{f \circ g} = \mathbf{J}_f \cdot \mathbf{J}_g$$

This matrix multiplication is the vectorized form of the chain rule, and it is what makes backpropagation efficient. The sum over $j$ in the scalar form becomes a dot product, and repeating it for all $(i, k)$ pairs gives a matrix product. Gradient computation can therefore be done with the same fast matrix operations that power the forward pass.

Before diving into the full backward pass, it's worth understanding Jacobian matrices carefully, because they formalize how gradients flow through vector-valued operations.

In neural network backpropagation, the loss $\mathcal{L}$ is a scalar, and each operation in the network transforms a vector to another vector. At each node, we don't technically need the full Jacobian. We only need the vector-Jacobian product (VJP), which maps the upstream gradient vector to a downstream gradient vector.

To see why the full Jacobian is unnecessary, consider a layer that maps $\mathbb{R}^{512} \to \mathbb{R}^{512}$. The full Jacobian would be a $512 \times 512$ matrix of 262,144 entries. But we only need to propagate a gradient vector of length 512 backward. The VJP computes this directly in $O(512^2)$ time, the same cost as forming the Jacobian, but without ever explicitly materializing it. For larger layers, the savings are more dramatic, and specialized operations like batch normalization have structured Jacobians that admit even faster VJP computations.

In practice, backpropagation never explicitly forms full Jacobian matrices. Instead, it computes vector-Jacobian products. If the upstream gradient is a row vector $\mathbf{v}^T$ (with the same shape as the output), the downstream gradient is:

$$\mathbf{v}^T \cdot \mathbf{J}$$

This produces a row vector with the same shape as the input. The key insight is that computing this product is much cheaper than computing $\mathbf{J}$ explicitly, then multiplying. PyTorch's autograd system is built entirely around this VJP abstraction.

This VJP perspective also clarifies what backpropagation is computing at each step. The upstream gradient $\mathbf{v}$ represents "how the loss changes when the output of this node changes." The Jacobian $\mathbf{J}$ represents "how the output of this node changes when the input changes." Their product gives "how the loss changes when the input changes," which is exactly the downstream gradient we need to propagate.

The power of computational graphs is that each operation only needs to define its own local Jacobian (or VJP). The full gradient is then assembled by chaining these local gradients via the multivariate chain rule.

Here are the local gradients for the key operations in a neural network:

Addition: $\mathbf{z} = \mathbf{x} + \mathbf{y}$. The Jacobian with respect to each input is the identity matrix $\mathbf{I}$. Gradients pass through unchanged: $\frac{\partial \mathcal{L}}{\partial \mathbf{x}} = \frac{\partial \mathcal{L}}{\partial \mathbf{z}}$ and $\frac{\partial \mathcal{L}}{\partial \mathbf{y}} = \frac{\partial \mathcal{L}}{\partial \mathbf{z}}$. This is why residual connections are so effective: the skip path provides a direct gradient highway back to earlier layers. The gradient flows through without any attenuation from the addition node itself.

Multiplication by scalar $w$: $\mathbf{z} = w \cdot \mathbf{x}$. The gradient with respect to $\mathbf{x}$ is $w \cdot \frac{\partial \mathcal{L}}{\partial \mathbf{z}}$ (scale up or down), and the gradient with respect to $w$ is $\mathbf{x}^T \frac{\partial \mathcal{L}}{\partial \mathbf{z}}$.

Matrix multiplication: $\mathbf{z} = \mathbf{W}\mathbf{x}$. The Jacobian with respect to $\mathbf{x}$ is $\mathbf{W}^T$, giving $\frac{\partial \mathcal{L}}{\partial \mathbf{x}} = \mathbf{W}^T \frac{\partial \mathcal{L}}{\partial \mathbf{z}}$. The gradient with respect to $\mathbf{W}$ is the outer product $\frac{\partial \mathcal{L}}{\partial \mathbf{z}} \mathbf{x}^T$. Notice the transpose: gradients flow backward through the weight matrix transposed, which is why the backward pass through a linear layer involves $\mathbf{W}^T$ while the forward pass involves $\mathbf{W}$.

Elementwise nonlinearity: $\mathbf{a} = \sigma(\mathbf{z})$. The Jacobian is a diagonal matrix with $\sigma'(z_i)$ on the diagonal, so $\frac{\partial \mathcal{L}}{\partial \mathbf{z}} = \frac{\partial \mathcal{L}}{\partial \mathbf{a}} \odot \sigma'(\mathbf{z})$, where $\odot$ denotes elementwise multiplication. Because the Jacobian is diagonal, we never form a full matrix; we just do an elementwise multiply.

With computational graphs and the chain rule in hand, we can now describe the complete backpropagation algorithm. The algorithm has two phases: a forward pass that computes outputs and stores intermediate values, and a backward pass that computes gradients using those stored values.

During the forward pass, we evaluate the network layer by layer, storing all intermediate values. These stored values are critical because the backward pass needs them to compute local gradients.

For a 2-layer network with input $\mathbf{x}$, weights $\mathbf{W}^{(1)}$, $\mathbf{W}^{(2)}$, biases $\mathbf{b}^{(1)}$, $\mathbf{b}^{(2)}$, and activation $\sigma$, the forward pass computes:

$$\begin{aligned} \mathbf{z}^{(1)} &= \mathbf{W}^{(1)} \mathbf{x} + \mathbf{b}^{(1)} \\ \mathbf{a}^{(1)} &= \sigma(\mathbf{z}^{(1)}) \\ \mathbf{z}^{(2)} &= \mathbf{W}^{(2)} \mathbf{a}^{(1)} + \mathbf{b}^{(2)} \\ \hat{\mathbf{y}} &= \text{softmax}(\mathbf{z}^{(2)}) \\ \mathcal{L} &= -\sum_k y_k \log \hat{y}_k \end{aligned}$$

where:

• $\mathbf{z}^{(l)}$: the pre-activation values at layer $l$ (also called logits)
• $\mathbf{a}^{(l)}$: the post-activation values at layer $l$
• $\hat{\mathbf{y}}$: the predicted probability distribution
• $\mathbf{y}$: the one-hot encoded target
• $\mathcal{L}$: the cross-entropy loss (scalar)

After the forward pass, we have stored $\mathbf{z}^{(1)}, \mathbf{a}^{(1)}, \mathbf{z}^{(2)}, \hat{\mathbf{y}}$. These are the "cache" that the backward pass will consume. The reason we store $\mathbf{z}^{(l)}$ specifically (and not just $\mathbf{a}^{(l)}$) is that the activation derivative $\sigma'(\mathbf{z}^{(l)})$ needs the pre-activation value, not the post-activation value. For sigmoid, though, we can recover $\sigma'(z)$ from the activation value $a$ itself using $\sigma'(z) = a(1-a)$, so storing either is fine.

The backward pass applies the chain rule in reverse, starting from the loss $\mathcal{L}$ and working backward through each operation.

Gradient of loss with respect to output. We start with $\frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}}$. For softmax + cross-entropy, this has the elegant form:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} = \hat{\mathbf{y}} - \mathbf{y}$$

This says: the gradient at the output layer is simply the difference between prediction and target. Larger errors produce larger gradients. This elegant form is not a coincidence: the cross-entropy loss was designed to pair with softmax precisely because the combination yields this clean gradient formula. This pairing is sometimes called the "softmax-cross-entropy trick," and it avoids numerical instability that would arise from computing the softmax and cross-entropy gradients separately.

Gradient through the second layer. Using the matrix multiply gradient rules:

$$\begin{aligned} \frac{\partial \mathcal{L}}{\partial \mathbf{W}^{(2)}} &= \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} \cdot (\mathbf{a}^{(1)})^T \\ \frac{\partial \mathcal{L}}{\partial \mathbf{b}^{(2)}} &= \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} \\ \frac{\partial \mathcal{L}}{\partial \mathbf{a}^{(1)}} &= (\mathbf{W}^{(2)})^T \cdot \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} \end{aligned}$$

Gradient through the activation function. Passing through $\sigma$:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(1)}} = \frac{\partial \mathcal{L}}{\partial \mathbf{a}^{(1)}} \odot \sigma'(\mathbf{z}^{(1)})$$

where $\sigma'(\mathbf{z}^{(1)})$ is the elementwise derivative of the activation function evaluated at the stored pre-activation values.

Gradient through the first layer:

$$\begin{aligned} \frac{\partial \mathcal{L}}{\partial \mathbf{W}^{(1)}} &= \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(1)}} \cdot \mathbf{x}^T \\ \frac{\partial \mathcal{L}}{\partial \mathbf{b}^{(1)}} &= \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(1)}} \end{aligned}$$

The pattern is clear: at each layer, we multiply the upstream gradient by the transpose of each weight matrix to flow backward through layers, then multiply elementwise by the activation derivative at the nonlinearity, then pass the result further back. The upstream gradient arrives from the right (closer to the loss), we compute local contributions for this layer's parameters, and we send a modified gradient leftward to the previous layer.

The backward pass depends critically on the derivatives of activation functions. Different activations have very different gradient behaviors, which directly affects how well information flows backward through deep networks. We covered activation functions in depth in the Activation Functions chapter, so here we focus on their gradient properties specifically in the context of backpropagation.

The sigmoid function $\sigma(z) = \frac{1}{1+e^{-z}}$ has derivative:

$$\sigma'(z) = \sigma(z)(1 - \sigma(z))$$

where:

• $\sigma(z)$: the sigmoid output, a value in $(0, 1)$
• $\sigma'(z)$: the derivative, which equals $\sigma(z)(1-\sigma(z))$

The maximum value of $\sigma'(z)$ is $0.25$, achieved at $z = 0$. For $|z|$ large, the sigmoid saturates (output near 0 or 1), and the derivative approaches 0. This vanishing gradient problem is severe: each sigmoid layer multiplies the gradient by at most 0.25. In a 10-layer network, the gradient at the first layer could be multiplied by $0.25^{10} \approx 10^{-6}$ relative to the last layer.

There is another way to see why this is so problematic. The gradient of the loss with respect to the first layer's weights is computed as a product of gradients across all intermediate layers. If those gradients are small fractions (as they are for sigmoid in saturation), their product is tiny. This means the first-layer weights barely move during training. The network effectively cannot learn any useful feature representation in its early layers when trained with sigmoid activations in deep networks.

The hyperbolic tangent $\tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}}$ has derivative:

$$\tanh'(z) = 1 - \tanh^2(z)$$

This has maximum value 1.0 at $z = 0$ and decays toward 0 as $|z|$ grows. Tanh is better than sigmoid in that its maximum derivative is 1 (versus 0.25 for sigmoid), but it still saturates and can cause vanishing gradients in deep networks. The fact that tanh is zero-centered (outputs range from $-1$ to $+1$) makes it slightly easier to optimize than sigmoid in practice, since the gradients with respect to downstream weights are not all positive or all negative for a given input.

The ReLU function $\text{ReLU}(z) = \max(0, z)$ has derivative:

$$\text{ReLU}'(z) = \begin{cases} 1 & \text{if } z > 0 \\ 0 & \text{if } z \leq 0 \end{cases}$$

For positive inputs, ReLU passes the gradient through unchanged (gradient is 1). For negative inputs, ReLU blocks the gradient entirely (gradient is 0). This has two important consequences:

1. No vanishing gradient for active neurons: When a neuron is active (output > 0), the gradient passes through without attenuation. A 100-layer ReLU network can potentially propagate gradients all the way back to the first layer without shrinkage, as long as enough neurons remain active.
2. Dying ReLU problem: If a neuron consistently outputs negative values (typically due to a large negative bias or unfortunate weight initialization), it permanently stops learning because all gradients through it are zero. This can be detected during training by monitoring the fraction of neurons with zero activations, and mitigated with careful initialization or variants like Leaky ReLU.

The ReLU gradient's binary nature means it is not just non-vanishing for positive inputs: it also cannot explode. A gradient of 1 times 1 is still 1. This stability is a major practical advantage, and it explains why the shift from sigmoid to ReLU activations in the early 2010s dramatically improved the trainability of deep networks.

The GELU function (used in modern transformers) has a smooth, data-dependent gradient that avoids both sharp saturation and hard zeros. Its derivative involves the standard normal PDF and CDF, creating a smoother gradient landscape than ReLU. GELU's derivative is positive even for moderately negative inputs, which means neurons rarely die completely, though they do attenuate gradients in the negative regime.

The left plot confirms the textbook result: sigmoid gradients decay exponentially as they travel toward earlier layers, while ReLU gradients remain comparatively stable. The right plot shows how rapidly this ratio grows with depth for sigmoid, reaching many orders of magnitude even for moderate-depth networks.

Let's trace backpropagation through a complete two-layer network on a single training example to make the abstract formulas concrete. Walking through this computation by hand builds the kind of intuition that helps you debug real implementations.

We'll use a tiny network with:

• Input: $\mathbf{x} = [1, 2]^T$ (2 features)
• Layer 1: 3 hidden units with sigmoid activation
• Layer 2: 2 output units with softmax
• Target: class 0, so $\mathbf{y} = [1, 0]^T$

The weights and biases (chosen for illustration):

$$\mathbf{W}^{(1)} = \begin{pmatrix} 0.1 & 0.2 \\ 0.3 & 0.4 \\ 0.5 & 0.6 \end{pmatrix}, \quad \mathbf{b}^{(1)} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}$$ $$\mathbf{W}^{(2)} = \begin{pmatrix} 0.7 & 0.8 & 0.9 \\ 0.1 & 0.2 & 0.3 \end{pmatrix}, \quad \mathbf{b}^{(2)} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$$

Step 1, compute pre-activations in layer 1:

$$\mathbf{z}^{(1)} = \mathbf{W}^{(1)} \mathbf{x} = \begin{pmatrix} 0.1 \cdot 1 + 0.2 \cdot 2 \\ 0.3 \cdot 1 + 0.4 \cdot 2 \\ 0.5 \cdot 1 + 0.6 \cdot 2 \end{pmatrix} = \begin{pmatrix} 0.5 \\ 1.1 \\ 1.7 \end{pmatrix}$$

Step 2, apply sigmoid activation:

$$\mathbf{a}^{(1)} = \sigma(\mathbf{z}^{(1)}) = \begin{pmatrix} 0.622 \\ 0.750 \\ 0.846 \end{pmatrix}$$

Step 3, compute output logits:

$$\mathbf{z}^{(2)} = \mathbf{W}^{(2)} \mathbf{a}^{(1)} \approx \begin{pmatrix} 2.394 \\ 0.560 \end{pmatrix}$$

Step 4, apply softmax:

$$\hat{\mathbf{y}} = \text{softmax}(\mathbf{z}^{(2)}) \approx \begin{pmatrix} 0.853 \\ 0.147 \end{pmatrix}$$

Step 5, compute cross-entropy loss:

$$\mathcal{L} = -\log(\hat{y}_0) = -\log(0.853) \approx 0.159$$

This is a relatively small loss because the network already predicts the correct class (class 0) with high confidence (85.3%). In the backward pass, we'll see that this small loss produces correspondingly small gradients, meaning the network weights change very little for this example.

Step 1, output gradient. For softmax + cross-entropy:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} = \hat{\mathbf{y}} - \mathbf{y} = \begin{pmatrix} 0.853 - 1 \\ 0.147 - 0 \end{pmatrix} = \begin{pmatrix} -0.147 \\ 0.147 \end{pmatrix}$$

The network predicted class 0 with 85.3% confidence, and the target is class 0, so the gradient is small. If the prediction were wrong (say, 0.1 confidence on the correct class), the gradient would be $\hat{\mathbf{y}} - \mathbf{y} \approx [-0.9, 0.9]^T$, much larger. This is the self-correcting property of cross-entropy: bigger mistakes generate bigger gradients, which cause bigger parameter updates.

Step 2, gradients for layer 2 weights:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{W}^{(2)}} = \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} (\mathbf{a}^{(1)})^T = \begin{pmatrix} -0.147 \\ 0.147 \end{pmatrix} \begin{pmatrix} 0.622 & 0.750 & 0.846 \end{pmatrix}$$ $$= \begin{pmatrix} -0.091 & -0.110 & -0.124 \\ 0.091 & 0.110 & 0.124 \end{pmatrix}$$

Each entry $\frac{\partial \mathcal{L}}{\partial W^{(2)}_{ij}}$ is the product of the upstream gradient at output $i$ and the activation at hidden unit $j$. This outer product structure is fundamental: the weight gradient is proportional to how much that weight "contributed" (via the activation magnitude) multiplied by how much the output was "wrong" (the upstream gradient).

Step 3, propagate gradient back through layer 2:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{a}^{(1)}} = (\mathbf{W}^{(2)})^T \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(2)}} \approx \begin{pmatrix} -0.088 \\ -0.102 \\ -0.117 \end{pmatrix}$$

Step 4, gradient through sigmoid. We need $\sigma'(\mathbf{z}^{(1)}) = \mathbf{a}^{(1)} \odot (1 - \mathbf{a}^{(1)})$:

$$\sigma'(\mathbf{z}^{(1)}) = \begin{pmatrix} 0.622 \cdot 0.378 \\ 0.750 \cdot 0.250 \\ 0.846 \cdot 0.154 \end{pmatrix} = \begin{pmatrix} 0.235 \\ 0.188 \\ 0.130 \end{pmatrix}$$ $$\frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(1)}} = \frac{\partial \mathcal{L}}{\partial \mathbf{a}^{(1)}} \odot \sigma'(\mathbf{z}^{(1)}) \approx \begin{pmatrix} -0.021 \\ -0.019 \\ -0.015 \end{pmatrix}$$

Notice how the sigmoid derivative compressed the gradients further. The activations at this layer are relatively large (0.6 to 0.8), which puts them in the mildly-saturated region of the sigmoid, so the derivatives are between 0.13 and 0.24. The gradient coming out of the activation layer is roughly 3-4x smaller than what came in, which is the vanishing gradient effect in action even on this single-layer example.

Step 5, gradients for layer 1 weights:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{W}^{(1)}} = \frac{\partial \mathcal{L}}{\partial \mathbf{z}^{(1)}} \mathbf{x}^T$$

The worked example shows the key pattern: gradients start at the output, get multiplied by the transpose of each weight matrix to flow backward through layers, and get multiplied elementwise by activation derivatives at each nonlinearity. At each step, the gradient is the "responsibility signal" telling each parameter how much it contributed to the final error.

When a variable participates in multiple computations, its gradient accumulates. This is the gradient accumulation rule: the gradient of the loss with respect to any variable equals the sum of gradients from all downstream uses of that variable.

Consider a weight $w$ used in two computations: $z_1 = w \cdot x_1$ and $z_2 = w \cdot x_2$, with loss $\mathcal{L} = f(z_1, z_2)$. The total gradient is:

$$\frac{\partial \mathcal{L}}{\partial w} = \frac{\partial \mathcal{L}}{\partial z_1} \cdot x_1 + \frac{\partial \mathcal{L}}{\partial z_2} \cdot x_2$$

This addition arises directly from the multivariate chain rule: $w$ has two "paths" to the loss, and the total effect is the sum of both paths' contributions.

This principle applies beyond mini-batch training. It is also the key concept behind backpropagation through time (BPTT) in recurrent neural networks, where the same weight matrix is used at each time step. The gradient of the loss with respect to the shared weight is the sum of gradients from every time step where that weight was used. This is why BPTT in deep recurrent networks also suffers from vanishing gradients: the gradient from a token 50 time steps ago must flow through 50 weight matrices, and the product of their individual contributions tends to vanish. LSTMs and GRUs were specifically designed to mitigate this by adding gating mechanisms that create shorter gradient paths.

In practice, deep learning frameworks handle gradient accumulation automatically. PyTorch's .grad attribute accumulates gradients from all backward passes through a tensor, which is why optimizer.zero_grad() must be called before each training step. Failing to zero gradients leads to gradients from multiple batches summing together, a common source of bugs that produces unexpectedly fast apparent convergence followed by instability.

Understanding the computational cost of backpropagation helps explain why it's practical and what the tradeoffs are.

Forward pass cost: Each matrix multiplication $\mathbf{W} \in \mathbb{R}^{m \times n}$ operating on $\mathbf{x} \in \mathbb{R}^n$ costs $O(mn)$ operations. For a network with $L$ layers, the total forward pass cost is $O(\sum_l m_l n_l)$.

Backward pass cost: At each layer, we compute:

• Gradient with respect to weights: another $O(mn)$ operation (outer product)
• Gradient with respect to input: another $O(mn)$ operation (transposed multiply)

The backward pass costs roughly twice the forward pass. This is a fundamental result: computing gradients via backpropagation is only about 2x more expensive than the forward pass, regardless of the number of parameters. In practice, the backward pass is often 2-3x more expensive due to additional memory access patterns and accumulation, but the asymptotic complexity is the same.

Memory cost: The backward pass must store all intermediate activations from the forward pass, because it needs $\mathbf{a}^{(l)}$ to compute gradients at layer $l$. For a mini-batch of size $B$, this costs $O(B \cdot \sum_l m_l)$ memory. This memory requirement explains gradient checkpointing: one can trade computation for memory by not storing all activations, recomputing them as needed during the backward pass. The tradeoff is worth it when GPU memory is the bottleneck, which is common for large language models.

Without backpropagation: Computing gradients naively requires one forward pass per parameter (finite differences), costing $O(P)$ forward passes where $P$ is the number of parameters. For GPT-3's 175 billion parameters, this is completely infeasible. Backpropagation computes all gradients in $O(1)$ forward passes (roughly two), making gradient-based training of large neural networks possible.

The remarkable efficiency of backpropagation: It is worth pausing to appreciate just how extraordinary the 2x constant is. The gradient of a function with $P$ outputs with respect to $N$ inputs generally requires $O(NP)$ work to compute naively. Backpropagation achieves $O(N + P)$ work (essentially just one forward plus one backward pass) specifically because our "function" is a scalar loss and the chain rule lets us share computation across the $N$ outputs. This is the reverse-mode automatic differentiation perspective: it is efficient when there are many inputs and few outputs, which is exactly the case for neural network training where the loss is a single scalar.

Real implementations work with mini-batches rather than single examples. With a batch of $B$ examples, the input is $\mathbf{X} \in \mathbb{R}^{B \times n}$, and all forward and backward computations generalize naturally.

For a linear layer $\mathbf{Z} = \mathbf{X} \mathbf{W}^T + \mathbf{b}$:

• Forward: $\mathbf{Z} = \mathbf{X} \mathbf{W}^T$ (matrix multiplication, $\mathbf{Z} \in \mathbb{R}^{B \times m}$)
• Backward gradient w.r.t. weights: $\frac{\partial \mathcal{L}}{\partial \mathbf{W}} = \left(\frac{\partial \mathcal{L}}{\partial \mathbf{Z}}\right)^T \mathbf{X}$ (result shape $m \times n$)
• Backward gradient w.r.t. input: $\frac{\partial \mathcal{L}}{\partial \mathbf{X}} = \frac{\partial \mathcal{L}}{\partial \mathbf{Z}} \mathbf{W}$ (result shape $B \times n$)
• Backward gradient w.r.t. bias: $\frac{\partial \mathcal{L}}{\partial \mathbf{b}} = \sum_{i=1}^{B} \frac{\partial \mathcal{L}}{\partial \mathbf{z}_i}$ (sum over batch dimension)

The bias gradient requires summing over the batch dimension because the same bias vector is added to all examples. This is another example of gradient accumulation: the bias participates in $B$ computations (one per example), so its gradient is the sum over all $B$ paths.

Vectorized backpropagation improves performance and naturally averaged (or summed, depending on how the loss is normalized) over the batch. A batch gradient descent step moves the parameters in a direction that reduces the average loss across all examples in the batch. With large batches, this average is a better approximation of the true gradient over the full dataset, leading to more stable optimization. With small batches (stochastic gradient descent), the noisier gradient estimate can help escape sharp local minima, which is why SGD with a small batch size often generalizes better than full-batch gradient descent.

Let's implement backpropagation from scratch using only NumPy. We'll build a small framework with Value nodes that record their computational history and can propagate gradients backward.

How do we verify that our backpropagation implementation is correct? The standard technique is numerical gradient checking, which approximates gradients using finite differences and compares them to the analytical gradients from backprop.

For a small $\epsilon > 0$, the centered finite difference approximation of $\frac{\partial \mathcal{L}}{\partial w}$ is:

$$\frac{\partial \mathcal{L}}{\partial w} \approx \frac{\mathcal{L}(w + \epsilon) - \mathcal{L}(w - \epsilon)}{2\epsilon}$$

where:

• $\mathcal{L}(w + \epsilon)$: the loss when we perturb weight $w$ by $+\epsilon$
• $\mathcal{L}(w - \epsilon)$: the loss when we perturb weight $w$ by $-\epsilon$
• $2\epsilon$: the denominator giving the centered difference (more accurate than one-sided)

The centered difference approximates the true derivative with $O(\epsilon^2)$ error, compared to $O(\epsilon)$ for the one-sided difference. Using $\epsilon \approx 10^{-5}$ typically gives sufficient accuracy. Choosing $\epsilon$ too small causes numerical catastrophe: floating-point subtraction of nearly equal numbers loses precision, making the approximation noisy. Choosing $\epsilon$ too large increases the approximation error from the Taylor series truncation. The sweet spot is around $10^{-5}$ for float64 arithmetic.

To check whether the analytical and numerical gradients agree, compute the relative error:

$$\text{relative error} = \frac{|\mathbf{g}_{\text{analytical}} - \mathbf{g}_{\text{numerical}}|}{\max(|\mathbf{g}_{\text{analytical}}|, |\mathbf{g}_{\text{numerical}}|) + 10^{-8}}$$

where:

• $\mathbf{g}_{\text{analytical}}$: the gradient computed by backpropagation
• $\mathbf{g}_{\text{numerical}}$: the gradient approximated by finite differences
• The $10^{-8}$ in the denominator prevents division by zero for near-zero gradients

Guidelines for acceptable relative error:

• Below $10^{-5}$: excellent, backprop is almost certainly correct
• Between $10^{-5}$ and $10^{-3}$: possibly suspicious, check further
• Above $10^{-2}$: likely a bug in backpropagation

All relative errors are well below $10^{-4}$, confirming our backpropagation implementation is correct.

Gradient checking is slow: it requires two forward passes per parameter, making it $O(P)$ in the number of parameters. Use it for debugging, not production training. The standard workflow is:

1. Run gradient check on a small version of your model with a handful of inputs
2. Once the check passes, disable it and train normally
3. If training behaves unexpectedly, re-enable the check

It is important to do gradient checking at the beginning of a project, before you have spent significant time training. Finding a gradient bug after 10 hours of training is much more costly than finding it after 10 minutes. Experienced practitioners routinely run gradient checks as a sanity step before launching any serious training run.

There are a few situations where gradient checking gives misleading results. First, if your model uses random dropout during training, you need to disable dropout before running the check, since the numerical approximation and the analytical gradient will use different dropout masks. Second, if your loss function has flat regions (like a hinge loss away from the margin), the numerical gradient might be zero while the analytical gradient is nonzero, giving large relative error that is not actually a bug. Third, very small gradients can have large relative errors due to the denominator normalization, even when the absolute difference is tiny.

Even with a correct conceptual understanding, backprop implementations often have subtle bugs. Here are the most common ones:

Forgetting to zero gradients. In PyTorch, calling loss.backward() accumulates gradients into .grad. If you call backward twice without zeroing, the gradients from both passes sum together. Always call optimizer.zero_grad() before each backward pass. This is perhaps the most common beginner mistake in PyTorch, and it can produce training curves that initially look plausible (the loss decreases) but eventually diverge as the accumulated gradients from many mini-batches grow unboundedly large.

Wrong shapes in matrix multiplications. The gradient of the loss with respect to weights $\mathbf{W}$ has shape equal to $\mathbf{W}$. A common mistake is computing $\frac{\partial \mathcal{L}}{\partial \mathbf{z}} \cdot \mathbf{x}$ (a scalar times vector) instead of the correct outer product $\frac{\partial \mathcal{L}}{\partial \mathbf{z}} \mathbf{x}^T$ (which produces a matrix matching $\mathbf{W}$'s shape). NumPy's broadcasting can silently hide shape errors, returning something that looks plausible but has the wrong semantics.

Missing the activation derivative. When computing the gradient at pre-activation values $\mathbf{z}$, you must multiply by the derivative of the activation function. Missing this step gives gradients for the post-activation values $\mathbf{a}$ instead of $\mathbf{z}$, which is incorrect. This bug is easy to introduce when refactoring code: it is tempting to skip the activation backward when the activation is ReLU (since the derivative is simple), but the elementwise multiply by the ReLU mask must still happen.

Using the same array for input and output. In-place operations can corrupt the stored cache needed for the backward pass. Always use separate arrays for input and output in each layer. A related pitfall is mutating a parameter in the forward pass (e.g., normalizing it) before storing it in the cache, so the backward pass receives the modified value rather than the original.

Summing bias gradient incorrectly in batch mode. With a batch of examples, the bias gradient is the sum over examples (not the mean, not a single example's gradient). This is because each example contributes independently to the loss through the bias. If you use the mean instead of the sum, your bias updates will be underscaled by $1/B$, effectively giving the bias a much smaller effective learning rate than the weights.

Not accounting for broadcasting. When a smaller tensor is broadcast to match a larger one during the forward pass, the backward pass must sum over the broadcast dimensions. PyTorch handles this automatically, but manual implementations often forget it. For example, if you add a bias vector of shape $(m,)$ to an activation matrix of shape $(B, m)$, the forward pass broadcasts the bias across the batch dimension. The bias gradient in the backward pass must then sum over the $B$ dimension to produce a gradient of shape $(m,)$ matching the original bias.

Incorrect handling of numerical stability. Operations like log(softmax(x)) should be computed as log_softmax(x) rather than log(softmax(x)) to avoid the intermediate softmax values going to zero. If you compute softmax first and then take the log, tiny probabilities can underflow to zero, and log(0) produces negative infinity. PyTorch provides F.log_softmax and F.cross_entropy precisely to avoid these numerical issues.

PyTorch's autograd engine implements the same algorithm described above, but in a highly optimized and general form. Understanding how it works demystifies the .backward() call.

Every PyTorch tensor with requires_grad=True participates in building a dynamic computation graph. When you perform operations on such tensors, PyTorch records the operation and its inputs, building a directed acyclic graph (DAG) of Function objects.

Each Function object knows:

• What operation was performed
• What the inputs were (stored for the backward pass)
• How to compute the VJP (vector-Jacobian product) for that operation

The graph is dynamic: it is rebuilt from scratch on each forward pass. This allows Python control flow (if statements, loops) to affect the computation graph, making PyTorch's approach more flexible than TensorFlow's original static graph approach. In a static graph framework like early TensorFlow, you described the computation once and then ran it many times. The graph had to be fixed before execution began. PyTorch's approach of building the graph as you run the Python code means you can use standard Python debugging tools, inspect intermediate values, and change the network structure from one forward pass to the next.

The grad_fn attribute of each tensor stores the function that created it. When you inspect a tensor and see grad_fn=<AddmmBackward0>, that is PyTorch telling you this tensor was created by a matrix multiply and add operation, and it knows exactly how to propagate gradients back through that operation.

When you call loss.backward(), PyTorch:

1. Starts from the loss node with initial gradient $\frac{\partial \mathcal{L}}{\partial \mathcal{L}} = 1$
2. Traverses the computation graph in reverse topological order
3. For each Function node, calls its backward() method with the upstream gradient
4. The backward() method computes VJPs and passes them to inputs
5. Gradients accumulate in each leaf tensor's .grad attribute

The gradients match to floating-point precision, confirming that our manual NumPy implementation computes the same gradients as PyTorch's optimized autograd engine.

PyTorch's gradient accumulation behavior enables several useful training patterns:

Gradient accumulation for large batches. When GPU memory is too small for the desired batch size, you can split a large batch into smaller micro-batches, call backward() on each, and only call optimizer.step() after all micro-batches. The gradients accumulate across micro-batches, producing the same gradient update as if you had run a single large batch. This technique is widely used in large language model training, where batch sizes of thousands or millions of tokens are desirable but a single token sequence might be tens of thousands of characters:

The gradient norms are essentially identical, confirming that gradient accumulation faithfully simulates a larger batch size. We'll explore gradient accumulation further when we discuss distributed training in later chapters.

The vanishing gradient problem shaped much of neural network architecture design over the last decade. Understanding the problem from first principles helps explain why modern architectures look the way they do.

The most influential solution to vanishing gradients is the residual connection, introduced by He et al. (2015) in the ResNet paper. Instead of computing $\mathbf{a}^{(l)} = f(\mathbf{z}^{(l)})$, a residual block computes:

$$\mathbf{a}^{(l)} = f(\mathbf{z}^{(l)}) + \mathbf{a}^{(l-1)}$$

The addition of the input $\mathbf{a}^{(l-1)}$ creates a shortcut path that bypasses the nonlinear transformation. From the backpropagation perspective, the gradient through an addition node is simply the identity: the gradient flows through unchanged. This means there is always a direct gradient highway from the loss all the way back to the earliest layers, regardless of how many layers intervene. The network can learn to use the residual path as heavily as needed, and in the early stages of training when the weights are small, the residual connections dominate and keep gradients flowing.

Transformer models, which power modern language models, use residual connections around both the self-attention sublayer and the feed-forward sublayer in every block. This is not coincidence: without them, even relatively shallow transformers (12-24 layers) would suffer from vanishing gradients and be very difficult to train.

Weight initialization also affects gradient flow. If weights are initialized too small, activations shrink layer by layer in the forward pass, and the gradients are correspondingly small in the backward pass. If weights are too large, activations grow and saturate, again killing gradients.

Xavier (Glorot) initialization sets the variance of weights to $\frac{2}{n_{\text{in}} + n_{\text{out}}}$, which keeps the variance of activations and gradients approximately constant across layers for tanh networks. He initialization sets variance to $\frac{2}{n_{\text{in}}}$, better suited for ReLU networks because ReLU cuts the variance roughly in half (by zeroing negative activations). Getting initialization right is often the difference between a network that starts training smoothly and one that immediately diverges or stalls.

Batch normalization (Ioffe and Szegedy, 2015) normalizes the pre-activation values at each layer to have zero mean and unit variance over the batch. By keeping activations in a reasonable range, it prevents both saturation (which kills gradients) and explosion (which causes instability). It also provides a direct regularization effect and allows higher learning rates, substantially speeding up training. The gradient through a batch normalization layer is more complex than through a plain activation (it involves the batch statistics), but it has the property of preventing activations from drifting into the extreme values where activation derivatives vanish.

Layer normalization, the variant used in transformers, applies the same idea but normalizes across the feature dimension rather than the batch dimension, making it better suited for variable-length sequences and training with batch size 1.

The key implementation choices when working with backpropagation are:

• Learning rate ($\eta$): Controls how large each gradient step is. Too large causes divergence; too small causes slow convergence. We'll cover learning rate selection in detail in the Stochastic Gradient Descent chapter.
• Numerical precision: Float32 is the standard for training (good speed/memory tradeoff); float64 is used for gradient checking (higher precision). Mixed precision training (float16 for forward, float32 for gradients) is common in large model training.
• Gradient clipping threshold: Prevents exploding gradients in deep or recurrent networks by capping the gradient norm. Covered in the Gradient Clipping chapter.
• Gradient checkpoint interval: The memory vs. compute tradeoff for very deep networks. Checkpointing at every $k$ layers reduces memory by a factor of $k$ at the cost of recomputing $k-1$ forward passes.

Backpropagation has some practical limitations that shape how neural networks are designed and trained.

The vanishing gradient problem remains the most persistent challenge. When gradients are multiplied across many sigmoid or tanh layers, they can shrink to numerically zero before reaching the earliest layers. The introduction of ReLU activations, residual connections (which add skip paths that provide direct gradient highways), batch normalization, and careful weight initialization all exist primarily to address this problem. As we'll see in subsequent chapters, architectural choices like residual connections in ResNets and layer normalization in Transformers are fundamentally about keeping gradient flow healthy.

Backpropagation requires storing all intermediate activations for the backward pass. For large models and large batches, this memory requirement can exceed GPU capacity. Gradient checkpointing (also called activation recomputation) trades computation for memory by not storing all intermediate values and recomputing them during the backward pass. Modern large language models routinely use gradient checkpointing during training. The memory savings can be substantial: for a 100-layer transformer with long sequences, checkpointing can reduce activation memory by 10-15x at the cost of roughly 30-40% additional computation time.

Backpropagation computes exact gradients with respect to the loss, which is both a strength and a limitation. The gradients are correct for the exact computation performed, but they contain no information about the loss landscape beyond the current point. Second-order methods like Newton's method use the Hessian (matrix of second derivatives) to better inform the update direction, but computing and storing the full Hessian for even moderately sized networks is infeasible. Methods like K-FAC and Shampoo approximate the Hessian in structured ways to get some second-order benefits without full Hessian computation, but they remain significantly more expensive than first-order methods like Adam.

Non-differentiable operations break backpropagation's clean mathematical story. Operations like argmax, sorting, or branching on data values have zero gradient almost everywhere. Modern deep learning works around this through soft approximations (softmax instead of argmax), straight-through estimators (for binary operations), and reinforcement learning signals (for discrete decisions). The straight-through estimator is particularly elegant: during the forward pass, you apply a non-differentiable operation; during the backward pass, you simply pretend the operation was the identity function and pass gradients through unchanged. This is mathematically unjustified but empirically works well for many applications including training binary neural networks and quantization-aware training.

Biological plausibility is another limitation, though it matters more for neuroscience than for engineering. Backpropagation requires exact knowledge of the forward weights to compute the backward pass (the $\mathbf{W}^T$ in the gradient formula). Neurons in biological brains do not have access to the exact synaptic weights of the neurons they synapse onto, making backpropagation unlikely to be the mechanism by which biological networks learn. Alternative learning rules like forward-forward, target propagation, and perturbation-based methods are active research areas trying to achieve the efficiency of backpropagation while respecting biological constraints. For practical deep learning, though, biological plausibility is not a concern.

Despite these limitations, backpropagation fundamentally changed what was computationally possible. Before efficient gradient computation, training more than two or three layer networks was practically infeasible. The algorithm enabled CNNs for image recognition, RNNs for sequence modeling, Transformers for language AI, and the billion-parameter models that power today's AI systems. In the next chapter, we'll use these gradients with stochastic gradient descent and its variants to train networks efficiently.

Backpropagation is the algorithm that makes training neural networks with gradient descent tractable:

• Computational graphs decompose neural network computations into simple operations, each with a known local gradient. The graph structure makes the chain rule applicable layer by layer.
• The chain rule shows how to compute the gradient of a composite function by multiplying local gradients from output to input. The multiplicative structure is both the source of backprop's power (gradients compose cleanly) and its vulnerability (vanishing and exploding gradients).
• The backward pass traverses the computational graph in reverse, computing gradients for all parameters in time proportional to a single forward pass, roughly 2x the cost of the forward pass.
• Jacobian matrices formalize gradient flow through vector operations. In practice, backprop never forms full Jacobians but computes vector-Jacobian products directly, which is much more efficient.
• Gradient flow through activations depends on the activation: sigmoid and tanh saturate (vanishing gradients), while ReLU passes gradients through unchanged for positive inputs. Modern architectures use ReLU or GELU activations and add residual connections to keep gradients flowing through many layers.
• Numerical gradient checking verifies implementations by comparing analytical gradients against finite-difference approximations. Relative errors below $10^{-5}$ indicate a correct implementation.
• Common bugs include forgetting to zero gradients, wrong shapes in matrix products, missing activation derivatives, incorrect bias gradient summation over batches, and numerical instability from computing log-softmax naively.
• PyTorch autograd implements the same VJP-based algorithm, building a dynamic computation graph as code runs and traversing it in reverse on the .backward() call. The dynamic graph enables full Python control flow in models.
• Gradient accumulation allows simulating large batch sizes by accumulating gradients across multiple smaller passes before calling optimizer.step().
• Architectural solutions to vanishing gradients include residual connections (which create direct gradient highways), careful weight initialization (Xavier for tanh, He for ReLU), and normalization layers (batch normalization, layer normalization) that keep activations in healthy ranges.