Backpropagation by Hand: A Two-Layer MLP in NumPy with a Gradient Check

Automatic differentiation has made it possible to train networks without ever writing a gradient, which is also why so many people are unsure what a gradient check is for. This post goes the other way: derive the backward pass of a two-layer perceptron on paper, implement it in about forty lines of NumPy, verify it against finite differences, and only then train it. The network is deliberately small (2 inputs, 16 hidden units, 1 output), but nothing in the derivation changes for larger ones.

The computational graph. Forward arrows compute values left to right; backward arrows (blue) carry the upstream gradient right to left, multiplying by one local Jacobian per node.

The forward pass

With a batch $X \in R^{N \times 2}$ , labels $y \in {0, 1}^{N}$ , weights $W_{1} \in R^{2 \times 16}$ , $W_{2} \in R^{16 \times 1}$ and biases $b_{1}, b_{2}$ :

\begin{aligned} Z_{1} & = X W_{1} + b_{1}, & H & = \tanh (Z_{1}), \\ Z_{2} & = H W_{2} + b_{2}, & P & = σ (Z_{2}) = \frac{1}{1 + e^{- Z_{2}}}, \\ L & = - \frac{1}{N} \sum_{i = 1}^{N} [y_{i} \log p_{i} + (1 - y_{i}) \log (1 - p_{i})] . \end{aligned}

import numpy as np

def init(n_in, n_h, n_out, rng):
    return {"W1": rng.normal(0, 1 / np.sqrt(n_in), (n_in, n_h)), "b1": np.zeros(n_h),
            "W2": rng.normal(0, 1 / np.sqrt(n_h), (n_h, n_out)), "b2": np.zeros(n_out)}

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

def forward(P, X):
    Z1 = X @ P["W1"] + P["b1"]
    H = np.tanh(Z1)
    Z2 = H @ P["W2"] + P["b2"]
    return Z1, H, Z2, sigmoid(Z2)

def loss(Yhat, y):
    y = y.reshape(-1, 1)                       # match Yhat's (N, 1) shape, see the warning below
    eps = 1e-12
    return -np.mean(y * np.log(Yhat + eps) + (1 - y) * np.log(1 - Yhat + eps))

The backward pass, one node at a time

Backpropagation is the chain rule applied in reverse topological order. Each node receives the gradient of the loss with respect to its output ("upstream gradient"), multiplies by its own local Jacobian, and passes the result to its inputs. For matrix-shaped nodes the rule that avoids ever forming a Jacobian is: if $Y = X W$ and the upstream gradient is $\partial L / \partial Y$ , then $\partial L / \partial X = (\partial L / \partial Y) W^{⊤}$ and $\partial L / \partial W = X^{⊤} (\partial L / \partial Y)$ . Shapes make it hard to get wrong: each gradient has the shape of the thing it is a gradient of.

The first step is the one worth memorising. The derivative of binary cross-entropy through the sigmoid collapses to a difference:

\frac{\partial L}{\partial Z_{2}} = \frac{1}{N} (P - y) .

Why the sigmoid and the cross-entropy cancel

For one sample, $ℓ = - [y \log p + (1 - y) \log (1 - p)]$ with $p = σ (z)$ . Using $σ^{'} (z) = σ (z) (1 - σ (z)) = p (1 - p)$ :

\frac{\partial ℓ}{\partial z} = \frac{\partial ℓ}{\partial p} \frac{\partial p}{\partial z} = (- \frac{y}{p} + \frac{1 - y}{1 - p}) p (1 - p) = - y (1 - p) + (1 - y) p = p - y .

The $1 / N$ comes from the mean over the batch. The same cancellation happens for softmax with categorical cross-entropy, which is why both pairs are the default choice: the gradient never saturates even when $p$ does.

From there the chain is mechanical:

\begin{aligned} \frac{\partial L}{\partial W_{2}} & = H^{⊤} \frac{\partial L}{\partial Z_{2}}, & \frac{\partial L}{\partial b_{2}} & = \sum_{i} \frac{\partial L}{\partial Z_{2}} |_{i}, \\ \frac{\partial L}{\partial H} & = \frac{\partial L}{\partial Z_{2}} W_{2}^{⊤}, & \frac{\partial L}{\partial Z_{1}} & = \frac{\partial L}{\partial H} ⊙ (1 - H^{2}), \\ \frac{\partial L}{\partial W_{1}} & = X^{⊤} \frac{\partial L}{\partial Z_{1}}, & \frac{\partial L}{\partial b_{1}} & = \sum_{i} \frac{\partial L}{\partial Z_{1}} |_{i} . \end{aligned}

The $(1 - H^{2})$ factor is $\tanh^{'} (z) = 1 - \tanh^{2} (z)$ , evaluated using the already computed activation, which is why the forward pass keeps $H$ around.

def backward(P, X, y, cache):
    Z1, H, Z2, Yhat = cache
    N = X.shape[0]
    dZ2 = (Yhat - y[:, None]) / N              # (N, 1)
    g = {"W2": H.T @ dZ2, "b2": dZ2.sum(0)}
    dH = dZ2 @ P["W2"].T                       # (N, 16)
    dZ1 = dH * (1 - H ** 2)
    g["W1"] = X.T @ dZ1
    g["b1"] = dZ1.sum(0)
    return g

The gradient check

A gradient check compares the analytic gradient with a numerical one. For every scalar parameter $θ_{k}$ , perturb it by $\pm ϵ$ and compute a central difference:

\frac{\partial L}{\partial θ_{k}} \approx \frac{L (θ_{k} + ϵ) - L (θ_{k} - ϵ)}{2 ϵ}, rel. err. = \frac{| g_{analytic} - g_{numeric} |}{| g_{analytic} | + | g_{numeric} |} .

The recipe that avoids false alarms:

Use float64 and a tiny network (here 2–5–1 on 8 samples) so the check runs in milliseconds.
Use central differences with $ϵ \approx 10^{- 5}$ ; one-sided differences have error $O (ϵ)$ instead of $O (ϵ^{2})$ .
Use a relative error; absolute errors mean nothing when gradients span several orders of magnitude.
Check before adding any regularization or dropout, which make the loss non-deterministic or non-smooth.

rng = np.random.default_rng(0)
Xs, ys = moons(8, 0.1, rng)                    # a tiny two-moons sample
Ps = init(2, 5, 1, rng)
g = backward(Ps, Xs, ys, forward(Ps, Xs))
eps = 1e-5
for k in Ps:
    num = np.zeros_like(Ps[k])
    it = np.nditer(Ps[k], flags=["multi_index"])
    while not it.finished:
        idx, old = it.multi_index, Ps[k][it.multi_index]
        Ps[k][idx] = old + eps; lp = loss(forward(Ps, Xs)[3], ys)
        Ps[k][idx] = old - eps; lm = loss(forward(Ps, Xs)[3], ys)
        Ps[k][idx] = old
        num[idx] = (lp - lm) / (2 * eps)
        it.iternext()
    rel = np.abs(g[k] - num) / np.maximum(np.abs(g[k]) + np.abs(num), 1e-12)
    print(f"{k}: max relative error = {rel.max():.3e}")

The output of this run:

Parameter	Max relative error
`W1`	6.3 × 10⁻¹⁰
`b1`	7.9 × 10⁻¹⁰
`W2`	1.6 × 10⁻¹⁰
`b2`	3.0 × 10⁻¹⁰

Anything around $10^{- 7}$ or better in double precision means the analytic gradient is right. Errors around $10^{- 2}$ usually mean a transposed matrix or a missing factor; an error of exactly $1.0$ means one of the two gradients is zero or has the wrong sign.

DANGER

The first version of this script had a relative error of exactly 1.0 on W1, b1 and W2 while b2 passed. The cause was not in the backward pass: the loss function received Yhat of shape (N, 1) and y of shape (N,), and NumPy silently broadcast them to an (N, N) matrix before taking the mean. The network still "trained" to 99 % accuracy, with a loss that went up. A gradient check catches the kind of bug that accuracy curves hide.

Training

With the gradient verified, training is full-batch gradient descent with learning rate 0.5 for 3000 steps on 400 two-moons points, with another 400 points held out.

X, y = moons(400, 0.15, rng)
P = init(2, 16, 1, rng)
for step in range(3000):
    cache = forward(P, X)
    g = backward(P, X, y, cache)
    for k in P:
        P[k] -= 0.5 * g[k]

Training and validation loss and validation accuracy over 3000 steps — Loss and validation accuracy during training. The plateau around 0.22 between steps 100 and 300 is the network fitting a nearly linear boundary before the hidden units bend it.

Step	Train BCE	Val BCE	Val accuracy
0	0.6465	0.4704	0.840
100	0.2282	0.2124	0.917
500	0.1067	0.1032	0.975
1000	0.0341	0.0346	0.990
2000	0.0174	0.0283	0.990
2999	0.0135	0.0304	0.990

The validation loss starts creeping up after step 2000 while the accuracy stays flat: the network is becoming over-confident on points near the boundary rather than misclassifying new ones, which is the first, mild form of overfitting.

Decision function of the trained network over the two-moons validation set — The learned decision function on the validation set. Sixteen tanh units are enough to bend the boundary around both moons; the few misclassified points sit in the overlap of the noisy crescents.

What generalizes from this exercise

Every layer type has the same contract: given the upstream gradient, return the gradient with respect to inputs and parameters. Frameworks implement exactly this per operation.
Shape bookkeeping is most of the work, and most of the bugs. Writing the expected shape next to each line, as in the backward function, is cheap insurance.
A gradient check is the unit test of a backward pass. It costs minutes to write and finds bugs that no learning curve will.

Backpropagation by Hand: A Two-Layer MLP in NumPy with a Gradient Check

The forward pass

The backward pass, one node at a time

The gradient check

DANGER

Training

What generalizes from this exercise

Further reading

Cite this post

Keep reading