Backpropagation: The Algorithm That Makes Deep Learning Possible

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 arguably 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 review the chain rule that makes gradient computation tractable. You'll trace through forward and backward passes step by step, understand gradient accumulation, and implement backpropagation from scratch. By the end, you'll see exactly how modern deep learning frameworks like PyTorch 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 of the gradient tells us how sensitive the function is to changes in the corresponding parameter. To minimize the loss, we move in the opposite direction, the direction of steepest decrease:

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 (typically 0.001 to 0.1) that controls how large each step is
• $\nabla_\mathbf{w} \mathcal{L}$: the gradient of the loss with respect to parameters, a vector with the same shape as $\mathbf{w}$, where each entry $\frac{\partial \mathcal{L}}{\partial w_i}$ indicates how much the loss would increase if we increased $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's activations, which depend on the second layer's activations, which depend on the first layer's activations, which depend on the input and the first layer's weights. This chain of dependencies seems hopelessly complex to differentiate.

Backpropagation solves this problem elegantly. By applying the chain rule in reverse order, from output back to input, we can compute all gradients in time proportional to a single forward pass. This efficiency is what makes training deep networks 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. A computational graph makes these explicit. Consider the simple expression:

This decomposes into two operations:

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

where:

• $x, y, z$: the input variables to the function
• $a$: an intermediate value that stores the result of the first operation
• $f$: the final output of the computation

The graph shows data flowing from inputs (x, y, z) through operations (+, ×) to the output (f). The intermediate value $a = x + y$ is stored at the addition node. This decomposition is the key insight: complex expressions are just chains of simple operations.

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

• $\mathbf{x}$: the input vector (or the output from the previous layer)
• $\mathbf{W}$: the weight matrix, with shape (output neurons × input neurons)
• $\mathbf{b}$: the bias vector, with one value per output neuron
• $\sigma$: the activation function (e.g., ReLU, sigmoid)
• $\mathbf{a}$: the output activations

Let's break this into its constituent operations:

1. $\mathbf{z}_1 = \mathbf{W}\mathbf{x}$ (matrix multiplication: transforms input through learned weights)
2. $\mathbf{z}_2 = \mathbf{z}_1 + \mathbf{b}$ (bias addition: shifts the output)
3. $\mathbf{a} = \sigma(\mathbf{z}_2)$ (activation: introduces nonlinearity)

Each node in the graph represents an operation that we know how to differentiate. The power of computational graphs is that they make the structure of the computation explicit, which is exactly what we need for backpropagation.

Computational graphs provide two crucial benefits for gradient computation:

1. Modularity: Each node has a well-defined local gradient. We don't need to derive the gradient of the entire network at once; we just need the gradient of each primitive operation.

2. Efficiency: By storing intermediate values during the forward pass, we can reuse them during the backward pass. This avoids redundant computation.

The forward pass computes values flowing left to right. The backward pass computes gradients flowing right to left. This symmetry is the essence of backpropagation.

We've seen how computational graphs break complex neural network computations into simple, modular operations. Now we face the central question: how do we compute the gradient of the loss with respect to every parameter in this graph?

The answer lies in the chain rule from calculus. To truly understand backpropagation, we need to see the chain rule not as a formula to memorize, but as a way of thinking about how influence propagates through a system of nested computations. Every operation in a neural network transforms its inputs in some way, and the chain rule tells us exactly how to trace the effect of any input change through all those transformations to the final output.

Before diving into formulas, let's build intuition with a physical analogy. Imagine a factory assembly line where raw materials pass through several processing stations, each transforming the product in some way, until a final item emerges. If you want to know how changing the raw material affects the final product, you need to trace that change through every station in sequence.

The same logic applies to neural networks. An input $x$ passes through layer after layer of transformations until it produces a loss value $\mathcal{L}$. To understand how $x$ affects $\mathcal{L}$, we trace the influence through each transformation, and the chain rule gives us the precise mathematical tool to do this.

Let's start with the simplest case: two functions arranged in a pipeline. First, $g$ transforms input $x$ into some intermediate value, then $f$ transforms that intermediate value into the final output. We write this composition as $f(g(x))$.

The question we want to answer is: if we nudge $x$ by a tiny amount, how much does the final output change?

Think of it as a two-stage amplification process:

1. Stage 1: The nudge to $x$ causes $g(x)$ to change. If $g$ is very sensitive to $x$ (steep slope), a small change in $x$ produces a large change in $g$. If $g$ is insensitive (flat slope), the change is small.

2. Stage 2: This change in $g(x)$ then causes $f$ to change. Again, the magnitude depends on how sensitive $f$ is to its input.

The total effect is the product of these two sensitivities, as each stage amplifies (or attenuates) the signal from the previous stage:

where:

• $\frac{df}{dx}$: the total derivative, how much the final output $f$ changes per unit change in the original input $x$. This is what we ultimately want to know.
• $\frac{df}{dg}$: the local gradient of the outer function, how sensitive $f$ is to changes in its immediate input $g$. We can compute this knowing only the function $f$, without caring about where $g$ came from.
• $\frac{dg}{dx}$: the local gradient of the inner function, how sensitive $g$ is to changes in $x$. We can compute this knowing only the function $g$.

Why multiplication? If $g$ doubles any change in $x$ (i.e., $\frac{dg}{dx} = 2$), and $f$ triples any change in $g$ (i.e., $\frac{df}{dg} = 3$), then a unit change in $x$ ultimately causes a $2 \times 3 = 6$ unit change in $f$. The effects compound multiplicatively.

This modularity is the key insight: we can break a complex derivative into a product of simple, local derivatives. Each function in the chain only needs to know how to differentiate itself with respect to its immediate inputs. It doesn't need to understand the entire pipeline.

Let's make this concrete with an example. Consider the composite function $f(x) = \sin(x^2)$, which we can decompose as:

• Inner function: $g(x) = x^2$
• Outer function: $f(g) = \sin(g)$

To find $\frac{df}{dx}$, we compute each local gradient and multiply:

1. $\frac{dg}{dx} = 2x$ (derivative of $x^2$)
2. $\frac{df}{dg} = \cos(g)$ (derivative of $\sin(g)$)
3. $\frac{df}{dx} = \cos(g) \cdot 2x = \cos(x^2) \cdot 2x = 2x\cos(x^2)$

The analytical gradient computed via the chain rule matches the numerical approximation to within $10^{-9}$, well within floating-point precision. This verification technique (comparing analytical gradients to numerical approximations) will become crucial when we implement backpropagation and need to verify our code is correct.

The single-variable chain rule handles a simple pipeline where one thing leads to another in a straight line. But neural networks are more complex: a single variable can influence the output through multiple intermediate pathways. This is where the multivariate chain rule becomes essential.

The Problem: Multiple Pathways of InfluenceLink Copied

Consider a concrete scenario: in a neural network, a single weight $w$ might affect multiple neurons in the next layer, each of which contributes to the final loss. How do we account for all these influences?

Think of it like a river that splits into multiple tributaries before reaching the ocean. If you add a drop of dye at the source, it flows through all the tributaries. To measure the total amount of dye reaching the ocean, you must sum the contributions from each path.

The same principle applies to gradients: when a variable affects the output through multiple paths, we sum the contributions from each path. If changing $w$ increases the loss through one path and decreases it through another, the net effect is the sum of these opposing influences.

The FormulaLink Copied

Mathematically, if a variable $x$ affects the loss $L$ through several intermediate variables $y_1, y_2, \ldots, y_n$, the total gradient is:

Let's unpack each component:

• $L$: the final loss, a single scalar value that we're trying to minimize
• $x$: the variable we want to understand the influence of (e.g., a weight or bias)
• $y_1, y_2, \ldots, y_n$: all the intermediate values that directly depend on $x$
• $\frac{\partial L}{\partial y_i}$: how sensitive the loss is to changes in $y_i$ (this is computed by backpropagating from the output)
• $\frac{\partial y_i}{\partial x}$: how sensitive $y_i$ is to changes in $x$ (this is computed locally)

Each term $\frac{\partial L}{\partial y_i} \cdot \frac{\partial y_i}{\partial x}$ represents one "pathway of influence," the effect of $x$ on $L$ through the specific route via $y_i$. The sum aggregates all these pathways.

Why This Matters for Neural NetworksLink Copied

This summation principle has profound implications for neural network training:

1. Gradient accumulation: When the same weight is used multiple times (as in recurrent networks processing sequences, or in architectures with shared weights), its gradient naturally accumulates contributions from each use.

2. Fan-out patterns: When a layer's output feeds into multiple downstream computations, the gradient flowing back must aggregate all these influences.

3. Batch processing: When we process multiple examples simultaneously, each example contributes to the gradient, and we sum (or average) these contributions.

Now we can see how the chain rule maps onto computational graphs, and this mapping is what transforms the abstract mathematics into a practical algorithm.

From Math to AlgorithmLink Copied

Recall that a computational graph represents a neural network's computation as a directed graph: nodes are operations, and edges show how data flows from inputs to outputs. During the forward pass, we compute values flowing from left to right. During the backward pass, we compute gradients flowing from right to left, from the loss back toward the inputs.

The key insight is that each node only needs to know its local gradient. A node doesn't need to understand the entire network; it just needs to know two things:

1. Forward: How to compute its output from its inputs
2. Backward: How to compute the gradient of its output with respect to its inputs

This locality is what makes backpropagation modular and scalable.

The Backward Pass FormulaLink Copied

For any node $v$ in the graph, the gradient of the loss with respect to $v$ follows directly from the multivariate chain rule:

where:

• $L$: the loss function (a scalar value we're trying to minimize)
• $v$: the current node whose gradient we're computing
• $\text{children}(v)$: the set of all nodes $u$ that directly receive $v$ as input (i.e., nodes downstream from $v$)
• $\frac{\partial L}{\partial u}$: the gradient of the loss with respect to each child node (already computed during backpropagation)
• $\frac{\partial u}{\partial v}$: the local gradient, how much node $u$'s output changes when $v$ changes

The Algorithm in Plain LanguageLink Copied

This formula encodes a simple, elegant algorithm. For each node $v$, working backward from the loss:

1. Look downstream: Find all nodes $u$ that use $v$ as an input (these are $v$'s "children" in graph terminology)

2. Collect incoming gradients: Each child $u$ sends back a gradient signal $\frac{\partial L}{\partial u}$, which tells us how sensitive the loss is to $u$'s output

3. Scale by local sensitivity: Multiply each incoming gradient by $\frac{\partial u}{\partial v}$, the local gradient that describes how much $u$'s output changes when $v$ changes

4. Sum the contributions: Add up all these scaled gradients to get the total gradient at $v$

Why This Works: The Recursive StructureLink Copied

The beauty of this approach is its recursive structure. We start at the loss node, where $\frac{\partial L}{\partial L} = 1$ (the loss is perfectly sensitive to itself). Then we work backward through the graph, layer by layer:

• Local information is sufficient: Each node computes its local gradient using only the operation type and values stored during the forward pass
• Gradients flow backward: Each node receives gradient signals from its children (which have already been computed) and passes gradient signals to its parents
• No global coordination needed: The algorithm is inherently parallel, and nodes at the same depth can compute their gradients simultaneously

This is the essence of backpropagation: a systematic, efficient way to apply the chain rule to computational graphs.

When a variable affects the loss through multiple paths, we sum the gradients from each path. This is why the chain rule involves a sum in the multivariate case.

With the chain rule understood, we're ready to trace through an actual neural network. The forward pass propagates input data through the network, computing and storing intermediate values at each node. These stored values are essential for the backward pass; without them, we couldn't compute the local gradients.

Think of the forward pass as leaving breadcrumbs: as we compute each intermediate value, we save it for later use when gradients flow backward.

Let's trace through a concrete 2-layer network with the following architecture:

• Input: $\mathbf{x} \in \mathbb{R}^2$ (2 features)
• Hidden layer: 3 neurons with ReLU activation
• Output layer: 1 neuron with sigmoid activation (for binary classification)

This small network is simple enough to trace by hand, yet complex enough to illustrate all the key concepts.

Now let's perform the forward pass, storing all intermediate values:

Notice how the ReLU activation zeroed out the negative pre-activation value in position 2 (z1[2] = -0.1 became a1[2] = 0). This sparsity is characteristic of ReLU networks, and it will have important implications for the backward pass, as we'll see shortly.

With the forward pass complete, we've arrived at a prediction. But how good is it? We need a way to quantify the error, and this is the role of the loss function.

For binary classification, we use binary cross-entropy (BCE), which measures how well the predicted probability matches the true label:

where:

• $y$: the true label, either 0 or 1
• $\hat{y}$: the predicted probability from the sigmoid output, a value between 0 and 1
• $\log$: the natural logarithm
• $\mathcal{L}$: the loss value (always non-negative)

Understanding the FormulaLink Copied

This formula has an elegant interpretation that becomes clear when we consider the two cases:

When $y = 1$ (positive class): The second term vanishes (since $1-y = 0$), leaving $\mathcal{L} = -\log(\hat{y})$.

• If $\hat{y} \approx 1$ (correct prediction): $\log(1) = 0$, so loss is small
• If $\hat{y} \approx 0$ (wrong prediction): $\log(0) \rightarrow -\infty$, so loss is large

When $y = 0$ (negative class): The first term vanishes (since $y = 0$), leaving $\mathcal{L} = -\log(1-\hat{y})$.

• If $\hat{y} \approx 0$ (correct prediction): $\log(1) = 0$, so loss is small
• If $\hat{y} \approx 1$ (wrong prediction): $\log(0) \rightarrow -\infty$, so loss is large

The logarithm creates an asymmetric penalty: confident wrong predictions are punished severely, while confident correct predictions are rewarded with near-zero loss.

The loss of approximately 0.59 is moderate, not terrible, but clearly room for improvement. Our prediction of ~0.56 is below the target of 1.0, meaning the network is underconfident about the positive class. The loss quantifies this gap: a perfect prediction of 1.0 would yield a loss near 0, while our current prediction produces a loss that will drive gradient updates to increase the output probability. Now we need to compute gradients to reduce this loss.

Now comes the heart of backpropagation. We've completed the forward pass, storing all intermediate values. We've computed the loss. Now we need to answer the fundamental question: how should we adjust each weight to reduce the loss?

The backward pass answers this by computing $\frac{\partial \mathcal{L}}{\partial w}$ for every weight $w$ in the network. Armed with these gradients, we can update weights in the direction that decreases the loss.

The backward pass is like unwinding a ball of yarn. We start at the end (the loss) and work our way back to the beginning (the inputs), carefully tracking how each intermediate value influenced the final result. At each step, we apply the chain rule to decompose the complex derivative into simpler, local derivatives.

Every backward pass begins at the loss function. The loss is our starting point because it's what we're trying to minimize. It's the "error signal" that will propagate backward through the network.

Before we can backpropagate through the network layers, we need to know: how does the loss change when the prediction changes? This is the gradient $\frac{\partial \mathcal{L}}{\partial \hat{y}}$, the sensitivity of the loss to our prediction.

Let's derive this step by step, starting from the BCE formula:

Step 1: Write out the loss function:

This formula has two terms, but only one is ever "active" for a given example. When $y = 1$ (positive class), the second term vanishes because $(1-y) = 0$. When $y = 0$ (negative class), the first term vanishes because $y = 0$.

Step 2: Differentiate each term with respect to $\hat{y}$:

For the first term, we use the fact that $\frac{d}{dx}\log(x) = \frac{1}{x}$:

For the second term, we need the chain rule because we're differentiating $\log(1-\hat{y})$:

Step 3: Combine the terms:

Interpreting the GradientLink Copied

This formula reveals something important about how the loss "pushes" the prediction:

The loss function naturally provides stronger correction signals for more confident wrong predictions. This is exactly what we want: big mistakes get big corrections.

Now we need to push the gradient one step further back, through the sigmoid activation function. This is our first example of backpropagating through a nonlinearity, and it reveals a beautiful mathematical property that makes the combined BCE + sigmoid gradient remarkably simple.

The sigmoid function $\sigma(z) = \frac{1}{1 + e^{-z}}$ squashes any input into the range (0, 1). But what's its derivative? How sensitive is the sigmoid output to changes in its input?

Deriving the sigmoid derivative:

Step 1: Rewrite sigmoid using exponent notation for easier differentiation:

Step 2: Apply the chain rule. Let $u = 1 + e^{-z}$, so $\sigma = u^{-1}$:

where:

• $u = 1 + e^{-z}$: a substitution variable to simplify the differentiation
• $\frac{\partial \sigma}{\partial u} = -u^{-2}$: the derivative of $u^{-1}$ with respect to $u$
• $\frac{\partial u}{\partial z} = -e^{-z}$: the derivative of $1 + e^{-z}$ with respect to $z$ (the constant 1 vanishes, and $\frac{d}{dz}e^{-z} = -e^{-z}$)

The two negative signs cancel, giving us a positive result.

Step 3: Here's where it gets elegant. We can rewrite this messy expression entirely in terms of $\sigma(z)$ itself:

• Note that $\sigma(z) = \frac{1}{1 + e^{-z}}$
• And $1 - \sigma(z) = 1 - \frac{1}{1 + e^{-z}} = \frac{e^{-z}}{1 + e^{-z}}$
• Therefore: $\sigma(z) \cdot (1 - \sigma(z)) = \frac{1}{1 + e^{-z}} \cdot \frac{e^{-z}}{1 + e^{-z}} = \frac{e^{-z}}{(1 + e^{-z})^2}$

This gives us the remarkably simple result:

where:

• $\sigma(z)$: the sigmoid output, a value between 0 and 1
• $1 - \sigma(z)$: the complement of the sigmoid output
• $\sigma(z)(1 - \sigma(z))$: the product, which is maximized when $\sigma(z) = 0.5$

The Geometry of the Sigmoid DerivativeLink Copied

This formula has a beautiful geometric interpretation:

• Maximum at the midpoint: The derivative peaks at $\sigma(z) = 0.5$, where $\sigma'(z) = 0.5 \times 0.5 = 0.25$. This is where the sigmoid is steepest.
• Vanishes at the extremes: As $\sigma(z) \rightarrow 0$ or $\sigma(z) \rightarrow 1$, the derivative approaches zero. The sigmoid flattens out.

This behavior has profound implications for training: when activations saturate near 0 or 1, gradients become tiny, making learning slow or impossible. This is the vanishing gradient problem for sigmoid activations.

The visualization shows why sigmoid activations are problematic in deep networks. When the pre-activation $z$ is far from zero (either very positive or very negative), the derivative approaches zero, causing gradients to vanish during backpropagation. Only inputs near $z = 0$ produce meaningful gradients.

Now comes the mathematical payoff that makes backpropagation through classification networks elegant. When we chain the BCE gradient with the sigmoid gradient using the chain rule, something remarkable happens:

where:

• $z_2$: the pre-activation value (logit) before the sigmoid, i.e., the raw output of the linear layer
• $\hat{y} = \sigma(z_2)$: the sigmoid output (predicted probability)
• $\frac{\partial \mathcal{L}}{\partial \hat{y}}$: the gradient of the loss with respect to the prediction (computed in the previous section)
• $\frac{\partial \hat{y}}{\partial z_2}$: the sigmoid derivative (computed above)

Substituting our expressions for $\frac{\partial \mathcal{L}}{\partial \hat{y}} = -\frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}}$ and $\frac{\partial \hat{y}}{\partial z_2} = \hat{y}(1-\hat{y})$:

Let's expand this carefully:

The complex-looking derivatives collapse into something beautifully simple:

where:

• $\hat{y}$: the predicted probability (sigmoid output)
• $y$: the true label (0 or 1)
• $\hat{y} - y$: the prediction error, which is positive when we over-predict and negative when we under-predict

This is remarkably simple: the gradient at the pre-sigmoid layer is just the prediction minus the target, the prediction error!

This elegant simplification is not a coincidence. Sigmoid and cross-entropy were designed to work together, and the mathematical complexity cancels out perfectly. The result is a gradient that's both easy to compute and deeply intuitive:

The gradient magnitude is proportional to how wrong we are. A prediction of 0.9 when the target is 0 produces gradient 0.9, while a prediction of 0.6 produces gradient 0.6. Bigger mistakes get bigger corrections.

Now we can implement the complete backward pass. The code below traces gradients from the loss back through both layers, computing gradients for all weights and biases.

Let's trace through what just happened, connecting each computation to the chain rule:

1. Starting at the loss: We computed $dz_2 = \hat{y} - y$, the gradient at the output layer. This single number tells us how wrong our prediction was, and it's the "error signal" that will propagate backward.

2. Gradient for output weights ($dW_2$): We asked: "How does each weight in $W_2$ affect $z_2$?" Since $z_2 = W_2 \cdot a_1 + b_2$, each weight $W_{2,ij}$ is multiplied by the corresponding activation $a_{1,j}$. By the chain rule:

3. Gradient for hidden activations ($da_1$): To continue backpropagating, we need to know how $a_1$ affected the loss. Since $z_2 = W_2 \cdot a_1$, changing $a_1$ by a small amount changes $z_2$ by $W_2$ times that amount:

4. Through the ReLU ($dz_1$): The ReLU function $\text{ReLU}(z) = \max(0, z)$ has a simple derivative:

This acts as a gate: gradients pass through unchanged for positive pre-activations, but are completely blocked for negative ones. The gradient becomes $dz_1 = da_1 \odot \mathbf{1}_{z_1 > 0}$, where $\odot$ denotes element-wise multiplication.

5. Gradient for hidden weights ($dW_1$): Following the same pattern as the output layer:

The ReLU Gate and "Dying Neurons"Link Copied

Notice how the ReLU gradient acts as a gate. Where $z_1 > 0$, the gradient passes through unchanged. Where $z_1 \leq 0$, the gradient is blocked (zeroed). In our example, the third neuron had $z_1[2] = -0.1 < 0$, so its gradient was blocked entirely.

This is the "dying ReLU" phenomenon: neurons that never activate don't receive gradient updates. If a neuron's pre-activation is always negative (perhaps due to poor initialization or a large negative bias), it will never learn. This is why variants like Leaky ReLU (which allow small gradients for negative inputs) are sometimes preferred.

Each gradient has a specific shape that matches its corresponding parameter:

The gradient shapes always match the parameter shapes. This must be true because we subtract gradients from parameters during the update step, and subtraction requires matching shapes.

We've derived and implemented the backward pass analytically. But how do we know we got it right? A single sign error or transposed matrix can silently produce wrong gradients, leading to a network that trains poorly or not at all.

The solution is gradient checking: verify analytical gradients against numerical approximations. The idea is simple: perturb each parameter slightly and measure how the loss changes. If our analytical gradient is correct, it should match this numerical estimate.

The Centered Difference FormulaLink Copied

The centered difference formula gives an accurate numerical approximation:

where:

• $w$: a single scalar parameter we're checking
• $\epsilon$: a small perturbation, typically $10^{-7}$
• $\mathcal{L}(w + \epsilon)$: the loss computed with $w$ increased by $\epsilon$
• $\mathcal{L}(w - \epsilon)$: the loss computed with $w$ decreased by $\epsilon$

Why centered difference? The one-sided formula $\frac{\mathcal{L}(w + \epsilon) - \mathcal{L}(w)}{\epsilon}$ has $O(\epsilon)$ error. The centered difference cancels first-order error terms, giving $O(\epsilon^2)$ accuracy, typically 5-7 orders of magnitude better for the same $\epsilon$.

The analytical and numerical gradients match to within $10^{-10}$, well within the expected precision for $\epsilon = 10^{-7}$. This confirms our backpropagation implementation is correct.

When a variable is used multiple times in a computation, its gradient accumulates contributions from all uses. This is the multivariate chain rule in action.

Consider a weight matrix $\mathbf{W}$ used to process a batch of examples. Each example contributes to the total loss, and each contribution produces its own gradient for $\mathbf{W}$. Since the total loss is the average (or sum) of individual losses, the total gradient is the sum of individual gradients:

where:

• $m$: the batch size (number of examples processed together)
• $\mathcal{L}_i$: the loss contribution from example $i$
• $\frac{\partial \mathcal{L}_i}{\partial \mathbf{W}}$: the gradient of the loss with respect to $\mathbf{W}$ from example $i$ alone, a matrix with the same shape as $\mathbf{W}$
• $\frac{\partial \mathcal{L}}{\partial \mathbf{W}}$: the total gradient, summed across all examples

This summation happens automatically in matrix operations when processing batches. When we compute $dW = \frac{1}{m} \cdot dZ \cdot A^T$, the matrix multiplication implicitly sums over the batch dimension, and the $\frac{1}{m}$ factor converts the sum to an average.

Let's see how gradient accumulation works with a batch of examples:

The predictions show mixed accuracy on this random batch, with some examples classified correctly while others are not. This is expected since we're using the same small network trained on XOR, which doesn't generalize to random inputs. The key observation is the gradient shape: dW2 has shape (1, 3), matching W2 exactly. The matrix multiplication dz2 @ a1.T (shapes 1×4 and 4×3) automatically sums over the batch dimension, and the $\frac{1}{m}$ factor converts this sum to an average, making the effective learning rate independent of batch size.

Let's visualize how gradients flow backward through our network:

The gradient magnitudes can vary significantly across layers. In deep networks, this can lead to vanishing gradients (gradients become tiny in early layers) or exploding gradients (gradients become huge). Techniques like batch normalization, residual connections, and careful initialization help mitigate these issues.

Let's train a deeper network and monitor how gradient magnitudes evolve across layers and epochs:

The gradient magnitude plot reveals a key pattern: gradients tend to be smaller in earlier layers (closer to the input). This is the vanishing gradient effect in action. Each layer's backward pass multiplies gradients by weight matrices and activation derivatives, which can compound to shrink gradients exponentially. In this 4-layer network, the effect is modest, but in networks with 50+ layers, it becomes severe without architectural interventions like residual connections.

Understanding the computational cost of backpropagation helps explain why it's so efficient and why training is more expensive than inference.

For a single layer with $n$ inputs and $m$ outputs: