The Calculus of Learning: Understanding Gradients and the Chain Rule

Training a neural network is an optimization problem dressed in engineering clothing. The training loop isn’t writing logic or querying data — it’s navigating a high-dimensional surface, looking for a configuration of weights that produces the lowest possible error.

That surface is the loss landscape. Every weight and bias in the model contributes one dimension to this space. A small network with a few thousand parameters already lives in a space too high-dimensional to visualize. A large language model navigates something with billions of dimensions.

Moving downhill in that space requires knowing the slope — not approximately, not by sampling, but exactly, at every point, with respect to every single parameter at once. The mechanism that computes this is the gradient vector. The theorem that makes gradient computation tractable in nested functions is the chain rule. Together, they are backpropagation.

The Gradient Vector

The goal during training is to minimize the loss function $L(\theta)$, where $\theta$ represents every weight and bias in the model and $L$ outputs a single scalar measuring the prediction error.

In single-variable calculus, the derivative $\frac{dy}{dx}$ gives you the slope at a point on a curve. If it’s positive, the function is rising; move left to descend. In multivariable calculus — which is what we need here — there isn’t a single slope. There’s a separate partial derivative for every parameter, measuring how the loss changes when that one parameter is nudged while everything else is held fixed:

$\frac{\partial L}{\partial w_i}$

Collecting these into a single vector gives the gradient:

$\nabla L = \begin{bmatrix} \frac{\partial L}{\partial w_1} \\ \frac{\partial L}{\partial w_2} \\ \vdots \\ \frac{\partial L}{\partial w_n} \end{bmatrix}$

The gradient points in the direction of steepest ascent in the loss landscape. To reduce the loss, we move in the opposite direction, scaled by the learning rate $\alpha$:

$\theta_{\text{new}} = \theta_{\text{old}} - \alpha \nabla L$

This is gradient descent. The learning rate controls step size — too large and updates overshoot, potentially making loss worse; too small and training stalls.

The Chain Rule

The awkward part: $L$ is not a direct function of $w_1$. In a layered network, the loss depends on the output of the final layer, which depends on the output of the second-to-last layer, all the way back. A weight in layer one influences the final loss only indirectly, through every transformation between them.

A 3-layer network computes:

$y = f_3(f_2(f_1(x)))$

To calculate $\frac{\partial L}{\partial w_1}$, we need to follow the signal from $w_1$ forward through three composed functions to the output, then trace it backward to the loss. This is exactly what the chain rule handles.

For a composed function $y = g(h(x))$:

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

The derivative of the outer function times the derivative of the inner function. This extends to arbitrarily deep compositions: each additional layer of nesting adds one more multiplicative factor. Backpropagation is the systematic application of this rule from output back to input — computing local gradients at each layer and multiplying them together.

Step-by-Step Derivation

A concrete example makes this precise. We’ll use a single input $x$, a single hidden neuron with a sigmoid activation, and mean squared error (MSE) loss. This is the smallest network that shows all the moving parts.

Forward Pass

Linear combination: $z = w \cdot x + b$

Sigmoid activation: $a = \sigma(z) = \frac{1}{1 + e^{-z}}$

MSE loss (the $\frac{1}{2}$ factor is a notational convenience that cancels the exponent during differentiation): $L = \frac{1}{2}(a - y)^2$

Backward Pass: Computing $\frac{\partial L}{\partial w}$

By the chain rule:

$\frac{\partial L}{\partial w} = \frac{\partial L}{\partial a} \cdot \frac{\partial a}{\partial z} \cdot \frac{\partial z}{\partial w}$

$\frac{\partial L}{\partial a}$ — Loss sensitivity to the activation output:

$\frac{\partial}{\partial a}\left[\frac{1}{2}(a-y)^2\right] = (a - y)$

$\frac{\partial a}{\partial z}$ — Sigmoid’s derivative:

The sigmoid has a well-known self-referential derivative: $\frac{d\sigma}{dz} = \sigma(z)(1 - \sigma(z)) = a(1 - a)$

$\frac{\partial z}{\partial w}$ — Linear term’s derivative:

$z = w \cdot x + b \implies \frac{\partial z}{\partial w} = x$

Assembled:

$\frac{\partial L}{\partial w} = (a - y) \cdot a(1 - a) \cdot x$

For completeness, the bias gradient follows the same chain but with $\frac{\partial z}{\partial b} = 1$:

$\frac{\partial L}{\partial b} = (a - y) \cdot a(1 - a)$

The bias sees the same error signal as the weight, minus the input scaling factor. In code, this means both $w$ and $b$ need to be updated on every pass — which the implementation below does correctly.

The Vanishing Gradient Problem

Deriving the backward pass manually makes a structural problem in deep networks visible. Look at the sigmoid’s derivative: $a(1 - a)$.

The sigmoid output $a$ lives in $(0, 1)$. The product $a(1-a)$ is maximized when $a = 0.5$, giving $0.25$. This is an upper bound — the gradient that passes through a single sigmoid layer is at most a quarter of what entered. In practice it’s often far less, because sigmoid saturates: when $|z|$ is large, $a$ approaches 0 or 1, and the derivative approaches 0. A saturated neuron passes essentially no gradient at all.

Stack 10 sigmoid layers with favorable (non-saturated) inputs:

$0.25^{10} \approx 9.5 \times 10^{-7}$

That’s the gradient that reaches layer one after propagating back through ten layers under the best-case assumption. In realistic saturated conditions, it’s orders of magnitude smaller. The early layers receive no meaningful update signal, and training stalls. This is the vanishing gradient problem.

ReLU ($f(x) = \max(0, x)$) addresses this because its derivative is 1 for any positive input — the gradient passes through unchanged rather than being multiplied by a fraction. This is why ReLU became the default activation for hidden layers in deep networks.

The tradeoff is the dead neuron problem: if a neuron’s input is consistently negative, it outputs zero, its gradient is zero, and its weights never update. Once dead, it stays dead. Variants like Leaky ReLU ($f(x) = \max(0.01x, x)$) and GELU address this by allowing a small nonzero gradient for negative inputs.

It’s also worth noting that activation choice alone doesn’t fully solve vanishing gradients. Weight initialization matters too — Xavier and He initialization scale the initial weights to keep activations in a reasonable range from the start, preventing saturation before training even begins.

Implementation: Bare-Metal Backpropagation

The code below implements everything from scratch — no autograd, no frameworks. Forward pass, chain rule, weight and bias update. The training target is $y = 1.0$ with input $x = 2.0$; we’re watching the network adjust a single weight and bias to push its output toward 1.

One note on scope: this runs on a single sample, not a mini-batch. Real training uses batches of data to estimate the gradient more reliably and take advantage of parallelism. The math is identical; the batch version averages the gradients across all samples in the batch before updating.

import math

def sigmoid(z: float) -> float:
    return 1.0 / (1.0 + math.exp(-z))

x = 2.0           # Input
y = 1.0           # True target
w = 0.5           # Initial weight
b = 0.0           # Initial bias
lr = 0.1          # Learning rate

print(f"Initial weight: {w:.4f}, bias: {b:.4f}\n")

for epoch in range(5):
    # Forward pass
    z = (w * x) + b
    a = sigmoid(z)
    loss = 0.5 * (a - y) ** 2

    # Backward pass — chain rule
    dL_da = a - y
    da_dz = a * (1.0 - a)      # sigmoid derivative
    dz_dw = x
    dz_db = 1.0

    dL_dw = dL_da * da_dz * dz_dw
    dL_db = dL_da * da_dz * dz_db

    # Gradient descent step
    w -= lr * dL_dw
    b -= lr * dL_db

    print(
        f"Epoch {epoch+1} | Loss: {loss:.5f} | "
        f"dL/dw: {dL_dw:.5f} | w: {w:.5f} | b: {b:.5f}"
    )

Initial weight: 0.5000, bias: 0.0000

Epoch 1 | Loss: 0.03532 | dL/dw: -0.10499 | w: 0.51050 | b: 0.00525
Epoch 2 | Loss: 0.03328 | dL/dw: -0.10091 | w: 0.52059 | b: 0.01029
Epoch 3 | Loss: 0.03138 | dL/dw: -0.09696 | w: 0.53028 | b: 0.01513
Epoch 4 | Loss: 0.02960 | dL/dw: -0.09314 | w: 0.53960 | b: 0.01980
Epoch 5 | Loss: 0.02794 | dL/dw: -0.08945 | w: 0.54854 | b: 0.02427

The gradient is negative throughout, meaning increasing $w$ would decrease the loss — so the update subtracts a negative value, pushing $w$ upward. As $w$ increases and the network’s output approaches the target, the error $(a - y)$ shrinks, which pulls the gradient magnitude down with it. Smaller gradient, smaller update, slower movement — the system is decelerating as it approaches a lower-error region.

Note that the loss is decreasing but nowhere near zero after five epochs. That’s expected: a single neuron with sigmoid activation has limited expressive capacity, and five gradient steps on one sample isn’t training — it’s an illustration. A real network reaches useful loss values through thousands of steps over large datasets, with an optimizer (Adam, AdamW) that adapts the learning rate per parameter rather than using a fixed global rate.

A Note on “Finding the Minimum”

It’s tempting to describe gradient descent as finding the minimum of the loss function, but the loss landscape of a non-convex network has no global minimum in any practical sense — it has an enormous number of local minima, saddle points, and flat plateaus. What gradient descent actually does is find a minimum that generalizes well enough to be useful.

Empirically, this works better than it should theoretically. Large networks seem to have many local minima of similar quality, and the saddle point problem turns out to be less severe in very high dimensions than in low dimensions. But “gradient descent finds the minimum” is an oversimplification. It finds a low-loss region, and whether that region generalizes to unseen data depends on factors well beyond the optimizer: architecture, regularization, data quality, and learning rate scheduling.

Understanding backpropagation at the calculus level matters because it makes those factors interpretable. Vanishing gradients aren’t a mysterious failure mode — they’re $0.25^n$ applied to saturated activations. Exploding gradients are the inverse: gradients with magnitude greater than 1 multiplied through many layers. Gradient clipping, batch normalization, careful initialization — these are all direct responses to what the math predicts will go wrong.