Backpropagation — How Neural Networks Learn

A manager wants to know why a project was delivered late. They start at the final delay (the loss) and trace backwards. The deployment was late because testing was late. Testing was late because development was late. Development was late because requirements were unclear. At each step they answer: how much did this step contribute to the final delay? That is the chain rule — each step's contribution multiplied together to reach the root cause.

Backpropagation traces the prediction error backwards through the network. At each layer it asks: how much did this layer's weights contribute to the final error? The answer — the gradient — tells each weight exactly how to change to reduce the error next time.

import numpy as np

# ── Chain rule on a simple example before the full network ────────────
# f(x) = sigmoid(2x + 1)  — two operations chained
# We want df/dx

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

def sigmoid_deriv(z):
    s = sigmoid(z)
    return s * (1 - s)

x = 0.5

# Forward pass — compute and cache intermediate values
z = 2 * x + 1        # linear step:   z = 2x + 1
f = sigmoid(z)        # sigmoid step:  f = sigmoid(z)

# Chain rule: df/dx = (df/dz) × (dz/dx)
df_dz = sigmoid_deriv(z)   # how much f changes with z
dz_dx = 2.0                 # how much z changes with x  (derivative of 2x+1 is 2)
df_dx = df_dz * dz_dx       # chain rule

# Verify numerically
h       = 1e-5
df_dx_numerical = (sigmoid(2*(x+h)+1) - sigmoid(2*(x-h)+1)) / (2*h)

print(f"f(x)         = sigmoid(2×{x} + 1) = {f:.6f}")
print(f"df/dz        = sigmoid'(z)        = {df_dz:.6f}")
print(f"dz/dx        = d(2x+1)/dx         = {dz_dx:.6f}")
print(f"df/dx (chain)= {df_dz:.6f} × {dz_dx:.6f}   = {df_dx:.6f}")
print(f"df/dx (numer)= {df_dx_numerical:.6f}  ← matches ✓")

# ── Now with three chained functions ─────────────────────────────────
# L = MSE(relu(Wx + b), y)  — loss through linear → relu → MSE

W, b, y = 0.8, -0.2, 1.0
x_in    = 0.5

# Forward
z_val  = W * x_in + b        # linear
a_val  = max(0, z_val)        # relu
loss   = (a_val - y) ** 2     # MSE (single sample, no mean)

# Backward — chain rule at each step
dL_da  = 2 * (a_val - y)                        # ∂L/∂a
da_dz  = 1.0 if z_val > 0 else 0.0              # ∂a/∂z (relu derivative)
dz_dW  = x_in                                    # ∂z/∂W
dz_db  = 1.0                                     # ∂z/∂b

dL_dW  = dL_da * da_dz * dz_dW   # chain: ∂L/∂W
dL_db  = dL_da * da_dz * dz_db   # chain: ∂L/∂b

print(f"
Three-function chain: L = MSE(relu(Wx+b), y)")
print(f"  Forward:   z={z_val:.3f}  a={a_val:.3f}  L={loss:.3f}")
print(f"  ∂L/∂a    = 2(a−y)          = {dL_da:.4f}")
print(f"  ∂a/∂z    = relu'(z)         = {da_dz:.4f}")
print(f"  ∂z/∂W    = x               = {dz_dW:.4f}")
print(f"  ∂L/∂W    = {dL_da:.3f} × {da_dz:.3f} × {dz_dW:.3f} = {dL_dW:.4f}")

import numpy as np

np.random.seed(42)

# ── Full backprop from scratch — every gradient explicit ───────────────

def relu(z):       return np.maximum(0, z)
def relu_d(z):     return (z > 0).astype(float)
def sigmoid(z):    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))
def sigmoid_d(z):  s = sigmoid(z); return s * (1 - s)

def init(n_in, n_out):
    return np.random.randn(n_in, n_out) * np.sqrt(2.0 / n_in)

# Network: 3 inputs → 4 hidden (relu) → 2 hidden (relu) → 1 output (sigmoid)
# Task: binary classification
W1, b1 = init(3, 4), np.zeros((1, 4))
W2, b2 = init(4, 2), np.zeros((1, 2))
W3, b3 = init(2, 1), np.zeros((1, 1))

def forward(X):
    """Forward pass — cache everything needed for backprop."""
    Z1 = X  @ W1 + b1;  A1 = relu(Z1)
    Z2 = A1 @ W2 + b2;  A2 = relu(Z2)
    Z3 = A2 @ W3 + b3;  A3 = sigmoid(Z3)   # probability output
    return A3, {'X': X, 'Z1': Z1, 'A1': A1,
                'Z2': Z2, 'A2': A2, 'Z3': Z3, 'A3': A3}

def binary_cross_entropy(y_pred, y_true):
    eps = 1e-8
    return -np.mean(y_true * np.log(y_pred + eps)
                    + (1 - y_true) * np.log(1 - y_pred + eps))

def backward(y_true, cache, lr=0.01):
    """
    Backprop — compute all gradients, update all weights.
    Returns gradients for inspection.
    """
    global W1, b1, W2, b2, W3, b3
    n  = y_true.shape[0]
    A3 = cache['A3']

    # ── Loss gradient ──────────────────────────────────────────────────
    # BCE loss: ∂L/∂A3 = -(y/A3 - (1-y)/(1-A3)) / n
    # For sigmoid output with BCE this simplifies beautifully to:
    dA3 = (A3 - y_true) / n          # (n, 1)

    # ── Layer 3 backward ──────────────────────────────────────────────
    # Z3 = A2 @ W3 + b3,  A3 = sigmoid(Z3)
    # dZ3 = dA3 × sigmoid'(Z3)   ← through activation
    # But sigmoid'(Z3) = A3(1-A3), and dA3 already equals (A3-y)/n
    # so dZ3 = dA3 (the simplification from combining BCE + sigmoid)
    dZ3 = dA3                         # (n, 1)
    dW3 = cache['A2'].T @ dZ3         # (2, 1)  ← gradient for W3
    db3 = dZ3.sum(axis=0, keepdims=True)
    dA2 = dZ3 @ W3.T                  # (n, 2)  ← pass to layer 2

    # ── Layer 2 backward ──────────────────────────────────────────────
    dZ2 = dA2 * relu_d(cache['Z2'])   # (n, 2)  ← through relu
    dW2 = cache['A1'].T @ dZ2         # (4, 2)  ← gradient for W2
    db2 = dZ2.sum(axis=0, keepdims=True)
    dA1 = dZ2 @ W2.T                  # (n, 4)  ← pass to layer 1

    # ── Layer 1 backward ──────────────────────────────────────────────
    dZ1 = dA1 * relu_d(cache['Z1'])   # (n, 4)
    dW1 = cache['X'].T  @ dZ1         # (3, 4)  ← gradient for W1
    db1_g = dZ1.sum(axis=0, keepdims=True)

    # ── Gradient descent update ────────────────────────────────────────
    W3 -= lr * dW3;  b3 -= lr * db3
    W2 -= lr * dW2;  b2 -= lr * db2
    W1 -= lr * dW1;  b1 -= lr * db1_g

    return {'dW1': dW1, 'dW2': dW2, 'dW3': dW3}

# ── Train on Brex-style loan default data ─────────────────────────────
np.random.seed(42)
n = 1000
X_raw = np.column_stack([
    np.random.normal(0, 1, n),   # credit score (standardised)
    np.random.normal(0, 1, n),   # emi-to-income ratio
    np.random.normal(0, 1, n),   # employment years
])
y_raw = ((X_raw[:, 0] < -0.3) & (X_raw[:, 1] > 0.4)).astype(float).reshape(-1, 1)

print("Training 3-layer network with manual backprop:")
print(f"{'Epoch':>7} {'Loss':>10} {'Accuracy':>10}")
print("─" * 32)

for epoch in range(300):
    # Mini-batch
    idx = np.random.choice(n, 64, replace=False)
    Xb, yb = X_raw[idx], y_raw[idx]

    # Forward + backward
    pred, cache = forward(Xb)
    grads = backward(yb, cache, lr=0.05)

    if (epoch + 1) % 60 == 0:
        pred_all, _ = forward(X_raw)
        loss = binary_cross_entropy(pred_all, y_raw)
        acc  = ((pred_all > 0.5) == y_raw).mean()
        print(f"  {epoch+1:>5}  {loss:>10.6f}  {acc:>10.4f}")

# ── Gradient magnitudes — sanity check ────────────────────────────────
print(f"
Gradient magnitudes (should be non-zero, not exploding):")
for name, g in grads.items():
    print(f"  {name}: mean={np.abs(g).mean():.6f}  max={np.abs(g).max():.6f}")

import numpy as np

np.random.seed(42)

def sigmoid(z):   return 1 / (1 + np.exp(-np.clip(z, -500, 500)))
def sigmoid_d(z): s = sigmoid(z); return s * (1 - s)
def relu(z):      return np.maximum(0, z)
def relu_d(z):    return (z > 0).astype(float)
def tanh(z):      return np.tanh(z)
def tanh_d(z):    return 1 - np.tanh(z)**2

# ── Simulate gradient flow through 10 layers ──────────────────────────
n_layers  = 10
n_neurons = 32
X_test    = np.random.randn(16, n_neurons)

print("Gradient magnitude at each layer (starting from output):")
print(f"{'Layer':>8}  {'sigmoid':>12}  {'tanh':>12}  {'relu':>12}")
print("─" * 50)

for activation_name, act_fn, act_d in [
    ('sigmoid', sigmoid, sigmoid_d),
    ('tanh',    tanh,    tanh_d),
    ('relu',    relu,    relu_d),
]:
    # Simulate forward pass and track gradient magnitudes
    A = X_test.copy()
    Zs = []
    for _ in range(n_layers):
        W = np.random.randn(n_neurons, n_neurons) * 0.1
        Z = A @ W
        Zs.append(Z)
        A = act_fn(Z)

    # Simulate backward pass — gradient starts at 1.0 at output
    dA = np.ones_like(A)
    grad_norms = []
    for i in reversed(range(n_layers)):
        dZ = dA * act_d(Zs[i])
        grad_norms.insert(0, np.abs(dZ).mean())
        W  = np.random.randn(n_neurons, n_neurons) * 0.1
        dA = dZ @ W.T

    if activation_name == 'sigmoid':
        print("Sigmoid:")
        for i, g in enumerate(grad_norms):
            bar = '█' * int(min(g * 200, 30))
            print(f"  Layer {i+1:2d}: {g:.8f}  {bar}")
    elif activation_name == 'relu':
        print("
ReLU:")
        for i, g in enumerate(grad_norms):
            bar = '█' * int(min(g * 20, 30))
            print(f"  Layer {i+1:2d}: {g:.6f}  {bar}")

# ── Solutions to vanishing gradients ──────────────────────────────────
print("
Solutions to vanishing gradients:")
print("  1. Use ReLU (or LeakyReLU) instead of sigmoid/tanh in hidden layers")
print("  2. Use He initialisation: scale = sqrt(2/n_inputs)")
print("  3. Use Batch Normalisation — normalises activations between layers")
print("  4. Use residual connections (ResNets) — skip connections bypass layers")
print("  5. Use gradient clipping for exploding gradients: clip_grad_norm_(params, 1.0)")

import numpy as np

np.random.seed(42)

def relu(z):    return np.maximum(0, z)
def relu_d(z):  return (z > 0).astype(float)
def sigmoid(z): return 1/(1+np.exp(-np.clip(z,-500,500)))

# Small network for gradient checking
W1 = np.random.randn(3, 4) * 0.1;  b1 = np.zeros((1, 4))
W2 = np.random.randn(4, 1) * 0.1;  b2 = np.zeros((1, 1))

X_gc = np.random.randn(5, 3)
y_gc = np.random.randint(0, 2, (5, 1)).astype(float)

def forward_gc(W1, b1, W2, b2, X):
    Z1 = X  @ W1 + b1;  A1 = relu(Z1)
    Z2 = A1 @ W2 + b2;  A2 = sigmoid(Z2)
    return A2, {'X':X,'Z1':Z1,'A1':A1,'Z2':Z2,'A2':A2}

def loss_gc(W1, b1, W2, b2, X, y):
    pred, _ = forward_gc(W1, b1, W2, b2, X)
    eps = 1e-8
    return -np.mean(y*np.log(pred+eps) + (1-y)*np.log(1-pred+eps))

def backward_gc(cache, y):
    n   = y.shape[0]
    dA2 = (cache['A2'] - y) / n
    dW2 = cache['A1'].T @ dA2
    db2 = dA2.sum(axis=0, keepdims=True)
    dA1 = dA2 @ W2.T
    dZ1 = dA1 * relu_d(cache['Z1'])
    dW1 = cache['X'].T @ dZ1
    db1 = dZ1.sum(axis=0, keepdims=True)
    return dW1, db1, dW2, db2

# ── Analytical gradients ───────────────────────────────────────────────
_, cache = forward_gc(W1, b1, W2, b2, X_gc)
dW1_a, db1_a, dW2_a, db2_a = backward_gc(cache, y_gc)

# ── Numerical gradients ────────────────────────────────────────────────
h = 1e-5

def numerical_grad(param, param_name):
    grad = np.zeros_like(param)
    it = np.nditer(param, flags=['multi_index'])
    while not it.finished:
        idx = it.multi_index
        orig = param[idx]

        param[idx] = orig + h
        loss_plus  = loss_gc(W1, b1, W2, b2, X_gc, y_gc)
        param[idx] = orig - h
        loss_minus = loss_gc(W1, b1, W2, b2, X_gc, y_gc)
        param[idx] = orig

        grad[idx] = (loss_plus - loss_minus) / (2 * h)
        it.iternext()
    return grad

dW1_n = numerical_grad(W1, 'W1')
dW2_n = numerical_grad(W2, 'W2')
db1_n = numerical_grad(b1, 'b1')
db2_n = numerical_grad(b2, 'b2')

# ── Compare ───────────────────────────────────────────────────────────
def relative_error(a, n):
    return np.abs(a - n).max() / (np.abs(a).max() + np.abs(n).max() + 1e-10)

print("Gradient check results (relative error < 1e-5 = correct):")
for name, analytical, numerical in [
    ('dW1', dW1_a, dW1_n), ('db1', db1_a, db1_n),
    ('dW2', dW2_a, dW2_n), ('db2', db2_a, db2_n),
]:
    err = relative_error(analytical, numerical)
    status = '✓ correct' if err < 1e-4 else '✗ WRONG — bug in backprop!'
    print(f"  {name}: relative error = {err:.2e}  {status}")

import torch
import torch.nn as nn
import numpy as np

torch.manual_seed(42)

# ── Manual backprop vs PyTorch autograd — same network, same result ───

# Data
X = torch.FloatTensor([[0.5, -0.3, 0.8]])   # (1, 3)
y = torch.FloatTensor([[1.0]])               # (1, 1)

# ── Manual backprop (numpy) ────────────────────────────────────────────
np.random.seed(42)
W1_np = np.random.randn(3, 4) * 0.1
b1_np = np.zeros((1, 4))
W2_np = np.random.randn(4, 1) * 0.1
b2_np = np.zeros((1, 1))

X_np, y_np = X.numpy(), y.numpy()
Z1 = X_np @ W1_np + b1_np
A1 = np.maximum(0, Z1)
Z2 = A1  @ W2_np + b2_np
A2 = 1 / (1 + np.exp(-Z2))

dA2  = (A2 - y_np)
dW2_manual = A1.T  @ dA2
dA1  = dA2 @ W2_np.T
dZ1  = dA1 * (Z1 > 0).astype(float)
dW1_manual = X_np.T @ dZ1

# ── PyTorch autograd ───────────────────────────────────────────────────
torch.manual_seed(42)
W1_pt = torch.FloatTensor(W1_np.copy()).requires_grad_(True)
b1_pt = torch.zeros(1, 4, requires_grad=True)
W2_pt = torch.FloatTensor(W2_np.copy()).requires_grad_(True)
b2_pt = torch.zeros(1, 1, requires_grad=True)

# Forward pass — PyTorch records every operation
Z1_pt = X @ W1_pt + b1_pt
A1_pt = torch.relu(Z1_pt)
Z2_pt = A1_pt @ W2_pt + b2_pt
A2_pt = torch.sigmoid(Z2_pt)

loss_pt = nn.BCELoss()(A2_pt, y)

# Backward pass — ONE LINE replaces our entire backward() function
loss_pt.backward()

# ── Compare gradients ──────────────────────────────────────────────────
print("Manual vs PyTorch autograd gradients:")
print(f"  dW1 max difference: {np.abs(dW1_manual - W1_pt.grad.numpy()).max():.2e}")
print(f"  dW2 max difference: {np.abs(dW2_manual - W2_pt.grad.numpy()).max():.2e}")
print("  (should be near 0 — same math, different implementation)")

# ── torch.no_grad() — disable gradient tracking for inference ─────────
print("
torch.no_grad() disables the computational graph:")
with torch.no_grad():
    pred = torch.sigmoid(X @ W1_pt + b1_pt @ W2_pt + b2_pt)
    # No graph is built — faster, uses less memory
    # W1_pt.grad will NOT be updated after this block
    print(f"  Prediction (no grad): {pred.item():.4f}")

# ── optimizer.zero_grad() — why gradients must be cleared ─────────────
print("
Why optimizer.zero_grad() is required every step:")
print("  PyTorch ACCUMULATES gradients by default — .grad += new_grad")
print("  Without zero_grad(), gradients from step t are added to step t+1")
print("  This is intentional for gradient accumulation techniques")
print("  But for standard training: always call zero_grad() before backward()")

1. Print the gradient norm of every layer's weights right after one backward() call — iterate model.named_parameters() and print each parameter's grad norm.

2. Compare the first layer's norm to the last layer's. A ratio smaller than roughly 1e-3 points to vanishing gradients. A ratio larger than roughly 1e3, or any norm in the thousands, points to exploding gradients.

3. If gradients look healthy at every layer but the loss still will not move, the problem is almost never backprop itself — check the learning rate, the loss function choice, and whether labels are shaped and typed correctly.

import torch

# ── A custom autograd.Function — the one place backprop is hand-written ─
# Real use case: quantisation-aware training. round() has zero gradient
# almost everywhere, so naive autograd would silently zero out every
# upstream gradient. The "straight-through estimator" fixes this by
# treating round() as the identity function purely for backward purposes.

class StraightThroughRound(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x):
        return torch.round(x)

    @staticmethod
    def backward(ctx, grad_output):
        # Pass the incoming gradient straight through, unchanged —
        # pretend the forward op were the identity function.
        return grad_output

ste_round = StraightThroughRound.apply

x = torch.tensor([1.2, 2.7, 3.5], requires_grad=True)
y = ste_round(x)
loss = (y - torch.tensor([1.0, 3.0, 4.0])).pow(2).sum()
loss.backward()

print(f"Forward pass (rounded values): {y.tolist()}")
print(f"Gradient w.r.t. input:         {x.grad.tolist()}")
print("
Without the custom backward(), x.grad would be all zeros —")
print("torch.round() has zero derivative almost everywhere, and plain")
print("autograd would faithfully, and unhelpfully, report exactly that.")

They solve two different problems that happen to run back to back. Backprop answers "how much does each weight contribute to the error?" — it computes gradients, nothing more. Gradient descent answers "given these gradients, how should I update the weights?" You could compute gradients with backprop and then update weights with a completely different optimiser (Adam, RMSprop, SGD with momentum) — backprop does not change. It is the gradient-computation half of training, not the whole thing.

Mathematically backprop works at any depth — the chain rule does not care how many layers there are. Practically, depth is exactly what breaks it: each layer's derivative gets multiplied into the gradient, and with sigmoid/tanh those derivatives are consistently below 1, so a 20-layer network can vanish the gradient to numerical zero before it reaches the early layers. Every major deep learning breakthrough since 2015 — ResNets' skip connections, batch normalisation, careful initialisation — exists specifically to keep very deep networks trainable despite this, not because backprop "handles" depth for free.

Backprop only supplies a direction (the gradient); gradient descent then takes a step in that direction. Neural network loss surfaces are highly non-convex, so this process is only guaranteed to reach *a* local minimum, never provably the global one. In practice, modern research suggests that in large, overparameterised networks most local minima found this way perform comparably — but that is an empirical observation about the loss landscape, not a guarantee backprop itself provides.

Autograd computes every gradient for you — you will rarely, if ever, write dW = A_prev.T @ dZ by hand in production code. But "rarely write it" is not "never need to understand it": every one of the errors in this module (dead ReLUs, NaN losses, gradient explosions) is invisible to autograd itself — it faithfully computes whatever gradient the math produces, including a completely broken one. Debugging those requires knowing what backprop is actually doing under the hood, which is exactly why this module exists even though nobody hand-codes it day to day.

It is the identical multivariable chain rule from a calculus course — nothing about it is neural-network-specific. What backprop contributes is an efficient *algorithm* for applying that chain rule to a deeply nested composition of functions: instead of re-deriving the derivative of the whole network with respect to every weight from scratch (which would recompute the same intermediate terms over and over), it caches intermediate gradients during a single backward pass — the same idea as dynamic programming applied to calculus.

A good answer avoids equations entirely: "The network makes a prediction, we measure how wrong it was, and backpropagation is the process of tracing that error backward through the network to figure out exactly how much each individual connection contributed to the mistake — like assigning blame after a failed project by tracing the decision chain back to its source. Once we know each connection's share of the blame, we nudge it slightly to make that mistake less likely next time."

They happen because the chain rule multiplies each layer's local derivative into the gradient flowing backward, and sigmoid/tanh derivatives are bounded well below 1 (sigmoid tops out at 0.25) — across many layers that product shrinks toward zero, so early layers effectively stop learning. Fixes: switch to ReLU-family activations (derivative is exactly 1 for positive inputs, so nothing shrinks), add residual/skip connections so gradients have a direct path around the multiplication, and use batch or layer normalisation to keep activations in a well-behaved range at every layer.

Backprop is a gradient-computation algorithm: given the network and the loss, it produces ∂L/∂W for every weight, using the chain rule and a single backward pass. Gradient descent is an optimisation algorithm: given those gradients, it decides how to update the weights (W ← W − learning_rate × gradient), and there are many variants (SGD, Adam, RMSprop) that all consume the same backprop-computed gradients differently. Conflating them is the single most common mix-up on this topic.

PyTorch accumulates gradients into .grad by default instead of overwriting them — useful on purpose for cases like gradient accumulation across several mini-batches (simulating a larger batch size than fits in memory) or RNNs where you sometimes want gradients to add up across time steps. But for standard training, forgetting zero_grad() means each step's gradient gets added on top of the previous step's leftover gradient, silently corrupting every update after the first — one of the most common real bugs in PyTorch code, and worth mentioning unprompted in an interview.

Gradient checking: compute the gradient two independent ways and compare them. The analytical gradient comes from your backprop implementation. The numerical gradient comes from the definition of a derivative directly — nudge one parameter by a tiny epsilon (e.g. 1e-7), rerun the forward pass, and measure how much the loss changed: (loss(w+ε) − loss(w−ε)) / (2ε). Compute the relative error between the two gradients; below roughly 1e-5–1e-7 means the implementation is correct. This should be done once per new architecture, on a small toy network, then removed — it is far too slow to run during actual training.