MLP & Backpropagation — the Core Built by Hand

Python · NumPy · PyTorch  ·  Download MLP module

The Network, and How It Learns

A multilayer perceptron stacks linear maps and nonlinear activations: each layer computes z(l)=a(l−1)W(l)+b(l)z^{(l)} = a^{(l-1)}W^{(l)} + b^{(l)} then a(l)=ϕ(z(l))a^{(l)} = \phi(z^{(l)}). The nonlinearity is not decoration — without it two linear layers compose to x(W(1)W(2))x(W^{(1)}W^{(2)}), still a single linear map no more expressive than logistic regression. ReLU on the hidden layers keeps gradients healthy where sigmoid and tanh saturate and stall them; a sigmoid stays on the output only because we want a probability there.

z(l)=a(l−1)W(l)+b(l),a(l)=ϕ(z(l)),ϕ=ReLU  (hidden),  σ  (output)z^{(l)} = a^{(l-1)}W^{(l)} + b^{(l)}, \qquad a^{(l)} = \phi\big(z^{(l)}\big), \qquad \phi = \mathrm{ReLU} \;\text{(hidden)},\; \sigma \;\text{(output)}

Training needs ∂L/∂W(l)\partial L/\partial W^{(l)} for every layer, and backpropagation is nothing more exotic than the chain rule applied layer by layer. The sigmoid–cross-entropy pairing makes the output error collapse to the residual, and from there the error signal is pushed backwards by undoing the linear map and the activation in turn. That three-line recursion is the entire algorithm behind every network in this section.

δ(L)=a(L)−yn,δ(l)=(δ(l+1)W(l+1)⊤)⊙ReLU′(z(l)),∂L∂W(l)=a(l−1)⊤δ(l)\delta^{(L)} = \frac{a^{(L)} - y}{n}, \qquad \delta^{(l)} = \big(\delta^{(l+1)}W^{(l+1)\top}\big) \odot \mathrm{ReLU}'\big(z^{(l)}\big), \qquad \frac{\partial L}{\partial W^{(l)}} = a^{(l-1)\top}\delta^{(l)}

Proving It Correct

Correctness is proved, not assumed. Nudge one weight by ±ε\pm\varepsilon, measure the actual change in loss, compare it to the analytic gradient: agreement at 10⁻⁸ means the hand-written backward pass computes the true gradient. A single wrong transpose or a missing activation derivative shows up here as ~10⁻¹.

The check contains a trap worth keeping, because it is a classic way to fool yourself. The differenced objective must be exactly the function backward differentiates. Since backward adds λW\lambda W to each weight gradient — the derivative of λ2∥W∥2\tfrac{\lambda}{2}\lVert W\rVert^2 — that penalty has to appear in the differenced loss too. Omit it and the check reports a relative error of 0.7 on gradients that are perfectly correct, and, worse, it silently validates only the unregularised network rather than the one actually trained. The check here runs at λ=0\lambda = 0, 10−410^{-4} and 10−210^{-2}, passing at 10⁻⁸ or better throughout.

weight decay λmax relative error (analytic vs finite difference)with the penalty omitted from the objective
01.9×10⁻⁸1.9×10⁻⁸ — no penalty to omit
10⁻⁴ (used for training)7.7×10⁻⁸6.8×10⁻¹ — spurious failure
10⁻²1.5×10⁻⁹1.0 — spurious failure

What the Framework Actually Adds

Writing backward by hand does not scale, and PyTorch removes the burden: you write only the forward computation, autograd records it and replays it in reverse. The interesting question is what that actually changes. The PyTorch net comes out 0.016 AUC ahead of the from-scratch one — so is the hand-written maths wrong?

No, and the decomposition shows it. Rebuild the PyTorch net with our optimiser and our initialisation and it lands within 0.0011 AUC of the hand-written network. That is the real test of the mathematics, and it passes. The remaining gap is settings, and splitting it is instructive: +0.0098 from initialisation (PyTorch's default beating He here) and +0.0071 from Adam over plain SGD. Initialisation contributes at least as much as the optimiser — which is why it belongs among the mechanics rather than in a footnote. What a framework buys is not better gradients; it is autograd for any architecture, plus defaults you would otherwise have to discover yourself.

credit default, same architecture, 60 epochsOOS AUC
from-scratch — SGD lr 0.3, He init0.7495
PyTorch — SGD lr 0.3, He init — matched setup0.7484 — 0.0011 apart: the maths is identical
PyTorch — SGD lr 0.3, PyTorch default init0.7582 — +0.0098 from initialisation
PyTorch — Adam lr 10⁻³, PyTorch default init0.7653 — +0.0071 more from the optimiser

The Honest Tabular Verdict

Does the network win on tabular data? A clean and important no — and the comparison is worth making per-model rather than lumping the tree methods together. On credit default the MLP reaches AUC 0.765, comfortably above the logistic baseline's 0.715, but below both the random forest (0.775) and XGBoost (0.774). On California housing a single seed gives RMSE 0.518, which flatters it: across five seeds the range is 0.518–0.532 (mean 0.525, sd 0.005), which brackets the random forest's 0.523. So on regression the net and the forest are a tie once seed variation is admitted, and quoting one run either way is reading noise. Against XGBoost's 0.494 there is no ambiguity.

methodcredit AUC (higher better)California RMSE (lower better)
logistic / linear0.7150.737
neural net (MLP)0.7650.525 — mean of 5 seeds, range 0.518–0.532
random forest0.7750.523 — tied with the net
XGBoost0.7740.494

This is not a tuning failure but a well-replicated regularity (Grinsztajn, Oyallon & Varoquaux, 2022). Tabular columns have no spatial or temporal ordering for a network's inductive biases to exploit — permute the features and nothing changes, whereas a convolutional net's entire design assumes neighbouring pixels are related. Trees handle mixed scales, skew and irrelevant features natively, since each split is scale-free, while networks need careful standardisation and still struggle with uninformative inputs. The lesson is not that neural nets are bad; it is to match the model to the structure of the data — which is exactly what the rest of this section supplies.

Where this sits

The L2 weight decay used here is the Ridge penalty from Ridge, Lasso & Elastic Net, now applied to network weights. The tree baselines come from Random Forests and XGBoost, LightGBM & CatBoost on identical splits. And the uncertainty theme that closes this section — MC-dropout and deep ensembles — connects to BART and Gaussian Processes & Splines, where the same question is answered from the Bayesian side.

Notebook

Python A full data description, the forward pass with shapes made concrete, backpropagation derived and implemented by hand with a gradient check at three levels of weight decay, the training mechanics, the PyTorch rebuild decomposed into initialisation and optimiser effects, and the tabular scoreboard with seed variation (Python) Open →

Downloads

nn_mlp.py Multilayer perceptron with manual backpropagation — forward pass with cached activations, the δ-recursion by hand, He initialisation, mini-batch SGD with L2 weight decay, and a gradient checker valid at any penalty strength (NumPy) credit_default.csv Default of credit-card clients (UCI) — the classification task, shared across the whole machine-learning arc cali_housing.csv California housing — the regression task, shared across the whole machine-learning arc

MLP Module — Source Code

"""From-scratch multilayer perceptron with manual backpropagation.

No autograd: the forward pass caches activations, and `backward` implements the
chain rule by hand to produce the exact gradients that PyTorch would compute
automatically.  ReLU hidden units, sigmoid output, binary cross-entropy loss,
He initialisation, mini-batch SGD, optional L2 weight decay.

The point of the notebook is that backpropagation *is* the chain rule applied
layer by layer:

    output layer   delta_L = (a_L - y) / n                 (BCE + sigmoid)
    hidden layer   delta_l = (delta_{l+1} W_{l+1}^T) * ReLU'(z_l)
    gradients      dW_l = a_{l-1}^T delta_l,   db_l = sum_i delta_l

`gradient_check` verifies these analytic gradients against finite differences to
~1e-6 -- the standard proof that a from-scratch backprop implementation is
correct -- and the notebook confirms the trained network matches a PyTorch net
of the same architecture.
"""
import numpy as np


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


def relu_grad(z):
    return (z > 0).astype(float)


def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))


class MLP:
    """Multilayer perceptron for binary classification.

    Parameters
    ----------
    sizes : layer widths, e.g. [n_features, 32, 16, 1] (last must be 1).
    lr    : SGD learning rate.
    l2    : L2 weight-decay strength (0 = off).
    """

    def __init__(self, sizes, lr=0.05, l2=0.0, seed=0):
        self.sizes = sizes; self.lr = lr; self.l2 = l2
        rng = np.random.default_rng(seed)
        self.W = [rng.standard_normal((sizes[i], sizes[i + 1])) * np.sqrt(2.0 / sizes[i])
                  for i in range(len(sizes) - 1)]                       # He init
        self.b = [np.zeros(sizes[i + 1]) for i in range(len(sizes) - 1)]

    def forward(self, X):
        self.a = [X]; self.z = []
        A = X
        for i in range(len(self.W)):
            Z = A @ self.W[i] + self.b[i]; self.z.append(Z)
            A = sigmoid(Z) if i == len(self.W) - 1 else relu(Z)        # sigmoid out, ReLU hidden
            self.a.append(A)
        return A

    def backward(self, y):
        n = len(y); L = len(self.W)
        delta = (self.a[-1] - y.reshape(-1, 1)) / n                    # BCE+sigmoid combined gradient
        gW = [None] * L; gb = [None] * L
        for i in reversed(range(L)):
            gW[i] = self.a[i].T @ delta + self.l2 * self.W[i]
            gb[i] = delta.sum(0)
            if i > 0:
                delta = (delta @ self.W[i].T) * relu_grad(self.z[i - 1])   # chain rule into the previous layer
        return gW, gb

    @staticmethod
    def _bce(y, p):
        p = np.clip(p, 1e-7, 1 - 1e-7)
        return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))

    def fit(self, X, y, epochs=100, batch=128, seed=0, verbose=False):
        rng = np.random.default_rng(seed); n = len(y); self.loss_ = []
        for ep in range(epochs):
            idx = rng.permutation(n)
            for s in range(0, n, batch):
                bi = idx[s:s + batch]
                self.forward(X[bi]); gW, gb = self.backward(y[bi])
                for i in range(len(self.W)):
                    self.W[i] -= self.lr * gW[i]; self.b[i] -= self.lr * gb[i]
            self.loss_.append(self._bce(y, self.forward(X).ravel()))
            if verbose and ep % 10 == 0:
                print(f"epoch {ep:4d}  BCE {self.loss_[-1]:.4f}")
        return self

    def predict_proba(self, X):
        return self.forward(X).ravel()


def gradient_check(net, X, y, n_params=8, eps=1e-5, seed=0):
    """Compare analytic backprop gradients to finite-difference estimates on a
    random handful of weights.  Returns the max relative error (~1e-8 in float64).

    The differenced objective must be the SAME function `backward` differentiates.
    `backward` adds ``l2 * W`` to every weight gradient, which is the derivative of
    ``0.5 * l2 * sum(W**2)``, so that penalty belongs in the objective too -- omit it
    and the check reports a spurious failure (relative error ~0.7 at l2=1e-4) on a
    net whose gradients are perfectly correct.  Including it keeps the check valid
    for any net, not only the unregularised one.
    """
    def objective():
        loss = net._bce(y, net.forward(X).ravel())
        if net.l2:
            loss += 0.5 * net.l2 * sum(float((W ** 2).sum()) for W in net.W)
        return loss

    net.forward(X); gW, _ = net.backward(y)
    rng = np.random.default_rng(seed); errs = []
    for _ in range(n_params):
        li = rng.integers(len(net.W))
        i = rng.integers(net.W[li].shape[0]); j = rng.integers(net.W[li].shape[1])
        orig = net.W[li][i, j]
        net.W[li][i, j] = orig + eps; lp = objective()
        net.W[li][i, j] = orig - eps; lm = objective()
        net.W[li][i, j] = orig
        num = (lp - lm) / (2 * eps); ana = gW[li][i, j]
        errs.append(abs(num - ana) / max(1e-12, abs(num) + abs(ana)))
    return max(errs)

References