Backpropagation Derived Step by Step

Backpropagation is the algorithm that made deep learning possible. It was discovered several times — Linnainmaa's reverse-mode differentiation in 1970, Werbos' application to neural networks in 1974 — and brought to prominence by Rumelhart, Hinton and Williams in 1986. Every time you call loss.backward(), this algorithm runs. Today we derive it completely, by hand, so that it never feels like magic again.

The setting#

Consider a two-layer network for classification with input $\mathbf{x} \in \mathbb{R}^d$, a hidden layer of width $h$ and $K$ output classes:

We want $\frac{\partial L}{\partial\mathbf{W}_1}, \frac{\partial L}{\partial\mathbf{b}_1}, \frac{\partial L}{\partial\mathbf{W}_2}, \frac{\partial L}{\partial\mathbf{b}_2}$.

The key idea#

Define the error signal at each layer, $\boldsymbol{\delta} = \partial L/\partial\mathbf{z}$. Compute it at the output, then propagate it backwards layer by layer using the chain rule. Once you know $\boldsymbol{\delta}$ for a layer, the weight gradients of that layer follow immediately.

Step 1: output layer#

For softmax with cross-entropy, we showed earlier that

where $\mathbf{y}$ is the one-hot label.

Step 2: gradients of the output weights#

Since $\mathbf{z}_2 = \mathbf{W}_2\mathbf{a}_1 + \mathbf{b}_2$, each entry $z_{2,k} = \sum_j W_{2,kj}a_{1,j} + b_{2,k}$, so $\partial z_{2,k}/\partial W_{2,kj} = a_{1,j}$. Therefore

Pattern: weight gradient = (error at the layer's output) × (input to the layer)ᵀ.

Step 3: propagate to the hidden layer#

The hidden activation $\mathbf{a}_1$ influences every output through $\mathbf{W}_2$. Summing over those paths (the multivariable chain rule):

Then pass through the ReLU, whose derivative is 1 where $z_1 > 0$ and 0 elsewhere:

Step 4: gradients of the first layer#

Exactly the same pattern as Step 2:

The general algorithm#

For a network with layers $\mathbf{z}_l = \mathbf{W}_l\mathbf{a}_{l-1} + \mathbf{b}_l$, $\mathbf{a}_l = \phi(\mathbf{z}_l)$:

1. Forward pass — compute and store all $\mathbf{z}_l$ and $\mathbf{a}_l$.
2. Output error — $\boldsymbol{\delta}_L = \partial L/\partial\mathbf{z}_L$.
3. Backward pass — for $l = L, \dots, 1$:

For a mini-batch stored as rows of $\mathbf{X}$, the products become $\mathbf{A}_{l-1}^\top\boldsymbol{\Delta}_l$ summed (or averaged) over the batch.

Why it is efficient#

A naive approach would compute each of the $P$ parameter derivatives separately by perturbation, costing $O(P)$ forward passes. Backpropagation reuses the shared error signals $\boldsymbol{\delta}_l$, computing all gradients in one backward pass whose cost is about two to three times the forward pass. For a model with a billion parameters, that is the difference between impossible and routine. The price is memory: stored activations from the forward pass.

Implementation in NumPy#

import numpy as np

rng = np.random.default_rng(0)
N, d, h, K = 64, 10, 32, 3
X = rng.normal(size=(N, d)); y = rng.integers(0, K, N); Y = np.eye(K)[y]
W1 = rng.normal(0, np.sqrt(2 / d), (d, h)); b1 = np.zeros(h)
W2 = rng.normal(0, np.sqrt(2 / h), (h, K)); b2 = np.zeros(K)

def forward(W1, b1, W2, b2):
    Z1 = X @ W1 + b1; A1 = np.maximum(0, Z1)
    Z2 = A1 @ W2 + b2
    Z2s = Z2 - Z2.max(1, keepdims=True)
    P = np.exp(Z2s) / np.exp(Z2s).sum(1, keepdims=True)
    loss = -np.log(P[np.arange(N), y]).mean()
    return loss, (Z1, A1, P)

def backward(cache):
    Z1, A1, P = cache
    D2 = (P - Y) / N                        # dL/dZ2, averaged over the batch
    dW2, db2 = A1.T @ D2, D2.sum(0)
    D1 = (D2 @ W2.T) * (Z1 > 0)             # dL/dZ1
    dW1, db1 = X.T @ D1, D1.sum(0)
    return dW1, db1, dW2, db2

loss, cache = forward(W1, b1, W2, b2)
grads = backward(cache)

# Gradient check on W1 with central differences
num = np.zeros_like(W1); eps = 1e-5
for i in range(d):
    for j in range(h):
        W1[i, j] += eps; lp, _ = forward(W1, b1, W2, b2)
        W1[i, j] -= 2 * eps; lm, _ = forward(W1, b1, W2, b2)
        W1[i, j] += eps; num[i, j] = (lp - lm) / (2 * eps)
print("relative error:", np.linalg.norm(num - grads[0]) / np.linalg.norm(num + grads[0]))

# Train with plain gradient descent
params = [W1, b1, W2, b2]
for step in range(500):
    loss, cache = forward(*params)
    for p, g in zip(params, backward(cache)):
        p -= 0.5 * g
print("final training loss:", round(loss, 4))

The gradient check should report a relative error around $10^{-8}$ or smaller — confirmation that our derivation and code agree.