Backpropagation and Gradient Descent

Imagine standing on a hillside in thick fog, trying to reach the valley floor. You cannot see the valley, but you can feel the slope under your feet, so you take a small step downhill, feel again, and repeat. That is gradient descent, and it is how every modern neural network is trained. The difficulty is that a network’s “hillside” has not two directions but millions or billions, one for every weight, and you need the slope in all of them at once, thousands of times over. Backpropagation is the procedure that delivers all those slopes for roughly the cost of running the network twice. It is one of the most consequential algorithms of the past half-century, and it fits in about ten lines of NumPy.

Loss: the number we make smaller

Training needs a single number that says how wrong the network is: the loss. Everything else follows from choosing it. Two choices cover most of what you will meet.

The maths page explained where cross-entropy comes from (maximum likelihood). Here is why it is preferred for classification. Suppose the true class is 1 and the model’s probability for it slides from good to terrible:

import numpy as np

y_true = 1.0
for y_pred in [0.9, 0.5, 0.1, 0.01]:
    mse = (y_pred - y_true) ** 2
    ce = -np.log(y_pred)             # cross-entropy when the true class is 1
    print(f"predicted {y_pred:<5} squared error {mse:.3f}   cross-entropy {ce:.3f}")

predicted 0.9   squared error 0.010   cross-entropy 0.105
predicted 0.5   squared error 0.250   cross-entropy 0.693
predicted 0.1   squared error 0.810   cross-entropy 2.303
predicted 0.01  squared error 0.980   cross-entropy 4.605

Squared error can never exceed 1 here, so a confidently wrong answer (0.01) is barely punished more than a merely wrong one (0.1). Cross-entropy keeps climbing — towards infinity as the probability approaches zero — so confident mistakes produce large losses and, as we will see, large corrective gradients. The loss is a design decision: it is the only thing the network is told about what you want.

Gradient descent and the learning rate

Collect every weight and bias of a network into one long vector θ (theta). The loss is a function of θ, and its gradient ∇L(θ) points in the direction in which the loss rises fastest. Gradient descent steps the other way:

The method goes back to Augustin-Louis Cauchy in 1847. Its one dial is the learning rate, and the simplest possible “model” shows why it matters. Take a single weight w with loss w² (minimum at 0, gradient 2w), start at 1.5, and take ten steps:

def loss(w):                      # a one-weight "model" with its minimum at w = 0
    return w ** 2

def grad(w):
    return 2 * w

for lr in [0.1, 0.5, 0.9, 1.1]:
    w = 1.5
    for step in range(10):
        w = w - lr * grad(w)
    print(f"learning rate {lr}: w after 10 steps = {w:10.4f}   loss = {loss(w):.4g}")

learning rate 0.1: w after 10 steps =     0.1611   loss = 0.02594
learning rate 0.5: w after 10 steps =     0.0000   loss = 0
learning rate 0.9: w after 10 steps =     0.1611   loss = 0.02594
learning rate 1.1: w after 10 steps =     9.2876   loss = 86.26

Each step multiplies w by (1 − 2η). With η = 0.1 that factor is 0.8: steady, cautious progress. With η = 0.5 it is 0, so this particular bowl is solved in one step — a coincidence of its shape, not something real networks offer. With η = 0.9 the factor is −0.8: every step overshoots to the other side of the valley, and the weight zig-zags inward, ending exactly as far along as the cautious run. With η = 1.1 the factor is −1.2 and each overshoot is bigger than the last: the loss explodes from 2.25 to 86. Try it yourself:

Push the learning rate up until the ball starts bouncing from side to side, then further until it flies off. Real training curves show the same three regimes: crawling, zig-zagging and diverging.

Real loss surfaces are not neat bowls. They have long narrow valleys, flat plateaus and saddle points, and the best learning rate differs between directions. That is why practical training uses refinements — momentum, Adam, learning-rate schedules — which Training Well takes up. All of them still need the same raw ingredient: the gradient.

The chain rule through a computation graph

To get the gradient efficiently, think of the forward pass as a computation graph: each node is one simple operation, and the loss sits at the end. Here is the smallest interesting example — one input x = 2, one weight w = 0.5, one bias b = 0.1, a tanh neuron and squared error against a target y = 1. Tap along the graph: the forward pass runs left to right, then the gradients flow back right to left.

The same computation as code, with a numerical check at the end:

import numpy as np

x, y = 2.0, 1.0                  # one training example
w, b = 0.5, 0.1                  # the parameters

# forward: record every intermediate value
z = w * x + b                    # 1.1
a = np.tanh(z)                   # squashed prediction
L = (a - y) ** 2                 # squared error

# backward: local slope times the gradient arriving from above
dL_da = 2 * (a - y)
dL_dz = dL_da * (1 - a ** 2)     # tanh'(z) = 1 - tanh(z)^2
dL_dw = dL_dz * x
dL_db = dL_dz * 1.0

print(f"z={z:.4f}  a={a:.4f}  L={L:.4f}")
print(f"dL/da={dL_da:.4f}  dL/dz={dL_dz:.4f}  dL/dw={dL_dw:.4f}  dL/db={dL_db:.4f}")

h = 1e-6                         # check dL/dw by nudging w
numeric = ((np.tanh((w + h) * x + b) - y) ** 2 - (np.tanh((w - h) * x + b) - y) ** 2) / (2 * h)
print(f"numeric dL/dw={numeric:.4f}")

z=1.1000  a=0.8005  L=0.0398
dL/da=-0.3990  dL/dz=-0.1433  dL/dw=-0.2866  dL/db=-0.1433
numeric dL/dw=-0.2866

Two things make this backpropagation rather than just “using the chain rule”. First, the backward pass reuses work: ∂L/∂z is computed once and shared by both w and b. In a big network, the gradient arriving at a layer is shared by every weight in it, so one backward sweep gives all the gradients for about the cost of one more forward pass. Nudging each weight separately, as the numeric check does, would need two forward passes per weight. Second, the forward pass must store its intermediate values (here a), because the backward pass needs them. That stored state is why training uses far more memory than running a model.

The idea has a long, multi-national history. Seppo Linnainmaa described reverse-mode automatic differentiation — the general algorithm — in a Finnish master’s thesis in 1970; Paul Werbos proposed applying it to neural networks in his 1974 Harvard thesis; and David Rumelhart, Geoffrey Hinton and Ronald Williams’ 1986 Nature paper showed convincingly that it lets hidden layers learn useful internal representations, launching the revival told in Backprop and the Connectionist Revival.

Backprop for the two-layer network

Now the real thing. Our TwoLayerNet from Neural Networks from Scratch computes, for a batch: H = tanh(XW₁ + b₁), then P = softmax(HW₂ + b₂), then the mean cross-entropy. Walking backwards through that graph gives five formulas.

A trick that makes these formulas easy to remember: a gradient always has the same shape as the thing it is the gradient of. ∂L/∂W₂ must be (16, 2) like W₂. The ingredients are H, which is (200, 16), and ∂L/∂Z₂, which is (200, 2). The only way to combine them into (16, 2) is Hᵀ @ ∂L/∂Z₂. Shapes almost write backprop for you. Here is the class again, now with a loss and a backward method. Save it as net.py; the next three snippets import it.

# net.py
import numpy as np

def softmax(Z):
    Z = Z - Z.max(axis=1, keepdims=True)
    E = np.exp(Z)
    return E / E.sum(axis=1, keepdims=True)

class TwoLayerNet:
    def __init__(self, n_in, n_hidden, n_out, seed=42):
        rng = np.random.default_rng(seed)
        self.W1 = rng.normal(0, 1 / np.sqrt(n_in), size=(n_in, n_hidden))
        self.b1 = np.zeros(n_hidden)
        self.W2 = rng.normal(0, 1 / np.sqrt(n_hidden), size=(n_hidden, n_out))
        self.b2 = np.zeros(n_out)

    def forward(self, X):
        self.X = X
        self.H = np.tanh(X @ self.W1 + self.b1)
        self.P = softmax(self.H @ self.W2 + self.b2)
        return self.P

    def loss(self, X, y):                      # mean cross-entropy
        P = self.forward(X)
        return -np.log(P[np.arange(len(y)), y]).mean()

    def backward(self, y):                     # call after forward/loss
        n = len(y)
        dZ2 = self.P.copy()
        dZ2[np.arange(n), y] -= 1              # softmax + cross-entropy: P - onehot
        dZ2 /= n                               # because the loss is a mean
        self.dW2 = self.H.T @ dZ2              # (hidden, out)
        self.db2 = dZ2.sum(axis=0)
        dH = dZ2 @ self.W2.T                   # send the blame back a layer
        dZ1 = dH * (1 - self.H ** 2)           # through tanh
        self.dW1 = self.X.T @ dZ1              # (in, hidden)
        self.db1 = dZ1.sum(axis=0)

    def params(self):
        return [("W1", self.W1, "dW1"), ("b1", self.b1, "db1"),
                ("W2", self.W2, "dW2"), ("b2", self.b2, "db2")]

The backward method is the five formulas, line for line. Subtracting 1 at each example’s true class turns P into P − Y without building the one-hot matrix. params pairs each parameter with the name of its gradient, so a training loop can update them all uniformly.

Never trust a backward pass you have not checked. Hand-written gradients are easy to get subtly wrong — a missing transpose, a forgotten 1/n — and a wrong gradient often still trains, just badly, so the bug hides. The cure is a gradient check: compare every analytic gradient with a central-difference estimate.

import numpy as np
from sklearn.datasets import make_moons
from net import TwoLayerNet

X, y = make_moons(n_samples=200, noise=0.2, random_state=42)
net = TwoLayerNet(2, 16, 2)
net.loss(X, y)
net.backward(y)

h = 1e-5
for name, param, gname in net.params():
    analytic = getattr(net, gname)
    numeric = np.zeros_like(param)
    for i in np.ndindex(param.shape):          # nudge every single number
        old = param[i]
        param[i] = old + h; up = net.loss(X, y)
        param[i] = old - h; down = net.loss(X, y)
        param[i] = old
        numeric[i] = (up - down) / (2 * h)
    rel = np.abs(analytic - numeric).max() / (np.abs(analytic).max() + np.abs(numeric).max())
    print(f"{name}: max relative difference {rel:.2e}")

W1: max relative difference 7.74e-11
b1: max relative difference 1.60e-10
W2: max relative difference 3.83e-11
b2: max relative difference 2.14e-11

Relative differences around 10⁻¹⁰ mean the analytic and numeric gradients agree to about ten significant figures: the backward pass is right. As a rule of thumb, values below about 10⁻⁷ are reassuring; values near 10⁻² or larger mean a bug. Gradient checks are slow (two forward passes for each of the 82 parameters here), so you run them once on a small network to validate the code, then switch them off.

Training: full batches, mini-batches and SGD

With a correct gradient, training is a loop of four steps, repeated. One pass through the whole training set is called an epoch.

First, the simplest version: full-batch gradient descent, using all 200 two-moons points for every step.

import numpy as np
from sklearn.datasets import make_moons
from net import TwoLayerNet

X, y = make_moons(n_samples=200, noise=0.2, random_state=42)
net = TwoLayerNet(2, 16, 2)
lr = 0.5

for epoch in range(2001):
    loss = net.loss(X, y)                      # forward
    net.backward(y)                            # backward
    for _, param, gname in net.params():       # gradient descent step
        param -= lr * getattr(net, gname)
    if epoch % 400 == 0:
        acc = (net.P.argmax(axis=1) == y).mean()
        print(f"epoch {epoch:>4}  loss {loss:.4f}  accuracy {acc:.3f}")

epoch    0  loss 0.7100  accuracy 0.550
epoch  400  loss 0.1096  accuracy 0.965
epoch  800  loss 0.0775  accuracy 0.970
epoch 1200  loss 0.0703  accuracy 0.970
epoch 1600  loss 0.0658  accuracy 0.970
epoch 2000  loss 0.0609  accuracy 0.970

The loss starts at 0.71 — the untrained value from the previous page — and falls fast, then slowly. Accuracy jumps from a coin flip to 96.5% in the first 400 steps, then creeps to 97%. The loss keeps falling after accuracy stops improving: the network is growing more confident on points it already gets right. Note that param -= … updates the array in place, which is why the network’s own attributes change. These are training-set figures, though, and the Builder track taught you not to trust those alone.

Full-batch descent has a scaling problem: every single step reads the entire dataset. With millions of examples that is absurdly slow. The answer, which dates back to Herbert Robbins and Sutton Monro’s stochastic approximation (1951), is to estimate the gradient from a small random mini-batch of examples. Each estimate is noisy, but it is cheap, so you can take many more steps for the same compute. This is stochastic gradient descent (SGD) — strictly, SGD uses one example per step, but in practice “SGD” almost always means mini-batches.

import numpy as np
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from net import TwoLayerNet

X, y = make_moons(n_samples=1000, noise=0.2, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)

net = TwoLayerNet(2, 16, 2)
rng = np.random.default_rng(42)
lr, batch_size = 0.5, 32

for epoch in range(1, 31):
    order = rng.permutation(len(X_train))          # reshuffle every epoch
    for start in range(0, len(order), batch_size):
        idx = order[start:start + batch_size]      # one mini-batch
        net.loss(X_train[idx], y_train[idx])
        net.backward(y_train[idx])
        for _, param, gname in net.params():
            param -= lr * getattr(net, gname)
    if epoch in (1, 5, 10, 30):
        test_acc = (net.forward(X_test).argmax(axis=1) == y_test).mean()
        print(f"epoch {epoch:>2}  train loss {net.loss(X_train, y_train):.4f}  test accuracy {test_acc:.3f}")

epoch  1  train loss 0.3190  test accuracy 0.884
epoch  5  train loss 0.2776  test accuracy 0.904
epoch 10  train loss 0.2023  test accuracy 0.948
epoch 30  train loss 0.1054  test accuracy 0.980

Now there are 750 training points and batches of 32, so each epoch makes 24 updates instead of one. After a single epoch — one look at each training point — the training loss is already 0.32; full-batch descent needed about 35 complete passes over its data to reach the same loss. After 30 epochs (720 updates) the network classifies 98% of the 250 held-out test points correctly. Reshuffling every epoch matters: it stops the batches arriving in the same order and gives each step a fresh, unbiased sample. Batch size is a trade-off: larger batches give smoother gradients and use hardware more efficiently; smaller ones give more updates per epoch, and their noise is often thought to help find solutions that generalise.

When gradients vanish or explode

Look again at the backward pass: at each layer, the incoming gradient is multiplied by a weight matrix and by the activation’s slope. Through many layers, that is a long product. If the factors are mostly smaller than 1, the product shrinks towards zero — a vanishing gradient, and the early layers barely learn. If they are mostly larger than 1, it grows without limit — an exploding gradient, and training blows up into nan. The sigmoid is especially prone to vanishing, because its slope is never more than 0.25. Here is a 20-layer network with identical weight scales, once with sigmoids and once with ReLUs, measuring how large the weight gradient is at each layer:

import numpy as np

def gradient_sizes(activation, depth=20, width=32, seed=42):
    rng = np.random.default_rng(seed)
    x = rng.normal(size=(64, width))
    Ws, acts = [], [x]
    for _ in range(depth):                          # forward through many layers
        W = rng.normal(0, np.sqrt(1 / width), size=(width, width))
        z = acts[-1] @ W
        a = 1 / (1 + np.exp(-z)) if activation == "sigmoid" else np.maximum(0, z)
        Ws.append(W); acts.append(a)
    grad = np.ones_like(acts[-1])                   # pretend dLoss/dOutput = 1
    sizes = []
    for layer in range(depth, 0, -1):               # backward, layer by layer
        a = acts[layer]
        local = a * (1 - a) if activation == "sigmoid" else (a > 0)
        grad_z = grad * local
        sizes.append(np.abs(acts[layer - 1].T @ grad_z).mean())   # size of dLoss/dW
        grad = grad_z @ Ws[layer - 1].T
    return sizes[::-1]                              # layer 1 first

for act in ["sigmoid", "relu"]:
    s = gradient_sizes(act)
    print(f"{act:>7}: layer 20 {s[19]:.2e}  layer 10 {s[9]:.2e}  layer 1 {s[0]:.2e}")

sigmoid: layer 20 7.17e+00  layer 10 4.37e-06  layer 1 9.72e-13
   relu: layer 20 3.32e-02  layer 10 2.94e-02  layer 1 2.38e-02

In the sigmoid network, the gradient at layer 1 is about 10¹³ times smaller than at layer 20: the first layers effectively cannot learn. In the ReLU network the gradient stays the same order of magnitude all the way down, because the ReLU’s slope is exactly 1 wherever a unit is active. (Ignore the absolute values; they depend on the batch size and on our pretend loss. The ratio between layers is the point.)

This problem was identified by Sepp Hochreiter in his 1991 diploma thesis and analysed by Yoshua Bengio, Patrice Simard and Paolo Frasconi in 1994, and it is a big part of why deep networks were so hard to train before about 2010. Xavier Glorot and Yoshua Bengio’s 2010 study showed how sigmoid units saturate in deep networks and proposed scaling the initial weights to keep signals a steady size in both directions. Consensus Today’s standard remedies, roughly in the order they appeared, are: gated recurrent units such as the LSTM; careful initialisation; ReLU-family activations; normalisation layers; residual (skip) connections, which give the gradient a direct path around each block; and, for explosions, gradient clipping — capping the gradient’s length before each update. The page on training well covers them properly.

Frameworks do this for you

You will rarely write a backward pass by hand again. Frameworks such as PyTorch, JAX and TensorFlow implement automatic differentiation: as your forward code runs, they record the computation graph, and one call walks it backwards applying each operation’s local derivative — exactly the procedure from the chain above. Here is the one-neuron example from earlier, in PyTorch (version 2.14 was used As of 2026):

import torch

x, y = torch.tensor(2.0), torch.tensor(1.0)
w = torch.tensor(0.5, requires_grad=True)   # ask torch to track these
b = torch.tensor(0.1, requires_grad=True)

L = (torch.tanh(w * x + b) - y) ** 2        # forward: torch records the graph
L.backward()                                # backward: chain rule, automatically
print(f"L={L.item():.4f}  dL/dw={w.grad.item():.4f}  dL/db={b.grad.item():.4f}")

L=0.0398  dL/dw=-0.2866  dL/db=-0.1433

The same three numbers we computed by hand. requires_grad=True tells PyTorch to track every operation involving those tensors; L.backward() runs reverse-mode differentiation and leaves each gradient in .grad. The same two lines work unchanged for a model with billions of parameters.

So why learn the manual version? Because when training goes wrong — a loss that will not fall, gradients that turn into nan, a layer that never changes — you need to know what the framework is computing to find out why. Andrej Karpathy’s micrograd, a working autograd engine in about a hundred lines of Python, is the best next step: it shows how the recording and the backward walk are actually implemented.

Key terms on this page

Sources & further reading

• David E. Rumelhart, Geoffrey E. Hinton & Ronald J. Williams, “Learning representations by back-propagating errors”, Nature 323, 533–536 (1986).
• Xavier Glorot & Yoshua Bengio, “Understanding the difficulty of training deep feedforward neural networks”, AISTATS (2010).
• Yoshua Bengio, Patrice Simard & Paolo Frasconi, “Learning long-term dependencies with gradient descent is difficult”, IEEE Transactions on Neural Networks 5 (1994); Herbert Robbins & Sutton Monro, “A stochastic approximation method”, Annals of Mathematical Statistics 22 (1951).
• Ian Goodfellow, Yoshua Bengio & Aaron Courville, Deep Learning (MIT Press, 2016) — section 6.5 on back-propagation and chapter 8 on optimisation.
• Andrej Karpathy, micrograd (2020) and the lecture “The spelled-out intro to neural networks and backpropagation: building micrograd”, part of Neural Networks: Zero to Hero (2022).
• Grant Sanderson (3Blue1Brown), Neural Networks (video series, 2017 onwards) — “Gradient descent, how neural networks learn”, “What is backpropagation really doing?” and “Backpropagation calculus”.
• Adam Paszke et al., “PyTorch: An imperative style, high-performance deep learning library”, NeurIPS (2019).