From 2efeca99b83424965e51f614f6fdc9333190f215 Mon Sep 17 00:00:00 2001 From: Thomas Hart Date: Wed, 22 Jul 2026 16:29:03 +0000 Subject: [PATCH] feat: Add MLP with hand-derived backprop on toy datasets --- README.md | 19 +++ src/mlp.py | 300 ++++++++++++++++++++++++++++++++++++++++++++++ tests/test_mlp.py | 171 ++++++++++++++++++++++++++ 3 files changed, 490 insertions(+) create mode 100644 src/mlp.py create mode 100644 tests/test_mlp.py diff --git a/README.md b/README.md index 1675380..011b29a 100644 --- a/README.md +++ b/README.md @@ -25,12 +25,18 @@ Every deep learning framework is, at its core, a graph that records operations a - The `sigmoid(z) - y` gradient identity, the same residual form that least squares has - Numerically stable sigmoid (sign-branched) and log loss via `softplus(z) - y·z` (`logaddexp`) - The linear decision boundary as the level set where the logit is zero +- A multilayer perceptron: stacked affine layers with `tanh`/`relu` nonlinearities +- Backpropagation derived by hand as a delta recursion, vectorized over a minibatch (no autodiff) +- The `softmax - onehot` output-delta shortcut and `delta_prev = (delta @ Wᵀ) ⊙ act'(z)` for hidden layers +- Xavier vs He weight initialization chosen from the activation to keep signal variance stable across depth +- Learning a non-linearly-separable target (XOR, interleaved spirals) that a single linear model cannot fit ## What's implemented - **Gradient descent + autodiff-lite**: a scalar reverse-mode autograd engine (`Value`) with `+`, `*`, `**`, division, and `relu`/`tanh`/`exp`/`log`/`sigmoid`, plus an `SGD` optimizer (with momentum) and a `minimize` training loop. Gradients are checked against finite differences and against numpy for softmax cross-entropy. - **Linear regression, two ways**: `fit_normal_equation` solves ordinary least squares exactly by setting the gradient to zero and solving `(AᵀA + λR)θ = Aᵀy`, with optional ridge and a least-norm fallback when `AᵀA` is singular. `fit_sgd` fits the same model with minibatch gradient descent on standardized features, then maps the weights back to raw feature space. Both are checked to recover the true coefficients and to agree with each other, so the exact-vs-iterative tradeoff is measurable. - **Logistic regression + cross-entropy, decision boundary demo**: `src/logreg.py` trains a binary classifier by minibatch SGD, using the fact that the gradient of the mean cross-entropy with respect to the logits is exactly `sigmoid(z) - y`, the same residual form linear regression has. The sigmoid branches on the sign of the logit so `exp` never overflows, and the loss is computed as `softplus(z) - y·z` through `numpy.logaddexp` so a confidently-wrong prediction gives a large finite loss instead of `inf`. `decision_boundary` returns the straight line where the model sits at 50 percent for a two-feature problem, the level set `w·x + b = 0`. +- **Multilayer perceptron with hand-derived backprop**: `src/mlp.py` stacks affine layers with `tanh` or `relu` and a softmax head, and trains a multiclass classifier by minibatch SGD. The backward pass is written out by hand as one recursion on the per-layer delta rather than delegated to an autodiff engine: the output delta is the `softmax - onehot` residual, each hidden delta is `(delta_next @ W_nextᵀ) ⊙ act'(z)`, and the parameter gradients are `dW = a_prevᵀ @ delta` and `db = Σ delta`. Weights use He init for `relu` and Xavier for `tanh` so the signal variance holds across depth. The gradients are verified against central finite differences to a tight tolerance, and the model learns XOR and a three-arm spiral, targets a single hyperplane provably cannot separate. Ships with `make_xor` and `make_spiral` toy generators. ## Usage @@ -79,6 +85,19 @@ xs = np.linspace(-3, 3, 50) boundary = decision_boundary(model, xs) ``` +Train a small network on a target no line can separate: + +```python +from src.mlp import fit, accuracy, make_xor + +X, y = make_xor(n=400, seed=0) +model, history = fit(X, y, hidden=(16,), activation="tanh", epochs=300) + +print(accuracy(model, X, y)) # ~1.0; a linear model caps out at 0.75 here +print(history[0], history[-1]) # cross-entropy falls over training +print(model.predict_proba(X[:3])) # per-class probabilities that sum to 1 +``` + Differentiate an arbitrary scalar expression: ```python diff --git a/src/mlp.py b/src/mlp.py new file mode 100644 index 0000000..d7be74f --- /dev/null +++ b/src/mlp.py @@ -0,0 +1,300 @@ +"""Multilayer perceptron with backpropagation derived by hand. + +A single linear model draws one hyperplane, so it cannot separate classes that +are not linearly separable (XOR is the classic counterexample). Stacking affine +layers with a nonlinearity between them lets the network bend the decision +surface, and the weights are learned by backpropagation: one forward pass caches +every pre-activation and activation, then the chain rule is walked backward layer +by layer. + +The whole method is one recursion on the "delta" (the gradient of the loss with +respect to a layer's pre-activation z). With softmax + cross-entropy at the head, +the output delta collapses to `p - y_onehot`. From there each earlier delta is +`(delta_next @ W_next.T) * act'(z)`, and the parameter gradients fall out as +`dW = a_prev.T @ delta / m` and `db = mean(delta)`. Everything is vectorized over +a minibatch with numpy, but no autodiff: the gradients are written out by hand so +the matrix calculus stays visible. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from numpy.typing import NDArray + +Array = NDArray[np.float64] + + +@dataclass(frozen=True) +class _Activation: + forward: object # Callable[[Array], Array] + backward: object # Callable[[Array, Array], Array]: (z, grad_out) -> grad_in + + +def _tanh_forward(z: Array) -> Array: + return np.tanh(z) + + +def _tanh_backward(z: Array, grad_out: Array) -> Array: + t = np.tanh(z) + return grad_out * (1.0 - t * t) + + +def _relu_forward(z: Array) -> Array: + return np.maximum(z, 0.0) + + +def _relu_backward(z: Array, grad_out: Array) -> Array: + return grad_out * (z > 0.0) + + +ACTIVATIONS: dict[str, _Activation] = { + "tanh": _Activation(_tanh_forward, _tanh_backward), + "relu": _Activation(_relu_forward, _relu_backward), +} + + +def _softmax(z: Array) -> Array: + # subtract the row max first so exp never overflows on large logits + z = z - z.max(axis=1, keepdims=True) + ez = np.exp(z) + return ez / ez.sum(axis=1, keepdims=True) + + +@dataclass(frozen=True) +class MLP: + """A trained network: weights, biases, activation name, and label order. + + `weights[i]` maps layer i to layer i+1, so the last entry produces logits + over the classes. `classes` records the original label values in the order + the softmax columns correspond to, so `predict` can map argmax back to them. + """ + + weights: list[Array] + biases: list[Array] + activation: str + classes: Array + mu: Array + sigma: Array + + def logits(self, X: Array) -> Array: + act = ACTIVATIONS[self.activation] + a = (_as_2d(X) - self.mu) / self.sigma + for i, (w, b) in enumerate(zip(self.weights, self.biases, strict=True)): + z = a @ w + b + a = z if i == len(self.weights) - 1 else act.forward(z) + return a + + def predict_proba(self, X: Array) -> Array: + return _softmax(self.logits(X)) + + def predict(self, X: Array) -> Array: + idx = np.argmax(self.logits(X), axis=1) + return self.classes[idx] + + +def _as_2d(X: Array) -> Array: + arr = np.asarray(X, dtype=np.float64) + return arr.reshape(-1, 1) if arr.ndim == 1 else arr + + +def _check_xy(X: Array, y: Array) -> tuple[Array, Array]: + Xm = _as_2d(X) + yv = np.asarray(y).reshape(-1) + if Xm.shape[0] == 0: + raise ValueError("need at least one sample") + if Xm.shape[0] != yv.shape[0]: + raise ValueError(f"X has {Xm.shape[0]} rows but y has {yv.shape[0]} entries") + return Xm, yv + + +def _init_params( + sizes: list[int], activation: str, rng: np.random.Generator +) -> tuple[list[Array], list[Array]]: + """He init for relu, Xavier for tanh, so signal variance holds across depth. + + Both scale the initial weights by the fan-in; picking the wrong one for the + nonlinearity is a real way to make a deep net fail to train, so it is chosen + from the activation rather than left as a constant. + """ + weights: list[Array] = [] + biases: list[Array] = [] + for fan_in, fan_out in zip(sizes[:-1], sizes[1:], strict=True): + scale = np.sqrt((2.0 if activation == "relu" else 1.0) / fan_in) + weights.append(rng.standard_normal((fan_in, fan_out)) * scale) + biases.append(np.zeros(fan_out)) + return weights, biases + + +def cross_entropy(proba: Array, y_onehot: Array) -> float: + # clip keeps log finite when a class probability rounds to exactly 0 + p = np.clip(proba, 1e-12, 1.0) + return float(-np.mean(np.sum(y_onehot * np.log(p), axis=1))) + + +def accuracy(model: MLP, X: Array, y: Array) -> float: + Xm, yv = _check_xy(X, y) + return float(np.mean(model.predict(Xm) == yv)) + + +def fit( + X: Array, + y: Array, + hidden: tuple[int, ...] = (16, 16), + activation: str = "tanh", + lr: float = 0.1, + epochs: int = 200, + batch_size: int = 32, + seed: int = 0, +) -> tuple[MLP, list[float]]: + """Train an MLP classifier by minibatch SGD, returning it and the loss curve. + + Inputs are standardized to zero mean and unit variance so one learning rate + works across features, and the standardizer is stored on the model so it + consumes raw X at predict time. Labels may be any hashable values; they are + mapped to softmax columns and remembered on the model. Returns the trained + model and the per-epoch training cross-entropy. + """ + if activation not in ACTIVATIONS: + raise ValueError(f"unknown activation {activation!r}") + if lr <= 0.0: + raise ValueError("lr must be positive") + if epochs <= 0: + raise ValueError("epochs must be positive") + if batch_size <= 0: + raise ValueError("batch_size must be positive") + if any(h <= 0 for h in hidden): + raise ValueError("hidden layer sizes must be positive") + + Xm, yv = _check_xy(X, y) + n, d = Xm.shape + + classes, y_idx = np.unique(yv, return_inverse=True) + k = classes.shape[0] + if k < 2: + raise ValueError("need at least two classes") + y_onehot = np.eye(k)[y_idx] + + mu = Xm.mean(axis=0) + sigma = Xm.std(axis=0) + sigma = np.where(sigma == 0.0, 1.0, sigma) # constant features carry no info + Xs = (Xm - mu) / sigma + + act = ACTIVATIONS[activation] + rng = np.random.default_rng(seed) + sizes = [d, *hidden, k] + weights, biases = _init_params(sizes, activation, rng) + + history: list[float] = [] + for _ in range(epochs): + order = rng.permutation(n) + for start in range(0, n, batch_size): + idx = order[start : start + batch_size] + xb, yb = Xs[idx], y_onehot[idx] + _sgd_step(weights, biases, act, xb, yb, lr) + proba = _softmax(_forward_logits(weights, biases, act, Xs)) + history.append(cross_entropy(proba, y_onehot)) + + model = MLP( + weights=weights, + biases=biases, + activation=activation, + classes=classes, + mu=mu, + sigma=sigma, + ) + return model, history + + +def _forward_logits( + weights: list[Array], biases: list[Array], act: _Activation, x: Array +) -> Array: + a = x + last = len(weights) - 1 + for i, (w, b) in enumerate(zip(weights, biases, strict=True)): + z = a @ w + b + a = z if i == last else act.forward(z) + return a + + +def _backprop( + weights: list[Array], + biases: list[Array], + act: _Activation, + xb: Array, + yb: Array, +) -> tuple[list[Array], list[Array]]: + """Forward with caches, then the delta recursion, returning per-layer grads. + + `zs[i]` is the pre-activation of layer i and `activations[i]` its input, so + the backward pass reuses them without recomputation. The output delta is the + softmax-cross-entropy shortcut (p - y) / m, and each earlier delta multiplies + by the next weight and the local activation derivative. + """ + m = xb.shape[0] + activations: list[Array] = [xb] + zs: list[Array] = [] + a = xb + last = len(weights) - 1 + for i, (w, b) in enumerate(zip(weights, biases, strict=True)): + z = a @ w + b + zs.append(z) + a = z if i == last else act.forward(z) + activations.append(a) + + grad_w = [np.zeros_like(w) for w in weights] + grad_b = [np.zeros_like(b) for b in biases] + delta = (_softmax(zs[-1]) - yb) / m + for i in reversed(range(len(weights))): + grad_w[i] = activations[i].T @ delta + grad_b[i] = delta.sum(axis=0) + if i > 0: + delta = act.backward(zs[i - 1], delta @ weights[i].T) + return grad_w, grad_b + + +def _sgd_step( + weights: list[Array], + biases: list[Array], + act: _Activation, + xb: Array, + yb: Array, + lr: float, +) -> None: + grad_w, grad_b = _backprop(weights, biases, act, xb, yb) + for i in range(len(weights)): + weights[i] -= lr * grad_w[i] + biases[i] -= lr * grad_b[i] + + +def make_xor(n: int = 400, noise: float = 0.15, seed: int = 0) -> tuple[Array, Array]: + """XOR: four gaussian blobs at the corners, labeled by parity of the corner. + + The canonical not-linearly-separable toy. A single hyperplane gets at most + 75% here; an MLP with a hidden layer reaches ~100%. + """ + rng = np.random.default_rng(seed) + centers = np.array([[-1.0, -1.0], [1.0, 1.0], [-1.0, 1.0], [1.0, -1.0]]) + labels = np.array([0, 0, 1, 1]) + per = n // 4 + X = np.vstack([c + noise * rng.standard_normal((per, 2)) for c in centers]) + y = np.repeat(labels, per) + return X, y + + +def make_spiral( + n: int = 300, classes: int = 3, noise: float = 0.2, seed: int = 0 +) -> tuple[Array, Array]: + """Interleaved spiral arms, one class per arm. Needs a curved boundary.""" + rng = np.random.default_rng(seed) + per = n // classes + X = np.zeros((per * classes, 2)) + y = np.zeros(per * classes, dtype=np.int64) + for c in range(classes): + r = np.linspace(0.0, 1.0, per) + theta = np.linspace(c * 4.0, (c + 1) * 4.0, per) + theta = theta + noise * rng.standard_normal(per) + X[c * per : (c + 1) * per] = np.c_[r * np.sin(theta), r * np.cos(theta)] + y[c * per : (c + 1) * per] = c + return X, y diff --git a/tests/test_mlp.py b/tests/test_mlp.py new file mode 100644 index 0000000..eee1325 --- /dev/null +++ b/tests/test_mlp.py @@ -0,0 +1,171 @@ +import numpy as np +import pytest + +from src.mlp import ( + ACTIVATIONS, + MLP, + _backprop, + _forward_logits, + _softmax, + accuracy, + cross_entropy, + fit, + make_spiral, + make_xor, +) + + +def test_softmax_rows_sum_to_one_and_survive_huge_logits(): + z = np.array([[1000.0, -1000.0, 0.0], [0.0, 0.0, 0.0]]) + p = _softmax(z) + assert np.allclose(p.sum(axis=1), 1.0) + assert np.all(np.isfinite(p)) + assert p[0, 0] == pytest.approx(1.0) + assert np.allclose(p[1], 1.0 / 3.0) + + +@pytest.mark.parametrize("activation", ["tanh", "relu"]) +def test_backprop_matches_finite_differences(activation): + rng = np.random.default_rng(3) + x = rng.standard_normal((6, 3)) + y = np.eye(4)[rng.integers(0, 4, size=6)] + weights = [rng.standard_normal((3, 5)), rng.standard_normal((5, 4))] + biases = [rng.standard_normal(5), rng.standard_normal(4)] + act = ACTIVATIONS[activation] + + def loss() -> float: + logits = _forward_logits(weights, biases, act, x) + return cross_entropy(_softmax(logits), y) + + grad_w, grad_b = _backprop(weights, biases, act, x, y) + + eps = 1e-6 + for layer in range(len(weights)): + for idx in np.ndindex(weights[layer].shape): + orig = weights[layer][idx] + weights[layer][idx] = orig + eps + hi = loss() + weights[layer][idx] = orig - eps + lo = loss() + weights[layer][idx] = orig + numeric = (hi - lo) / (2.0 * eps) + assert numeric == pytest.approx(grad_w[layer][idx], abs=1e-5) + for j in range(biases[layer].shape[0]): + orig = biases[layer][j] + biases[layer][j] = orig + eps + hi = loss() + biases[layer][j] = orig - eps + lo = loss() + biases[layer][j] = orig + numeric = (hi - lo) / (2.0 * eps) + assert numeric == pytest.approx(grad_b[layer][j], abs=1e-5) + + +def test_learns_xor_that_a_linear_model_cannot(): + X, y = make_xor(n=400, seed=0) + model, history = fit(X, y, hidden=(16,), activation="tanh", epochs=300, seed=0) + assert accuracy(model, X, y) > 0.97 + assert history[-1] < history[0] + + +def test_learns_three_class_spiral(): + X, y = make_spiral(n=600, classes=3, seed=1) + model, history = fit( + X, y, hidden=(32, 32), activation="relu", lr=0.3, epochs=400, seed=1 + ) + assert accuracy(model, X, y) > 0.9 + assert len(history) == 400 + + +def test_loss_curve_trends_down(): + X, y = make_xor(n=200, seed=2) + _, history = fit(X, y, hidden=(8,), epochs=150, seed=2) + # compare averaged windows so a single noisy epoch can't fail it + assert np.mean(history[-20:]) < np.mean(history[:20]) + + +def test_predict_recovers_original_label_values(): + X, y = make_xor(n=200, seed=0) + y_labeled = np.where(y == 0, 7, 42) + model, _ = fit(X, y_labeled, hidden=(16,), epochs=200, seed=0) + preds = model.predict(X) + assert set(np.unique(preds)).issubset({7, 42}) + assert accuracy(model, X, y_labeled) > 0.95 + + +def test_predict_proba_is_a_distribution(): + X, y = make_spiral(n=300, classes=3, seed=0) + model, _ = fit(X, y, hidden=(16,), epochs=50, seed=0) + p = model.predict_proba(X) + assert p.shape == (X.shape[0], 3) + assert np.allclose(p.sum(axis=1), 1.0) + assert np.all(p >= 0.0) + + +def test_constant_feature_does_not_break_standardization(): + X, y = make_xor(n=160, seed=0) + X = np.hstack([X, np.full((X.shape[0], 1), 5.0)]) # zero-variance column + model, history = fit(X, y, hidden=(16,), epochs=200, seed=0) + assert np.all(np.isfinite(history)) + assert accuracy(model, X, y) > 0.95 + + +def test_accepts_1d_input_as_single_feature(): + rng = np.random.default_rng(0) + x = rng.standard_normal(120) + y = (x > 0).astype(np.int64) + model, _ = fit(x, y, hidden=(8,), epochs=200, seed=0) + assert accuracy(model, x, y) > 0.95 + + +def test_deterministic_under_seed(): + X, y = make_xor(n=200, seed=0) + m1, h1 = fit(X, y, hidden=(16,), epochs=50, seed=123) + m2, h2 = fit(X, y, hidden=(16,), epochs=50, seed=123) + assert h1 == h2 + pairs = zip(m1.weights, m2.weights, strict=True) + assert all(np.array_equal(a, b) for a, b in pairs) + + +def test_empty_input_raises(): + with pytest.raises(ValueError, match="at least one sample"): + fit(np.empty((0, 2)), np.array([]), epochs=1) + + +def test_mismatched_lengths_raise(): + with pytest.raises(ValueError, match="rows but"): + fit(np.zeros((4, 2)), np.zeros(3), epochs=1) + + +def test_single_class_raises(): + with pytest.raises(ValueError, match="two classes"): + fit(np.zeros((5, 2)), np.zeros(5), epochs=1) + + +@pytest.mark.parametrize( + "kwargs,match", + [ + ({"activation": "sigmoid"}, "unknown activation"), + ({"lr": 0.0}, "lr must be positive"), + ({"epochs": 0}, "epochs must be positive"), + ({"batch_size": 0}, "batch_size must be positive"), + ({"hidden": (16, 0)}, "hidden layer sizes"), + ], +) +def test_bad_hyperparameters_raise(kwargs, match): + X, y = make_xor(n=40, seed=0) + with pytest.raises(ValueError, match=match): + fit(X, y, **kwargs) + + +def test_model_is_frozen(): + model = MLP( + weights=[np.zeros((2, 2))], + biases=[np.zeros(2)], + activation="tanh", + classes=np.array([0, 1]), + mu=np.zeros(2), + sigma=np.ones(2), + ) + with pytest.raises(AttributeError): + model.activation = "relu" # type: ignore[misc]