Neural Networks & Backpropagation

The perceptron, introduced by Frank Rosenblatt in 1958, is the atomic unit of neural computation. It takes a vector of real-valued inputs x = (x₁, x₂, …, xₙ), forms a weighted sum z = w·x + b = Σᵢ wᵢxᵢ + b, and passes it through a threshold function to produce a binary output.

The geometric intuition is sharp: the weights w and bias b define a hyperplane in n-dimensional input space. Points on one side are classified as +1, points on the other as −1. Training the perceptron means moving that hyperplane until it correctly separates the two classes.

The perceptron learning rule is beautifully simple. For each misclassified example (x, y):

w ← w + η · y · x b ← b + η · y

where η is the learning rate. This nudges the decision boundary toward the misclassified point. Rosenblatt's convergence theorem guarantees that if the data is linearly separable, this procedure converges in a finite number of steps — one of the first provable learning guarantees in history.

But therein lies the rub. Minsky and Papert (1969) proved that a single perceptron cannot learn XOR, which is not linearly separable. This observation temporarily collapsed funding for neural network research. The fix — stacking perceptrons into layers — was obvious in hindsight but required the training algorithm we develop later in this module.

import numpy as np

class Perceptron:
    def __init__(self, n_features, lr=0.1):
        self.w = np.zeros(n_features)
        self.b = 0.0
        self.lr = lr

    def predict(self, x):
        return 1 if np.dot(self.w, x) + self.b >= 0 else -1

    def fit(self, X, y, epochs=100):
        for _ in range(epochs):
            for xi, yi in zip(X, y):
                if self.predict(xi) != yi:
                    self.w += self.lr * yi * xi
                    self.b += self.lr * yi

A multilayer perceptron (MLP) arranges neurons into layers: one input layer, one or more hidden layers, and an output layer. Each neuron in a layer receives every activation from the previous layer (hence fully connected or dense).

Why does depth help? Consider universal approximation (Cybenko, 1989; Hornik, 1991): a single hidden layer with enough neurons and a non-polynomial activation can approximate any continuous function on a compact domain to arbitrary precision. But this is an existence result — the required width may be exponential in the input dimension. Depth allows a network to represent the same function with exponentially fewer parameters by building hierarchical feature detectors. An image network learns edges, then textures, then object parts, then objects — each layer reusing lower-level abstractions.

Formally, for L layers with weight matrices W¹, W², …, Wᴸ and bias vectors b¹, …, bᴸ:

a⁰ = x (input) z^l = W^l · a^(l-1) + b^l (pre-activation) a^l = σ(z^l) (activation) ŷ = a^L (output)

where σ is an elementwise nonlinearity. Without σ, any number of stacked linear layers collapses to a single linear transformation — depth would be useless. The nonlinearity is what makes the composition genuinely more expressive.

# Minimal MLP in numpy, forward pass only
import numpy as np

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

def forward(X, weights):
    """weights: list of (W, b) tuples, one per layer."""
    a = X
    caches = []
    for W, b in weights[:-1]:
        z = a @ W.T + b     # (batch, out_features)
        a = relu(z)
        caches.append((a, z))
    W, b = weights[-1]
    logits = a @ W.T + b    # final layer, no activation here
    return logits, caches

The choice of activation function profoundly affects what a network can express and how easily it trains. Here are the canonical options and their tradeoffs:

Sigmoid σ(z) = 1/(1+e⁻ᶻ): maps ℝ → (0,1). Historically popular but problematic in deep networks. Its derivative σ'(z) = σ(z)(1−σ(z)) ≤ 0.25 everywhere, so gradients shrink by at least 4× each layer — leading to the vanishing gradient problem in deep nets.

Tanh tanh(z) = (eᶻ−e⁻ᶻ)/(eᶻ+e⁻ᶻ): maps ℝ → (−1,1). Zero-centered (unlike sigmoid), so gradients have mixed sign, leading to less zig-zagging in weight updates. Still suffers vanishing gradients.

ReLU (Rectified Linear Unit) max(0, z): piecewise linear, computationally trivial, derivative is 1 for z > 0. Gradients do not vanish for positive activations, enabling much deeper networks. Downside: dying ReLU — neurons with z ≤ 0 permanently output 0 and receive zero gradient; they "die."

Leaky ReLU max(αz, z) with small α (e.g. 0.01): fixes dying ReLU by allowing a small gradient for negative inputs.

GELU (Gaussian Error Linear Unit) z·Φ(z), where Φ is the standard normal CDF: smooth approximation of ReLU used in BERT, GPT, and most modern transformers. Empirically superior on language tasks.

Softmax (output layer for classification): converts a vector of logits to a probability distribution: softmax(z)ᵢ = eᶻⁱ / Σⱼ eᶻʲ. Always used with cross-entropy loss.

import numpy as np

def sigmoid(z):     return 1 / (1 + np.exp(-z))
def sigmoid_d(z):   s = sigmoid(z); return s * (1 - s)

def tanh_d(z):      return 1 - np.tanh(z)**2

def relu(z):        return np.maximum(0, z)
def relu_d(z):      return (z > 0).astype(float)

def softmax(z):
    # numerically stable: subtract max before exp
    e = np.exp(z - z.max(axis=-1, keepdims=True))
    return e / e.sum(axis=-1, keepdims=True)

Training a network means minimising a loss function L(ŷ, y) that measures the discrepancy between predictions ŷ and true labels y over the training set.

Mean Squared Error (MSE) L = (1/N)Σᵢ ‖yᵢ − ŷᵢ‖²: natural for regression. Its gradient ∂L/∂ŷ = −2(y−ŷ)/N is smooth everywhere.

Binary Cross-Entropy L = −(1/N)Σᵢ [yᵢ log(ŷᵢ) + (1−yᵢ)log(1−ŷᵢ)]: the right loss for binary classification with a sigmoid output. It is the negative log-likelihood under a Bernoulli model, so minimising it is maximum likelihood estimation.

Categorical Cross-Entropy L = −(1/N)Σᵢ Σₖ yᵢₖ log(ŷᵢₖ): the multi-class generalisation. Combined with softmax output, its gradient simplifies beautifully: ∂L/∂z = ŷ − y (the error vector), which is the numerical gift that makes classification training stable.

The key insight: the loss is a scalar function of all the network's parameters. Gradient descent requires computing ∂L/∂wᵢⱼ for every weight in the network — potentially millions of them. Doing this naively by finite differences (perturbing each weight one at a time) would cost O(P) forward passes for P parameters. Backpropagation does it in one backward pass at cost proportional to the forward pass. That is the algorithmic miracle at the heart of deep learning.

import numpy as np

def mse_loss(y_pred, y_true):
    diff = y_pred - y_true
    return np.mean(diff**2), 2 * diff / y_true.size  # (loss, grad)

def cross_entropy_loss(logits, y_true_one_hot):
    """logits: (batch, C); y_true_one_hot: (batch, C)"""
    probs = softmax(logits)          # defined earlier
    N = logits.shape[0]
    loss = -np.sum(y_true_one_hot * np.log(probs + 1e-12)) / N
    grad = (probs - y_true_one_hot) / N   # elegant closed form
    return loss, grad

Backpropagation is the chain rule of calculus, applied systematically to a computation graph. If z = f(g(x)), then dz/dx = (dz/dg)·(dg/dx). For vector-valued functions this generalises: if z = f(h) and h = g(x), then the Jacobian satisfies ∂z/∂x = (∂z/∂h)·(∂h/∂x).

In a neural network the loss L is a scalar, so we always work with gradients (row vectors or, equivalently, transposed column vectors). Define δ^l (the "error" at layer l) as the gradient of L with respect to the pre-activation z^l:

δ^L = ∂L/∂z^L (depends on loss + output activation) δ^l = (W^(l+1))ᵀ δ^(l+1) ⊙ σ'(z^l) (backpropagation recurrence)

Here ⊙ is elementwise multiplication and σ'(z^l) is the derivative of the activation at layer l.

The weight and bias gradients are then:

∂L/∂W^l = δ^l (a^(l-1))ᵀ ∂L/∂b^l = δ^l

This is the complete algorithm. Read these equations carefully: the δ terms are computed from the output layer backward, each one depending on the next. This is why it is called backpropagation.

The computational efficiency comes from the fact that each δ^l is computed once and reused for both the weight gradient at layer l and the δ at layer l−1. No redundant computation — just the chain rule unwound layer by layer.

# Notation: L layers indexed 1..L; a[0]=X (input)
# params = [(W1,b1), (W2,b2), ..., (WL,bL)]

def backward(grad_output, caches, params):
    """
    grad_output: dL/d(a^L) -- gradient from loss
    caches: [(a^l, z^l)] stored during forward pass
    params: [(W^l, b^l)]
    Returns: list of (dW^l, db^l) in forward order
    """
    grads = [None] * len(params)
    delta = grad_output  # start from output layer

    for l in reversed(range(len(params))):
        a_prev = caches[l-1][0] if l > 0 else X_input
        _, z_l  = caches[l]
        # chain rule: delta at this layer
        delta = delta * relu_d(z_l)         # elem-wise: sigma'
        dW = delta.T @ a_prev               # outer product (summed over batch)
        db = delta.sum(axis=0)
        grads[l] = (dW, db)
        # propagate delta to previous layer
        W_l, _ = params[l]
        delta = delta @ W_l
    return grads

Gradient descent updates each parameter in the direction of steepest descent of the loss:

θ ← θ − η · ∇_θ L

where η (the learning rate) controls the step size. The three main variants differ in how much data is used to estimate the gradient:

Batch Gradient Descent: uses the entire training set to compute the exact gradient. Stable but extremely slow for large datasets.

Stochastic Gradient Descent (SGD): updates parameters after each example. Noisy but fast; the noise can actually help escape sharp local minima.

Mini-batch SGD: the practical standard. Computes gradients on batches of 32–512 examples. Leverages GPU parallelism (matrix multiplications over a batch are efficient) while providing lower variance than single-sample SGD.

Momentum: accumulates a velocity vector v = β·v − η·∇L and sets θ ← θ + v. This dampens oscillations across dimensions (helpful in ravines) and accelerates along the gradient direction.

Adam (Adaptive Moment Estimation): the default optimizer in most deep learning. Maintains per-parameter estimates of the first moment (mean) m and second moment (uncentered variance) v of the gradients:

m ← β₁m + (1−β₁)g v ← β₂v + (1−β₂)g² θ ← θ − η · (m̂ / (√v̂ + ε))

where m̂ = m/(1−β₁ᵗ) and v̂ = v/(1−β₂ᵗ) are bias-corrected estimates. Standard defaults: β₁=0.9, β₂=0.999, ε=1e-8. Adam is largely insensitive to learning rate scale, making it forgiving to tune.

class Adam:
    def __init__(self, lr=1e-3, b1=0.9, b2=0.999, eps=1e-8):
        self.lr, self.b1, self.b2, self.eps = lr, b1, b2, eps
        self.m, self.v, self.t = {}, {}, 0

    def step(self, params, grads):
        self.t += 1
        for i, (theta, g) in enumerate(zip(params, grads)):
            self.m[i] = self.b1 * self.m.get(i, 0) + (1 - self.b1) * g
            self.v[i] = self.b2 * self.v.get(i, 0) + (1 - self.b2) * g**2
            m_hat = self.m[i] / (1 - self.b1**self.t)
            v_hat = self.v[i] / (1 - self.b2**self.t)
            params[i] -= self.lr * m_hat / (v_hat**0.5 + self.eps)

The choice of initial weights is not cosmetic — it determines whether training proceeds at all. Two failure modes:

1. Weights too large: activations saturate (especially with sigmoid/tanh), gradients vanish.
2. Weights too small (e.g. all zero): every neuron in a layer computes the same function and receives the same gradient — the symmetry breaking problem. The network can never differentiate neurons in the same layer.

Xavier / Glorot Initialisation (for tanh): initialises weights ~ Uniform(−√(6/(nᵢₙ+nₒᵤₜ)), √(6/(nᵢₙ+nₒᵤₜ))). The scaling ensures that the variance of activations and gradients is preserved across layers — derived by requiring Var(aˡ) = Var(aˡ⁻¹) under a linear approximation.

He Initialisation (for ReLU): weights ~ N(0, √(2/nᵢₙ)). The factor of 2 compensates for the fact that ReLU zeros out half of its inputs, halving the variance of the output.

Both initialisations are derived from the same principle: the forward signal (and its backward gradient) should neither explode nor vanish as it traverses the network. This is a special case of the general concern of dynamical isometry — keeping the singular values of the weight matrices close to 1.

import numpy as np

def xavier_init(fan_in, fan_out):
    limit = np.sqrt(6 / (fan_in + fan_out))
    return np.random.uniform(-limit, limit, (fan_out, fan_in))

def he_init(fan_in, fan_out):
    std = np.sqrt(2 / fan_in)
    return np.random.randn(fan_out, fan_in) * std

The backpropagation recurrence δ^l = (W^(l+1))ᵀ δ^(l+1) ⊙ σ'(z^l) means that gradients are multiplied by the weight matrices and activation derivatives at every layer. In a network of depth L, the gradient at the first layer involves a product of L such factors.

Vanishing gradients: if |W·σ'| < 1 on average, this product shrinks exponentially. With sigmoid activations, σ'(z) ≤ 0.25, so gradients decay at least as fast as 0.25ᴸ — effectively zero at layer 1 for L ≥ 10. Early layers learn almost nothing; the network stalls.

Exploding gradients: if |W·σ'| > 1, the product grows exponentially, causing weights to diverge. Common in recurrent networks processing long sequences.

Mitigations:

• Use ReLU or GELU (gradient = 1 in the active region).
• Careful weight initialisation (Xavier, He).
• Batch Normalisation: normalise pre-activations within each mini-batch to zero mean and unit variance (with learnable scale/shift γ, β). This prevents saturation and keeps intermediate distributions stable, dramatically accelerating training.
• Residual / skip connections (ResNets): add a shortcut z^(l+1) = F(z^l, W) + z^l. The '+z^l' term ensures the gradient has an additive path ∂L/∂z^l = ∂L/∂z^(l+1)·(∂F/∂z^l + I), where I is the identity. This direct path prevents the gradient from vanishing, enabling networks of hundreds of layers.

# Batch Normalisation (inference mode, simplified)
import numpy as np

def batch_norm(z, gamma, beta, eps=1e-5):
    """
    z: (batch, features)
    gamma, beta: learnable scale and shift, shape (features,)
    """
    mu    = z.mean(axis=0)              # per-feature mean
    sigma = z.std(axis=0)              # per-feature std
    z_hat = (z - mu) / (sigma + eps)   # normalise
    return gamma * z_hat + beta        # rescale + shift

With all the pieces in place, the full training procedure is a loop over the following steps:

1. Sample a mini-batch (X_batch, y_batch) from the training set.
2. Forward pass: compute activations a¹, a², …, aᴸ and cache z^l, a^l for each layer.
3. Compute loss L(aᴸ, y_batch) and its gradient ∂L/∂aᴸ.
4. Backward pass: propagate gradients δ^L, δ^(L−1), …, δ¹ using the chain rule recurrence. Compute ∂L/∂W^l and ∂L/∂b^l from the deltas.
5. Update parameters using the chosen optimizer (SGD, Adam, etc.).
6. Repeat until convergence (monitored via validation loss).

Key practical details:

• Shuffling the training set each epoch prevents the network from overfitting to ordering artifacts.
• Validation split: hold out 10–20% of data to monitor overfitting. Training loss always decreases; validation loss eventually increases if the model overfits.
• Regularisation: L2 weight decay (add λ‖W‖² to the loss, equivalent to shrinking weights by λη each step) and Dropout (randomly zero activations during training with probability p, scale by 1/(1-p) at test time) both reduce overfitting.
• The combination of BatchNorm + He init + ReLU + Adam is the practical default that makes most architectures trainable out of the box.

# Complete minimal MLP training loop in numpy
import numpy as np

def train(X, y, layers=[64, 32, 10], lr=1e-3, epochs=20, batch=64):
    # Initialise weights with He init
    dims = [X.shape[1]] + layers
    params = [(he_init(dims[i], dims[i+1]),
               np.zeros(dims[i+1])) for i in range(len(dims)-1)]
    optimizer = Adam(lr=lr)

    N = X.shape[0]
    for epoch in range(epochs):
        idx = np.random.permutation(N)
        total_loss = 0.0
        for start in range(0, N, batch):
            Xb = X[idx[start:start+batch]]
            yb = y[idx[start:start+batch]]

            # Forward
            logits, caches = forward(Xb, params)

            # Loss + grad
            loss, dlogits = cross_entropy_loss(logits, yb)
            total_loss += loss

            # Backward
            param_grads = backward(dlogits, caches, params)

            # Update
            flat_p = [p for pair in params for p in pair]
            flat_g = [g for pair in param_grads for g in pair]
            optimizer.step(flat_p, flat_g)

        print(f'Epoch {epoch+1}: loss = {total_loss/(N//batch):.4f}')
    return params