Neural networks and backpropagation from first principles

A fully connected neural network is a chain of simple functions: multiply by a matrix, add a vector, apply a fixed nonlinearity, repeat. Training it means adjusting every matrix and vector so a loss gets smaller. Backpropagation is the bookkeeping that tells us, for every one of those numbers at once, which direction reduces the loss. It is not a separate learning theory. It is the chain rule from calculus, applied in an order that reuses intermediate results instead of recomputing them.

This guide derives backpropagation for a multilayer perceptron (MLP) carefully enough that you could reproduce it on paper, then checks every claim with numbers. We work one tiny network by hand, verify the gradient with finite differences, and finish with a complete NumPy implementation that trains on a toy dataset. Frameworks such as PyTorch automate exactly this computation; knowing what they automate makes their failures much easier to diagnose.

Neurons and layers are affine maps followed by a nonlinearity

A single neuron receives an input vector $x\in\mathbb R^{n}$, computes a weighted sum plus a bias, and passes the result through an activation function $\sigma$:

The weighted sum $z$ is called the pre-activation. Geometrically, $w^\top x + b = 0$ is a hyperplane, and the neuron reports how far, and on which side, the input lies, squashed or clipped by $\sigma$.

A layer is many neurons reading the same input. Stack their weight vectors as the rows of a matrix and the whole layer becomes one affine map followed by an elementwise nonlinearity. For layer $l$ of a network with $L$ layers:

The nonlinearity is essential. Without it, composing layers gives $W^{2}(W^{1}x + b^{1}) + b^{2}$, which is still a single affine map. Depth adds expressive power only because $\sigma$ bends the space between layers.

Keep shapes visible. Almost every backpropagation bug is a shape or transpose mistake that happens to broadcast without an error. We use column vectors for a single example and stack a batch of $B$ examples as columns.

The forward pass is just the recursion above, evaluated from $l=1$ to $l=L$. The output $a^{L}$ is the prediction, often written $\hat y$.

flowchart LR
  X["Input x<br/>a0"] --> Z1["z1 = W1 a0 + b1"]
  Z1 --> A1["a1 = sigma of z1"]
  A1 --> Z2["z2 = W2 a1 + b2"]
  Z2 --> A2["a2 = sigma of z2"]
  A2 --> ZL["zL = WL a_L-1 + bL"]
  ZL --> P["Output y_hat"]
  P --> LOSS["Loss L of y_hat and y"]

Loss functions turn predictions into one number

Gradient-based training needs a single scalar to decrease. For regression, the usual choice is squared error. Including a factor of one half makes the derivative cleaner:

For classification with $K$ classes, the last layer produces $K$ unconstrained scores, or logits, $z^{L}$. The softmax turns them into a probability vector:

With a one-hot target $y$, the cross-entropy loss is the negative log probability assigned to the correct class:

If the correct class gets probability $0.9$, the loss is about $0.105$. If it gets $0.01$, the loss is about $4.605$. Confident mistakes are punished hard, which is exactly the signal a classifier needs. Over a batch, we average the per-example losses so that the scale of the gradient does not depend on batch size.

Computational graphs and the chain rule

Write the computation as a graph: each node is an intermediate value, and each edge says which values were used to compute which. The forward pass walks the graph from inputs to loss. To differentiate, we walk it backward.

The rule we need is the multivariable chain rule. If a scalar $\mathcal L$ depends on $u$ only through intermediate variables $v_1,\ldots,v_m$, then

In words: the sensitivity of the loss to $u$ is the sum, over every path from $u$ to the loss, of the local derivative along that path times the downstream sensitivity. Each node only needs two things: the gradient arriving from above, and its own local derivative. That locality is why backpropagation is efficient. Computing the gradient for all parameters costs a small constant multiple of one forward pass, instead of one forward pass per parameter as naive finite differences would.

flowchart LR
  subgraph FWD["Forward pass: compute and cache values"]
    direction LR
    F0["a0 = x"] --> F1["z1"] --> F2["a1"] --> F3["z2"] --> F4["a2 = y_hat"] --> F5["L"]
  end
  subgraph BWD["Backward pass: reuse cached values"]
    direction RL
    B5["dL/dL = 1"] --> B4["dL/da2"] --> B3["delta2 = dL/dz2"] --> B2["dL/da1 = W2^T delta2"] --> B1["delta1 = dL/dz1"]
  end
  B3 --> G2["dW2 = delta2 a1^T"]
  B1 --> G1["dW1 = delta1 a0^T"]

The forward pass must store, or cache, the pre-activations $z^{l}$ and activations $a^{l}$. The backward pass reads them. This is why training uses much more memory than inference: inference can discard each activation once the next layer has consumed it.

Deriving backpropagation in matrix form

One layer

Start with a single layer, $z = W a + b$, where $a$ is the layer input, and suppose something downstream has already given us $\delta = \partial \mathcal L / \partial z$. Entry $z_i = \sum_j W_{ij} a_j + b_i$ depends on $W_{ij}$ only through itself, and the local derivative is $a_j$. So

Collecting all $(i,j)$ entries gives an outer product, which has the same shape as $W$, as it must:

To continue backward, we also need the gradient with respect to the layer input. Input $a_j$ influences every $z_i$, so we sum over $i$:

The transpose is not a convention to memorize. It appears because the forward pass sends input $j$ into every output $i$, so the backward pass gathers from every $i$ back into $j$.

L layers: the delta recursion

Define the error signal of layer $l$ as $\delta^{l} = \partial \mathcal L / \partial z^{l}$. At the output layer, apply the chain rule through the final activation:

Here $\odot$ is elementwise multiplication. It appears because $\sigma$ acts elementwise, so its Jacobian is diagonal, and multiplying by a diagonal matrix is the same as multiplying elementwise.

For a hidden layer, $z^{l}$ affects the loss only through $a^{l} = \sigma(z^{l})$, which affects the loss only through $z^{l+1} = W^{l+1} a^{l} + b^{l+1}$. Componentwise,

In matrix form, this is the central equation of backpropagation:

Combined with the one-layer result, the parameter gradients for every layer are

Check the shapes: $(W^{l+1})^\top$ is $n_l \times n_{l+1}$ and $\delta^{l+1}$ is $n_{l+1} \times 1$, so their product is $n_l \times 1$, matching $z^{l}$. And $\delta^{l}(a^{l-1})^\top$ is $n_l \times n_{l-1}$, matching $W^{l}$.

Batch form

With a batch of $B$ examples stacked as columns, $Z^{l} = W^{l} A^{l-1} + b^{l}\mathbf 1^\top$, where the bias is broadcast across columns. The batch loss is the mean of per-example losses, so each parameter gradient is the mean of per-example gradients. Writing $\Delta^{l}$ for the matrix whose columns are the per-example deltas:

The matrix product $\Delta^{l}(A^{l-1})^\top$ is precisely the sum of the $B$ outer products, so a single matrix multiplication does the averaging work. The bias gradient sums $\Delta^{l}$ across columns, because the bias was copied to every column in the forward pass. In code, it is convenient to divide by $B$ once, at the output delta, and then let the recursion carry that factor through.

The softmax and cross-entropy gradient is p minus y

For a classifier, we could compute $\partial \mathcal L/\partial p$ and then multiply by the softmax Jacobian. That works but is numerically fragile, since $\partial \mathcal L/\partial p_k = -y_k / p_k$ explodes when $p_k$ is tiny. Combining the two steps gives a much simpler result.

First, the softmax Jacobian. Differentiate $p_k = \exp(z_k)/S$ with $S = \sum_j \exp(z_j)$, using the quotient rule. If $j = k$:

If $j \neq k$, only the denominator depends on $z_j$:

Both cases fit one formula using the Kronecker delta $\mathbb 1[j = k]$:

Now apply the chain rule to $\mathcal L = -\sum_k y_k \log p_k$:

The last step uses $\sum_k y_k = 1$, which holds for one-hot targets and for any target probability vector. So

The intuition is pleasant: the error signal is the gap between what the network predicted and what was true. If the network assigns probability $0.7$ to the correct class, that logit receives gradient $-0.3$ and is pushed up; every wrong class receives its own probability as gradient and is pushed down. The same simplification happens for a sigmoid output with binary cross-entropy, and for a linear output with squared error: each is the canonical pairing of an output function with its loss.

Activation derivatives, saturation, and initialization

The delta recursion multiplies by $\sigma'(z^{l})$ at every layer. The size of that factor therefore controls whether gradients survive the trip back through a deep network. The three classic activations and their derivatives are:

The sigmoid derivative peaks at $z = 0$ with value $1/4$ and decays quickly. When a neuron's pre-activation is large in magnitude, the neuron is saturated: its output barely changes, so its derivative is nearly zero, and almost no gradient passes through it. Even in the best case, a chain of sigmoid layers multiplies the signal by at most $0.25$ per layer, ignoring the weights.

import numpy as np

sigmoid = lambda z: 1 / (1 + np.exp(-z))
for z in [0.0, 2.0, 5.0, 10.0]:
    s = sigmoid(z)
    print(f"z = {z:5.1f}   sigmoid'(z) = {s * (1 - s):.6f}")

# Best case: every layer sits at z = 0, so each factor is 0.25.
for depth in [1, 5, 10, 20]:
    print(f"{depth:2d} layers: 0.25 ** depth = {0.25 ** depth:.3e}")

Running it shows $0.25$ at $z = 0$, about $0.105$ at $z = 2$, about $0.0066$ at $z = 5$, and about $4.5\times 10^{-5}$ at $z = 10$. Ten layers at the most favourable point already shrink the factor to about $9.5\times 10^{-7}$. Larger weights can partly compensate, but they also push pre-activations into saturation. This is the vanishing-gradient problem in its simplest form. Tanh is zero-centred and has a maximum derivative of $1$, which helps, but it saturates too. ReLU passes gradient with factor exactly $1$ wherever it is active, which is a large part of why it made deeper networks practical. Its failure mode is different: a unit whose pre-activation is negative for every input receives zero gradient forever, a so-called dead ReLU.

Why initialization scale matters

Initialization decides where the pre-activations start. Consider one unit, $z_i = \sum_{j=1}^{n_{\mathrm{in}}} W_{ij} a_j$, with weights drawn independently with mean zero and variance $\operatorname{Var}(W)$, independent of zero-mean inputs. Then

If we want the signal's scale to stay roughly constant from layer to layer, we need $n_{\mathrm{in}}\operatorname{Var}(W) \approx 1$ for an activation that behaves like the identity near zero, such as tanh. Applying the same argument to the backward pass gives $n_{\mathrm{out}}\operatorname{Var}(W) \approx 1$. Glorot and Bengio (2010) compromise between the two:

ReLU zeroes roughly half of a symmetric input distribution, so $\mathbb E[a^{2}] = \operatorname{Var}(z)/2$ for the previous layer's pre-activation. To compensate, He et al. (2015) double the variance:

These are heuristics derived under independence assumptions, not guarantees. They matter most in deep networks; in a two-layer toy problem, almost any reasonable small scale trains. One rule is not optional: do not initialize all weights to the same value. Every unit in a layer would then compute the same function and receive the same gradient, and they would never become different.

A fully worked numeric example

Here is one forward pass, one backward pass, and one stochastic gradient descent (SGD) step for a 2-2-1 network, small enough to follow with a calculator. The hidden layer uses a sigmoid; the output is linear with loss $\frac12(\hat y - y)^{2}$, so $\delta^{2} = \hat y - y$.

Forward

Backward

The output delta and the output-layer gradients:

Send the error back through $W^{2}$, then through the sigmoid, whose derivative is $a^{1}\odot(1 - a^{1})$:

Notice how much smaller $\delta^{1}$ is than $\delta^{2}$: the sigmoid derivative of roughly $0.23$ cut it by a factor of about four, as the previous section predicted.

One SGD step

With learning rate $\eta = 0.1$, each parameter moves by $-\eta$ times its gradient, $\theta \leftarrow \theta - \eta\,\partial\mathcal L/\partial\theta$:

Re-running the forward pass with the updated parameters gives a loss of about $0.245718$, down from $0.383832$. The code below reproduces every number, and also checks one weight gradient with a centred finite difference; the analytic and numeric values agree to about ten significant figures.

import numpy as np

sigmoid = lambda z: 1 / (1 + np.exp(-z))
x = np.array([[1.0], [2.0]])
y = 1.0
W1 = np.array([[0.1, 0.2], [-0.3, 0.4]])
b1 = np.array([[0.0], [0.1]])
W2 = np.array([[0.5, -0.6]])
b2 = np.array([[0.2]])
lr = 0.1


def loss(W1, b1, W2, b2):
    a1 = sigmoid(W1 @ x + b1)
    return 0.5 * ((W2 @ a1 + b2).item() - y) ** 2


# Forward
z1 = W1 @ x + b1
a1 = sigmoid(z1)
y_hat = (W2 @ a1 + b2).item()
print("z1", z1.ravel(), "a1", a1.ravel().round(6))
print("y_hat", round(y_hat, 6), "loss", round(loss(W1, b1, W2, b2), 6))

# Backward
d2 = np.array([[y_hat - y]])
dW2, db2 = d2 @ a1.T, d2
d1 = (W2.T @ d2) * a1 * (1 - a1)
dW1, db1 = d1 @ x.T, d1
print("dW2", dW2.round(6), "db2", db2.ravel().round(6))
print("delta1", d1.ravel().round(6))
print("dW1", dW1.round(6))

# One SGD step
W1n, b1n, W2n, b2n = W1 - lr * dW1, b1 - lr * db1, W2 - lr * dW2, b2 - lr * db2
print("new loss", round(loss(W1n, b1n, W2n, b2n), 6))

# Finite-difference check of dL/dW1[0, 1]
eps = 1e-6
Wp, Wm = W1.copy(), W1.copy()
Wp[0, 1] += eps
Wm[0, 1] -= eps
numeric = (loss(Wp, b1, W2, b2) - loss(Wm, b1, W2, b2)) / (2 * eps)
print(f"numeric {numeric:.10f}  analytic {dW1[0, 1]:.10f}")

Gradient checking with finite differences

A backpropagation implementation can be subtly wrong and still reduce the loss, because a gradient that is merely correlated with the true one can point downhill. The only way to know is to compare against an independent estimate. The centred difference

has truncation error proportional to $\epsilon^{2}$, compared with $\epsilon$ for the one-sided version. Making $\epsilon$ too small introduces floating-point cancellation, so in double precision a value around $10^{-5}$ to $10^{-6}$ is a sensible default. Compare with a relative error, because gradients vary widely in scale:

As a rough guide in float64, values around $10^{-7}$ or smaller indicate a correct implementation, while values near $10^{-2}$ or larger almost always mean a bug. Some practical rules make the check trustworthy:

• Use float64 for the check, even if training later uses lower precision.
• Turn off randomness such as dropout, and keep the batch fixed between the two evaluations.
• Check a small batch and a random sample of parameters in every layer, including biases. A bug confined to one layer can hide in an average.
• ReLU has a kink at zero. If a perturbation crosses it, the finite difference is meaningless for that coordinate. Occasional isolated failures with ReLU are worth inspecting before panicking.
• Gradient checking is slow, roughly two forward passes per parameter checked. Use it to validate the code, then switch it off.

A complete MLP in NumPy

The program below generates a two-moons dataset, which is not linearly separable, builds a 2-16-16-2 network with ReLU hidden layers and a softmax output, checks its gradients, and trains it with minibatch SGD. Every gradient follows the equations above: the output delta is $(P - Y)/B$, the hidden deltas follow the recursion with the ReLU derivative, and the weight gradient is $\Delta (A^{l-1})^\top$. Examples are columns, so the code and the math use the same shapes.

import numpy as np

rng = np.random.default_rng(0)


def make_moons(n, noise=0.15):
    """Two interleaving half-circles. Returns X of shape (2, n) and labels y of shape (n,)."""
    m = n // 2
    t1 = rng.uniform(0, np.pi, m)
    t2 = rng.uniform(0, np.pi, n - m)
    upper = np.stack([np.cos(t1), np.sin(t1)])
    lower = np.stack([1 - np.cos(t2), 0.5 - np.sin(t2)])
    X = np.concatenate([upper, lower], axis=1) + noise * rng.standard_normal((2, n))
    y = np.concatenate([np.zeros(m, dtype=int), np.ones(n - m, dtype=int)])
    return X, y


def init_params(sizes):
    """He initialization for ReLU layers: W ~ N(0, 2 / fan_in)."""
    params = []
    for n_in, n_out in zip(sizes[:-1], sizes[1:]):
        W = rng.standard_normal((n_out, n_in)) * np.sqrt(2.0 / n_in)
        b = np.zeros((n_out, 1))
        params.append([W, b])
    return params


def relu(z):
    return np.maximum(z, 0.0)


def softmax(z):
    z = z - z.max(axis=0, keepdims=True)  # stability: subtract the column max
    e = np.exp(z)
    return e / e.sum(axis=0, keepdims=True)


def forward(params, X):
    """Columns are examples. Returns the class probabilities and a cache for backprop."""
    a = X
    cache = [(None, a)]
    for i, (W, b) in enumerate(params):
        z = W @ a + b
        a = softmax(z) if i == len(params) - 1 else relu(z)
        cache.append((z, a))
    return a, cache


def loss_and_grads(params, X, y):
    B = X.shape[1]
    P, cache = forward(params, X)
    Y = np.eye(P.shape[0])[:, y]  # one-hot targets, shape (K, B)
    loss = -np.mean(np.log(P[y, np.arange(B)] + 1e-12))

    grads = [None] * len(params)
    delta = (P - Y) / B  # dL/dz at the output layer, already averaged over the batch
    for l in range(len(params) - 1, -1, -1):
        W, _ = params[l]
        a_prev = cache[l][1]
        grads[l] = [delta @ a_prev.T, delta.sum(axis=1, keepdims=True)]
        if l > 0:
            z_prev = cache[l][0]
            delta = (W.T @ delta) * (z_prev > 0)  # ReLU derivative
    return loss, grads


def grad_check(params, X, y, eps=1e-5, checks=8):
    _, grads = loss_and_grads(params, X, y)
    worst = 0.0
    for l, (W, b) in enumerate(params):
        for P, G in ((W, grads[l][0]), (b, grads[l][1])):
            for _ in range(checks):
                idx = tuple(rng.integers(s) for s in P.shape)
                old = P[idx]
                P[idx] = old + eps
                lp, _ = loss_and_grads(params, X, y)
                P[idx] = old - eps
                lm, _ = loss_and_grads(params, X, y)
                P[idx] = old
                numeric = (lp - lm) / (2 * eps)
                rel = abs(numeric - G[idx]) / max(1e-8, abs(numeric) + abs(G[idx]))
                worst = max(worst, rel)
    return worst


X, y = make_moons(400)
X_test, y_test = make_moons(200)
params = init_params([2, 16, 16, 2])

print(f"worst relative error in gradient check: {grad_check(params, X[:, :20], y[:20]):.2e}")
print(f"loss before training: {loss_and_grads(params, X, y)[0]:.4f}")

lr, batch = 0.5, 32
for epoch in range(301):
    order = rng.permutation(X.shape[1])
    for start in range(0, X.shape[1], batch):
        idx = order[start:start + batch]
        _, grads = loss_and_grads(params, X[:, idx], y[idx])
        for (W, b), (dW, db) in zip(params, grads):
            W -= lr * dW
            b -= lr * db
    if epoch % 100 == 0:
        loss, _ = loss_and_grads(params, X, y)
        acc = np.mean(forward(params, X)[0].argmax(axis=0) == y)
        print(f"epoch {epoch:3d}  loss {loss:.4f}  train accuracy {acc:.3f}")

test_acc = np.mean(forward(params, X_test)[0].argmax(axis=0) == y_test)
print(f"test accuracy: {test_acc:.3f}")

In my run with this seed, the worst relative error in the gradient check was about $6.6\times 10^{-10}$. The training loss started at about $1.06$, fell to about $0.28$ after the first epoch, and reached about $0.010$ after 300 epochs. Training accuracy ended near $99.8$ percent and accuracy on the 200 held-out points was $97.5$ percent. Your numbers will vary with the seed, the noise level, and the NumPy version. The point is not the score on a toy problem; it is that the gradient check passed before any training happened, so a falling loss can be trusted to mean the optimizer is following the true gradient.

A few design choices deserve a note. The softmax subtracts the column maximum before exponentiating; this leaves the probabilities unchanged but prevents overflow. The loss adds a tiny constant inside the logarithm only as a guard; the gradient never uses that logarithm, because it uses $p - y$ directly. And the update modifies the arrays in place with -=, so the list params sees the new values without being rebuilt.

flowchart TD
  S["Sample a minibatch"] --> F["Forward pass<br/>cache z and a for every layer"]
  F --> L["Compute mean loss"]
  L --> D["Output delta = P - Y over B"]
  D --> R["Backward recursion<br/>delta = W^T delta times relu prime"]
  R --> G["Gradients dW = delta a_prev^T and db = row sums"]
  G --> U["SGD update W -= lr dW"]
  U --> S

Common bugs and practice exercises

Most backpropagation bugs fall into a handful of families. Knowing them turns a mysterious plateau into a short checklist.

• Missing or extra transpose. Using $W^{l+1}\delta^{l+1}$ instead of $(W^{l+1})^\top\delta^{l+1}$ fails loudly when layers have different widths, but silently runs when they are square. Test with unequal layer sizes.
• Wrong cached value. The derivative must be evaluated at the pre-activation $z^{l}$ of the same forward pass. Writing $\sigma'(a^{l})$ instead of $\sigma'(z^{l})$ is a classic sigmoid mistake; for sigmoid, the shortcut $a^{l}(1 - a^{l})$ is correct precisely because it is written in terms of $a^{l}$.
• Batch averaging twice or not at all. If the loss is a mean, divide by $B$ exactly once. Dividing at the output delta and again in the weight gradient makes the effective learning rate shrink with batch size.
• Bias gradient summed over the wrong axis. With examples as columns, sum across columns. Broadcasting will happily produce an array of the wrong shape that still subtracts from $b$ without an error.
• Applying softmax twice, or softmax and then a loss that expects logits. Framework losses such as cross-entropy on logits include the softmax internally.
• Updating weights before finishing the backward pass. The recursion for layer $l$ needs the old $W^{l+1}$. Compute all gradients first, then update.
• Symmetric or badly scaled initialization. All-zero hidden weights never break symmetry. Very large weights saturate sigmoids or tanh units from the first step.
• Overflow in exp or log. Subtract the maximum logit before exponentiating and never take the log of a probability that can be exactly zero.

The habit that catches nearly all of these: run a gradient check on a tiny batch, then confirm that the network can drive the loss almost to zero on a handful of examples. If it cannot overfit eight points, something is wrong with the gradients, the data pipeline, or the learning rate, and adding data or depth will not fix it.

Practice

1. Redo the worked example by hand with a tanh hidden layer instead of a sigmoid. Before computing, predict whether $\delta^{1}$ will be larger or smaller, using the maximum derivative of each activation.
2. Show that for a sigmoid output $p = \operatorname{sigmoid}(z)$ with binary cross-entropy $-y\log p - (1-y)\log(1-p)$, the gradient with respect to $z$ is $p - y$.
3. Introduce a deliberate bug into the NumPy MLP, such as removing the transpose or dividing by $B$ twice, and see how the gradient-check error changes. Which bugs does training alone fail to reveal?
4. Replace He initialization with weights multiplied by $0.01$ and then by $3$ in a network with eight hidden layers. Print the standard deviation of each layer's activations at initialization and explain what you see using the variance argument.
5. Train the MLP on XOR, generated as the four corners of the unit square, with a single hidden layer of two units. Run several seeds. How often does it fail, and what do the hidden units look like when it does?

Further reading

The method was popularized for neural networks by Rumelhart, Hinton and Williams (1986). The Deep Learning book by Goodfellow, Bengio and Courville treats computational graphs and backpropagation in its sixth chapter, and the initialization papers by Glorot and Bengio and He et al. are short and readable once the variance argument above is familiar.