Backpropagation by hand, then in autograd

This lesson is for anyone who has trained a network with loss.backward() and wants to know what that line actually does. You will push one input through a two-layer network with pencil-sized numbers, differentiate the loss with respect to every weight using nothing but the chain rule, write the same arithmetic in numpy, prove it correct with finite differences, and then run the identical network through PyTorch autograd. The gradients agree to the last few bits of a double, and by the end you will understand why they have to.

The training loop itself was covered in the beginner lesson on loss and gradients, and the shape of a small network in your first neural network. Here the focus is only the gradient.

The network and the numbers

The network has two inputs, two hidden units with a tanh non-linearity, and one linear output. The loss is half the squared error against a scalar target. Half is there so the derivative is clean.

The forward computation, written as four lines:

z1 = W1 x + b1          W1 is 2 by 2, b1 has 2 entries
h  = tanh(z1)
ŷ  = W2 h + b2          W2 has 2 entries, b2 is a scalar
L  = ½ (ŷ − y)²

The concrete values, chosen so the first steps come out as round numbers:

Symbol	Value
x	[0.5, −1.0]
y	1.0
W1	rows [0.2, −0.4] and [0.6, 0.1]
b1	[0.1, −0.3]
W2	[0.8, −0.5]
b2	0.2

There are nine trainable numbers in total: four in W1, two in b1, two in W2 and one in b2. The job of backpropagation is to produce nine partial derivatives, one per parameter.

Forward pass

Work through the four lines with the values above.

The first hidden unit sees 0.2 × 0.5 + (−0.4) × (−1.0) + 0.1, which is 0.1 + 0.4 + 0.1, so z1 is 0.6 for that unit. The second sees 0.6 × 0.5 + 0.1 × (−1.0) − 0.3, which is 0.3 − 0.1 − 0.3, so −0.1. The pre-activation vector is z1 = [0.6, −0.1].

Applying tanh gives h = [0.5370, −0.0997]. The output is 0.8 × 0.5370 + (−0.5) × (−0.0997) + 0.2 = 0.4296 + 0.0498 + 0.2 = 0.6795.

The error is ŷ − y = −0.3205, and the loss is half of its square, 0.0514. Write the error down as e, because it is the first thing you need going backwards and it appears in every gradient.

The chain rule, backwards

The chain rule says that if L depends on ŷ and ŷ depends on W2, then the derivative of L with respect to W2 is the derivative of L with respect to ŷ multiplied by the derivative of ŷ with respect to W2. Backpropagation is nothing more than applying that rule once per operation, starting at the loss and moving towards the inputs, reusing every intermediate derivative instead of recomputing it.

Output layer

The derivative of ½ (ŷ − y)² with respect to ŷ is ŷ − y, which is e = −0.3205. This single number is the seed for everything else.

The output ŷ = W2 h + b2 is linear in W2, so the derivative of ŷ with respect to each entry of W2 is the matching entry of h. Multiply by the seed:

dL/dW2 = e · h  = −0.3205 × [0.5370, −0.0997] = [−0.1721, 0.0319]
dL/db2 = e      = −0.3205

Two of the nine gradients are done. The output is also linear in h, and the derivative of ŷ with respect to h is W2, so the gradient that must be handed back to the hidden layer is:

dL/dh = e · W2 = −0.3205 × [0.8, −0.5] = [−0.2564, 0.1603]

Hidden layer

Now h = tanh(z1), applied element by element. The derivative of tanh(z) is 1 − tanh(z)², and tanh(z1) is h, which you already have. So the local derivative of each hidden unit is 1 − h². For the two units that is 1 − 0.2884 = 0.7116 and 1 − 0.0099 = 0.9901. Multiply element-wise:

dL/dz1 = dL/dh ⊙ (1 − h²) = [−0.2564 × 0.7116, 0.1603 × 0.9901] = [−0.1825, 0.1587]

Finally z1 = W1 x + b1 is linear in W1 and b1, exactly as the output layer was linear in W2 and b2. Each row of W1 feeds one hidden unit, so the gradient of row i is dL/dz1 for unit i times the input vector x. That is an outer product.

dL/dW1 = dL/dz1 ⊗ x = rows [−0.1825 × 0.5, −0.1825 × −1.0] and [0.1587 × 0.5, 0.1587 × −1.0]
       = rows [−0.0912, 0.1825] and [0.0793, −0.1587]
dL/db1 = dL/dz1 = [−0.1825, 0.1587]

All nine numbers are in hand. Notice the pattern: every linear layer contributes gradient equal to the incoming gradient times its input, every non-linearity contributes an element-wise factor, and the incoming gradient for a layer is always the previous layer’s incoming gradient multiplied by that layer’s weights. Nothing was differentiated twice.

Computational graph of the two-layer network. Forward values flow left to right along the top row. The gradient of the loss with respect to each intermediate value flows right to left along the bottom row, and each box multiplies the incoming gradient by its own local derivative.

The same thing in numpy

The code below is the hand calculation, line for line. The backward section runs in the reverse order of the forward section, and each line uses only values that are already available.

import numpy as np

np.set_printoptions(precision=4, suppress=True)

x = np.array([0.5, -1.0])
y = 1.0
W1 = np.array([[0.2, -0.4],
               [0.6, 0.1]])
b1 = np.array([0.1, -0.3])
W2 = np.array([0.8, -0.5])
b2 = 0.2

# forward
z1 = W1 @ x + b1
h = np.tanh(z1)
y_hat = W2 @ h + b2
loss = 0.5 * (y_hat - y) ** 2

# backward, in reverse order of the forward pass
d_y_hat = y_hat - y
d_W2 = d_y_hat * h
d_b2 = d_y_hat
d_h = d_y_hat * W2
d_z1 = d_h * (1.0 - h ** 2)
d_W1 = np.outer(d_z1, x)
d_b1 = d_z1

print("z1    ", z1)
print("h     ", h)
print("y_hat ", round(float(y_hat), 4))
print("loss  ", round(float(loss), 4))
print("dW2   ", d_W2)
print("db2   ", round(float(d_b2), 4))
print("dz1   ", d_z1)
print("dW1   ", d_W1)
print("db1   ", d_b1)

The printed values match the hand calculation (numpy pads the columns, so spacing may differ):

z1     [ 0.6 -0.1]
h      [ 0.537  -0.0997]
y_hat  0.6795
loss   0.0514
dW2    [-0.1721  0.0319]
db2    -0.3205
dz1    [-0.1825  0.1587]
dW1    [[-0.0912  0.1825]
 [ 0.0793 -0.1587]]
db1    [-0.1825  0.1587]

Shapes are the usual source of bugs when this is scaled up to batches. The tensor shapes and broadcasting reference covers the rules, and the batched version later in this lesson shows how the outer product becomes a matrix product.

Gradient checking with finite differences

The hand derivation could be wrong in a sign or a transpose and the code would still run. The only independent check is to nudge each parameter by a tiny amount, measure how the loss changes, and compare that slope with the analytic gradient. The central difference (L(p + ε) − L(p − ε)) / 2ε is accurate to order ε², so with ε = 1e-5 in double precision you expect agreement to roughly ten decimal places.

import numpy as np

x = np.array([0.5, -1.0])
y = 1.0
W1 = np.array([[0.2, -0.4], [0.6, 0.1]])
b1 = np.array([0.1, -0.3])
W2 = np.array([0.8, -0.5])
b2 = np.array([0.2])


def forward_backward(W1, b1, W2, b2):
    z1 = W1 @ x + b1
    h = np.tanh(z1)
    y_hat = W2 @ h + b2[0]
    loss = 0.5 * (y_hat - y) ** 2
    d_y_hat = y_hat - y
    grads = {"W2": d_y_hat * h, "b2": np.array([d_y_hat])}
    d_z1 = d_y_hat * W2 * (1.0 - h ** 2)
    grads["W1"] = np.outer(d_z1, x)
    grads["b1"] = d_z1
    return loss, grads


params = {"W1": W1, "b1": b1, "W2": W2, "b2": b2}
loss, grads = forward_backward(**params)

eps = 1e-5
worst = 0.0
for name, p in params.items():
    numeric = np.zeros_like(p)
    for idx in np.ndindex(p.shape):
        saved = p[idx]
        p[idx] = saved + eps
        loss_plus, _ = forward_backward(**params)
        p[idx] = saved - eps
        loss_minus, _ = forward_backward(**params)
        p[idx] = saved
        numeric[idx] = (loss_plus - loss_minus) / (2 * eps)
    rel = np.abs(numeric - grads[name]) / (np.abs(numeric) + np.abs(grads[name]) + 1e-12)
    worst = max(worst, float(rel.max()))
    print(f"{name}: analytic {grads[name].ravel().round(6)}  numeric {numeric.ravel().round(6)}")

print("worst relative error:", f"{worst:.1e}")

Each line shows the analytic and numeric gradients agreeing to six decimal places, and the worst relative error prints as a number of order 1e-10. A sign error in the hidden layer would push that to order 1, and a missing (1 − h²) factor to order 0.1, so the check is not subtle.

The same network in PyTorch autograd

Now build the identical network from tensors and let PyTorch differentiate it. The only new ingredient is requires_grad=True on the four parameters, which tells autograd to record every operation that touches them.

import numpy as np
import torch

torch.set_printoptions(precision=4, sci_mode=False)

x = torch.tensor([0.5, -1.0], dtype=torch.float64)
y = torch.tensor(1.0, dtype=torch.float64)
W1 = torch.tensor([[0.2, -0.4], [0.6, 0.1]], dtype=torch.float64, requires_grad=True)
b1 = torch.tensor([0.1, -0.3], dtype=torch.float64, requires_grad=True)
W2 = torch.tensor([0.8, -0.5], dtype=torch.float64, requires_grad=True)
b2 = torch.tensor(0.2, dtype=torch.float64, requires_grad=True)

z1 = W1 @ x + b1
h = torch.tanh(z1)
y_hat = W2 @ h + b2
loss = 0.5 * (y_hat - y) ** 2

print("loss           ", round(loss.item(), 4))
print("grad_fn of loss", type(loss.grad_fn).__name__)
print("grad_fn of h   ", type(h.grad_fn).__name__)

loss.backward()

print("dW1", W1.grad)
print("db1", b1.grad)
print("dW2", W2.grad)
print("db2", b2.grad)

# the numpy backward pass from earlier, for a direct comparison
d_y_hat = y_hat.item() - 1.0
h_np = h.detach().numpy()
d_z1 = d_y_hat * np.array([0.8, -0.5]) * (1.0 - h_np ** 2)
print("W1 gradients match numpy:", np.allclose(W1.grad.numpy(), np.outer(d_z1, [0.5, -1.0]), atol=1e-12))

The loss prints as 0.0514, the two grad_fn names print as MulBackward0 and TanhBackward0, the four .grad tensors hold the same nine numbers as the numpy block, and the final line prints True.

Before backward() is called, every intermediate tensor already carries a grad_fn. That attribute is the record autograd kept while the forward pass ran. The loss was produced by a multiplication, hence MulBackward0, and h by tanh, hence TanhBackward0. Each of those objects knows its inputs and knows how to compute its own local derivative.

Why the numbers match

Autograd is not doing something cleverer than the hand calculation. It is doing exactly the hand calculation, with the bookkeeping automated.

During the forward pass, every operation on a tensor that requires a gradient appends a node to a graph. The node stores whatever the backward formula will need: for tanh it stores the output h, because the derivative is 1 − h². For a matrix product it stores the other operand. This is the same information you wrote down on paper while going forwards.

When backward() runs, it starts at the loss with a seed gradient of 1 and walks the graph in reverse topological order. At each node it calls that node’s backward function, which takes the gradient arriving from above and returns the gradient with respect to each input. The tanh node receives dL/dh and returns dL/dh ⊙ (1 − h²). The matrix product node receives dL/dz1 and returns the outer product with x for the weight input and the product with W1 for the vector input. When a gradient reaches a leaf tensor with requires_grad=True, it is added into that tensor’s .grad field.

Three consequences follow from this description and all three matter in practice.

It is reverse mode. The work is proportional to one backward sweep regardless of how many parameters there are, because the sweep carries a single vector of gradients and multiplies it by local Jacobians. Computing the full Jacobian would cost one sweep per output, which is why forward-mode differentiation is the wrong tool for a scalar loss over millions of weights.
Gradients accumulate. The .grad field is added to, not overwritten, because a parameter can feed several nodes and their contributions must sum. That is also why a training loop calls optimizer.zero_grad() or sets grad to None before each step. Forgetting it adds this batch’s gradient to last batch’s, and the loss curve becomes mysteriously bumpy. The PyTorch training loop cheat sheet has the standard ordering.
The graph is freed after use. By default backward() releases the saved tensors, so calling it twice on the same loss raises an error unless you pass retain_graph=True. Rebuilding the graph on every forward pass is what lets control flow such as if and for work naturally: the graph is whatever the code actually executed.

Because the tensors above were created in float64, the two implementations agree to about fifteen significant figures. In the default float32 they still agree, but only to about six, and that is expected rounding rather than an error in either method.

Training both side by side

The hand example used one input. Real training uses a batch, and the only change to the backward pass is that outer products become matrix products and bias gradients sum over the batch. The block below trains the numpy network and an autograd copy on the same 64 synthetic points, from the same initial weights, with plain gradient descent. If the numpy backward pass is right, the two loss columns are identical.

import numpy as np
import torch

rng = np.random.default_rng(0)
X = rng.uniform(-1.0, 1.0, size=(64, 2))
Y = X[:, 0] * X[:, 1] + 0.05 * rng.normal(size=64)

hidden = 8
W1 = rng.normal(scale=0.5, size=(hidden, 2))
b1 = np.zeros(hidden)
W2 = rng.normal(scale=0.5, size=hidden)
b2 = 0.0
lr = 0.1

# torch copies of the same starting point, in float64 so the arithmetic agrees
tW1 = torch.tensor(W1, requires_grad=True)
tb1 = torch.tensor(b1, requires_grad=True)
tW2 = torch.tensor(W2, requires_grad=True)
tb2 = torch.tensor(b2, dtype=torch.float64, requires_grad=True)
tX = torch.tensor(X)
tY = torch.tensor(Y)

for step in range(1001):
    # numpy: forward, backward, update
    Z1 = X @ W1.T + b1                      # (64, hidden)
    H = np.tanh(Z1)
    y_hat = H @ W2 + b2                     # (64,)
    loss = np.mean((y_hat - Y) ** 2)
    d_y_hat = 2.0 * (y_hat - Y) / len(Y)    # derivative of the mean squared error
    dW2 = H.T @ d_y_hat
    db2 = d_y_hat.sum()
    dZ1 = np.outer(d_y_hat, W2) * (1.0 - H ** 2)
    dW1 = dZ1.T @ X
    db1 = dZ1.sum(axis=0)
    W1 -= lr * dW1
    b1 -= lr * db1
    W2 -= lr * dW2
    b2 -= lr * db2

    # torch: same forward, gradients from autograd, same update
    t_loss = torch.mean((torch.tanh(tX @ tW1.T + tb1) @ tW2 + tb2 - tY) ** 2)
    t_loss.backward()
    with torch.no_grad():
        for p in (tW1, tb1, tW2, tb2):
            p -= lr * p.grad
            p.grad = None

    if step % 200 == 0:
        print(f"step {step:4d}  numpy {loss:.6f}  torch {t_loss.item():.6f}")

print("max |W1 difference|:", f"{np.abs(W1 - tW1.detach().numpy()).max():.1e}")

The two columns print the same six-decimal loss on every line, the loss falls steadily over the run, and the final weight difference is noise in the last bits of a double, far below 1e-10. The numpy code has no grad_fn, no graph and no backward(), yet it reaches the same weights because it encodes the same chain rule by hand.

Look at the line that computes d_y_hat. It differs from the single-example version only because the loss is now a mean, so the factor 2/N appears. Every other line is the earlier line with a batch dimension, which is exactly what the deep learning hub means when it says the backward pass of a layer is defined by two matrix products. On PyTorch those products are what nn.Linear registers as its backward function, and a model with fifty layers is this lesson’s two-layer calculation applied fifty times.

Exercise

Replace the hidden non-linearity with the sigmoid function σ(z) = 1 / (1 + exp(−z)) in both the numpy and the PyTorch versions. Derive the local derivative in terms of the output (it is σ(z) (1 − σ(z))), update the backward pass, and confirm with the finite-difference check that the worst relative error is still of order 1e-10. Then switch the loss to binary cross-entropy against a target of 1.0 and repeat. Finally, try the gradient check with all tensors in float32 and ε = 1e-5, and explain from the printed numbers why it no longer passes.

Summary

Backpropagation is the chain rule applied once per operation, from the loss back towards the inputs, reusing each intermediate gradient instead of recomputing it.
A linear layer’s weight gradient is the incoming gradient times the layer’s input, and the gradient it passes back is the incoming gradient times its weights. A non-linearity contributes an element-wise factor.
Central finite differences in float64 check an analytic gradient to about ten decimal places and catch sign and transpose errors immediately.
Autograd records each forward operation as a node that stores what its backward formula needs, then walks the graph in reverse and accumulates into .grad.
The numpy and PyTorch gradients agree because they perform the same arithmetic in the same order, so any difference is floating-point rounding.