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 CRED-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()")