Backpropagation — the chain rule walking backward

The last lesson left us needing the gradient ∇L and unable to find it: a real network has millions of weights, and differentiating that tangle by hand — or one partial at a time — is hopeless. Backpropagation is the way out. It gets called “the algorithm that powers deep learning,” and the name lends the field a mystery it does not deserve. Backprop is nothing but the chain rule of calculus applied to a chain of operations, computed backward so that each local derivative is done exactly once. Trace it through on a tiny network and the magic quietly evaporates.

The chain rule (one-line recap)

If y = f(g(x)), then:

dy/dx = f'(g(x)) · g'(x)

Computing a derivative through nested functions means multiplying local derivatives. That’s the entire idea.

For a longer chain y = f₃(f₂(f₁(x))):

dy/dx = f₃'(f₂(f₁(x))) · f₂'(f₁(x)) · f₁'(x)

A neural network is just a long chain (with branches) of such functions. The chain rule, applied to that graph, is backpropagation.

Why go backward? If you computed each gradient starting from the input and walking forward, you’d re-traverse the rest of the network for every weight — cost grows with depth. Going backward, each local derivative is computed exactly once and reused — dynamic programming on the graph.

A tiny network — by hand

Let’s do the smallest non-trivial example: one input, one hidden unit, one output. We’ll use a sigmoid activation and squared-error loss. The parameters are w1 and w2. We want ∂L/∂w1 and ∂L/∂w2.

Before grinding the algebra by hand, play with the mechanism. Below is the same idea on a single neuron — z = w·x + b, a = σ(z), L = (a − y)². Run the forward pass to fill in each value, then the backward pass to watch every gradient form as downstream gradient × the local derivative on the edge.

Forward pass — compute the values

For x = 1.0, y = 0.7, w1 = 0.5, w2 = 0.8:

z1 = w1 * x      = 0.5
h  = σ(z1)       = σ(0.5)     ≈ 0.6225
z2 = w2 * h      = 0.8 * 0.6225 ≈ 0.4980
L  = (z2 - y)²   = (0.4980 - 0.7)² ≈ 0.0408

Backward pass — chain rule, right to left

Start at the loss. Compute ∂L/∂z2, then propagate backward.

∂L/∂z2 = 2 * (z2 - y)                       (derivative of (z2 - y)²)
∂z2/∂h = w2                                 (z2 = w2 * h)
∂h/∂z1 = σ'(z1) = σ(z1) * (1 - σ(z1))       (sigmoid derivative)
∂z1/∂w1 = x                                 (z1 = w1 * x)
∂z2/∂w2 = h                                 (z2 = w2 * h)

Now chain them:

∂L/∂w2 = (∂L/∂z2) · (∂z2/∂w2)
       = 2(z2 - y) · h

∂L/∂w1 = (∂L/∂z2) · (∂z2/∂h) · (∂h/∂z1) · (∂z1/∂w1)
       = 2(z2 - y) · w2 · σ'(z1) · x

Each term in the chain is a local derivative. You compute them, then multiply. That’s it.

In code, step by step

import numpy as np

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

# Input + target
x = 1.0
y = 0.7

# Parameters (starting values)
w1 = 0.5
w2 = 0.8

# --- Forward pass ---
z1 = w1 * x
h  = sigmoid(z1)
z2 = w2 * h
L  = (z2 - y) ** 2

print(f"Forward:  z1={z1:.4f}  h={h:.4f}  z2={z2:.4f}  L={L:.5f}")

# --- Backward pass ---
dL_dz2 = 2.0 * (z2 - y)             # ∂L/∂z2
dz2_dh = w2                         # ∂z2/∂h
dh_dz1 = h * (1.0 - h)              # σ'(z1) = σ(1-σ)
dz1_dw1 = x                         # ∂z1/∂w1
dz2_dw2 = h                         # ∂z2/∂w2

dL_dw2 = dL_dz2 * dz2_dw2
dL_dw1 = dL_dz2 * dz2_dh * dh_dz1 * dz1_dw1

print(f"Gradients:  ∂L/∂w1={dL_dw1:.5f}   ∂L/∂w2={dL_dw2:.5f}")

# --- One gradient-descent step ---
lr = 0.5
w1 -= lr * dL_dw1
w2 -= lr * dL_dw2
print(f"Updated:  w1={w1:.4f}  w2={w2:.4f}")

# Did the loss go down? Re-run the forward pass.
z1 = w1 * x; h = sigmoid(z1); z2 = w2 * h
L_new = (z2 - y) ** 2
print(f"New loss: {L_new:.5f}  (was {L:.5f})")

Forward:  z1=0.5000  h=0.6225  z2=0.4980  L=0.04082
Gradients:  ∂L/∂w1=-0.07597   ∂L/∂w2=-0.25151
Updated:  w1=0.5380  w2=0.9258
New loss: 0.01335  (was 0.04082)

The numbers match the hand-derivation exactly, and the loss fell from 0.041 to 0.013 in a single step — one backprop pass genuinely moved the model toward the target. Loop this thousands of times and you have trained a network.

A real (tiny) MLP — batched, vectorized

Now scale it up: a 2-feature input, a 4-neuron hidden layer, and a 1-output regression head. We’ll vectorize so the math handles a whole batch at once.

import numpy as np

rng = np.random.default_rng(0)

# Tiny dataset: y = sum of two features + a bit of noise
n = 64
X = rng.normal(size=(n, 2))
y = (X[:, 0] + X[:, 1] + 0.1 * rng.normal(size=n)).reshape(-1, 1)

# Parameters: input -> hidden (4) -> output (1)
W1 = rng.normal(size=(2, 4)) * 0.5
b1 = np.zeros((1, 4))
W2 = rng.normal(size=(4, 1)) * 0.5
b2 = np.zeros((1, 1))

def relu(z):     return np.maximum(0, z)
def d_relu(z):   return (z > 0).astype(float)

lr = 0.05
for step in range(2000):
    # ---- forward ----
    Z1 = X @ W1 + b1
    H  = relu(Z1)
    Z2 = H @ W2 + b2
    err = Z2 - y
    L = (err ** 2).mean()

    # ---- backward (chain rule) ----
    # dL/dZ2 = 2 * err / n  (because loss is mean of err^2)
    dZ2 = 2 * err / n
    dW2 = H.T @ dZ2
    db2 = dZ2.sum(axis=0, keepdims=True)

    dH  = dZ2 @ W2.T
    dZ1 = dH * d_relu(Z1)
    dW1 = X.T @ dZ1
    db1 = dZ1.sum(axis=0, keepdims=True)

    # ---- step ----
    W1 -= lr * dW1; b1 -= lr * db1
    W2 -= lr * dW2; b2 -= lr * db2

    if step % 400 == 0:
        print(f"step {step:>4}  loss={L:.4f}")

# Final prediction sanity check on a fresh point
x_test = np.array([[1.0, 2.0]])
pred = relu(x_test @ W1 + b1) @ W2 + b2
print(f"\nPredict y for [1, 2]: {pred[0, 0]:.3f}  (truth ~ 3.0)")

step    0  loss=1.3148
step  400  loss=0.0095
step  800  loss=0.0092
step 1200  loss=0.0090
step 1600  loss=0.0090

Predict y for [1, 2]: 2.952  (truth ~ 3.0)

Two-line forward pass, four-line backward pass — and the loss falls from 1.31 to 0.009, with the trained net predicting 2.952 for [1, 2] against a true answer near 3.0. This is a neural network from scratch.

Notice the pattern. Each backward step:

1. Takes the gradient from the next layer (dZ2, then dH).
2. Splits it into a gradient for the layer’s weights and the gradient to pass on backward to the previous layer.
3. The weight gradient is (input to layer).T @ (gradient out).
4. The “pass-on” gradient is (gradient out) @ W.T, possibly times the activation derivative.

If you internalize that pattern, every layer’s backward formula in any framework will look familiar.

In one breath

Backprop is the chain rule applied to a network’s computation graph and walked backward, so each local derivative is computed once and reused (going forward would re-walk the network for every weight — backward is dynamic programming on the graph). The forward pass computes and stashes the intermediate values; the backward pass starts at the loss and multiplies local derivatives edge by edge, each gradient forming as downstream gradient × local derivative. The layer pattern repeats: the weight gradient is (layer input).T @ (gradient out) and the pass-on gradient is (gradient out) @ W.T times the activation derivative — every formula falling straight out of shape arithmetic. PyTorch’s .backward() runs exactly this for you; that is all “autograd” means.

Practice

A question to carry forward

Backprop hands us the gradient of a scalar loss — one number out, a vector of slopes back. But two cracks have already shown through. First, a network layer does not output one number; it outputs a vector, and “the derivative of a vector with respect to a vector” is something richer than a gradient. Second, and quieter: all through gradient descent we watched one coordinate sprint to its answer while another crawled, purely because the loss curved more steeply in one direction than the other — and the gradient, a first-order object, is blind to curvature entirely.

Here is the thread onward: when a function maps many inputs to many outputs, what is the full table of its first derivatives — the Jacobian — and what is the matrix of its second derivatives, the Hessian, that finally measures curvature? How does a Taylor series stitch slope and curvature into one local portrait of a function, and why does that portrait explain both why plain gradient descent struggles in a steep-sided ravine and how smarter optimizers find their way out?