Module 1: Machine learning foundations

An engineer has a function in mind and cannot write it down. The sag of a scaffold as a function of its material and geometry can be computed by a solver, but the solver took a person-year to build. The probability that a hazard is controllable, as a function of the situation, cannot be computed at all; it is estimated by people. Machine learning is the practice of obtaining such a function from examples of its inputs and outputs, when the examples are cheaper than the theory.

Formally, there is an unknown function f^\ast and a data distribution P over pairs (\mathbf{x}, y) of inputs \mathbf{x} \in \R^d and outputs y. The function f^\ast is the best possible prediction of y from \mathbf{x}; the outputs scatter around it because the inputs do not determine them completely. We see a sample \mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^{N} drawn from P, and we want a function f that predicts y from \mathbf{x} well on inputs we have not seen. The last clause is the whole subject. Predicting the examples we already have is a lookup table.

Expected risk and empirical risk

A loss \ell(f(\mathbf{x}), y) \ge 0 says how wrong one prediction is, for example the squared error (f(\mathbf{x}) - y)^2. What we want small is the expected risk, the average loss over the distribution the inputs will come from:

It cannot be computed, because P is unknown. What can be computed, for a model f_\theta with parameters \theta, is the empirical risk, the average loss over the sample:

For parameters fixed before the sample is drawn, \mathcal{L}(\theta) is an unbiased estimate of R(f_\theta): each term has expectation R(f_\theta). Fitting destroys that property. The optimiser chooses \hat\theta to make \mathcal{L} small on these N examples, noise included, so the training loss at \hat\theta is biased low. For an optimiser that finds the minimum the argument takes three steps. Write R(\theta) for R(f_\theta) and let \theta^\circ minimise it over the model family. Then \mathcal{L}(\hat\theta) \le \mathcal{L}(\theta^\circ) for every sample, because \hat\theta minimises \mathcal{L}; \E[\mathcal{L}(\theta^\circ)] = R(\theta^\circ), because \theta^\circ does not depend on the sample; and R(\theta^\circ) \le R(\hat\theta), because \theta^\circ is the best. Chained:

The gap between the training loss and the expected risk is the subject of Section 8, which measures how it grows with the flexibility of the model, and of Section 10, which estimates the expected risk without fooling yourself.

Three settings

• Supervised learning. y is given for every example: regression when y is continuous, classification when it is one of K labels. Most of this series.
• Unsupervised learning. No y. The task is to find structure: clusters, low-dimensional coordinates, densities. k-means and principal component analysis appear in Section 11; the neural forms are in Module 05.
• Reinforcement learning. No y either, but a reward that arrives after actions. Module 09 uses it to train language models.

The five ingredients

Every learning method, from a straight line to a language model, is specified by five choices. Write them down for any method you meet; a method you cannot state this way is one you do not understand yet.

1. Data. \mathcal{D}, and the distribution P it came from. The second is the one people forget. Everything learned is learned about P, so a model is only as good as the match between P and the inputs it will meet in service.
2. Model. A family of functions f_\theta indexed by parameters \theta, also called the hypothesis class. It fixes what can be learned at all: a line has two parameters and can only ever be a line; a language model has billions.
3. Loss. \ell \ge 0 and its average over the data, the empirical risk \mathcal{L}(\theta). It decides which mistakes count and how much; Section 5 shows that it also states an assumption about the noise.
4. Optimiser. A procedure that finds \theta making \mathcal{L} small, almost always a variant of gradient descent (Sections 3 and 4). It decides whether the \theta the loss asks for is actually found, and at what cost.
5. Evaluation. A measurement of how well f_\theta does on data it was not fitted to: an estimate of R, not the loss on \mathcal{D}, which the optimiser has already minimised.

Inputs are numbers

A model sees a feature vector \mathbf{x} \in \R^d. A measured quantity is already a number; a category, such as a material grade, a supplier or a defect type, is not. Encode a K-way category as one-hot: K binary columns, with a 1 in the column of the observed category and 0 in the others. Never code it as the integers 1, \dots, K: a linear model then treats the category as a quantity, which imposes an order and equal spacing between consecutive codes. The exception is a genuinely ordinal scale, such as a severity from minor to critical, whose order is real; even then, integer codes assume equal steps between the levels, a modelling choice to check rather than a fact.

One-hot codes also explain learned representations. A one-hot row vector \mathbf{e}_k^\top times a trainable K \times m table picks out row k of the table, so an embedding is such a table learned with the rest of the model, and it works the same way for any categorical feature (a part number, a supplier, a sensor ID) as for the tokens of a language model. Module 06, Section 4 shows the token embedding entering the transformer’s residual stream, and Module 06, Section 11 counts its V \times d table.

The running examples

Four examples recur in this module, so that each new idea lands on something already familiar:

• a synthetic linear regression with N = 200 rows and d = 3 features (Sections 2–4, Lab 1);
• noisy samples of \sin(2\pi x) fitted by polynomials (Sections 8–9, Lab 3);
• a classifier of breast tumours as malignant or benign, on the dataset that ships with scikit-learn (Sections 6–7, Labs 2 and 4);
• a hyperelastic material model fitted to a compression test of a soft hydrogel (Sections 5 and 12, Lab 5).

Linear regression is the simplest model that exercises every ingredient, and the only one in this module whose every property can be computed exactly. This section fits it in closed form, reads the answer as a projection, and says when the problem is ill-posed and how to compute the answer stably.

Set-up

Let f(\mathbf{x}) = \mathbf{w}^\top\mathbf{x} + b. Absorb the intercept b by appending a constant 1 to every input, so that f(\mathbf{x}) = \mathbf{w}^\top\mathbf{x} with \mathbf{w} \in \R^{d+1} and b its last entry. Stack the inputs as the rows of \mathbf{X} \in \R^{N\times(d+1)}, whose last column is all ones, and the outputs as \mathbf{y} \in \R^N. The loss is the mean squared error:

Section 5 shows why squared error is the right loss when the noise is Gaussian; here take it as given.

The gradient, derived

Expand the squared norm as an inner product:

The two cross terms \mathbf{w}^\top\mathbf{X}^\top\mathbf{y} and \mathbf{y}^\top\mathbf{X}\mathbf{w} are the same scalar, one the transpose of the other, which is where the 2 comes from. Two identities differentiate the pieces. For a symmetric matrix \mathbf{A}, \nabla_{\mathbf{w}}(\mathbf{w}^\top\mathbf{A}\mathbf{w}) = 2\mathbf{A}\mathbf{w}: writing the sum out, \partial_{w_k}\sum_{j,l}w_jA_{jl}w_l = \sum_l A_{kl}w_l + \sum_j w_jA_{jk} = 2(\mathbf{A}\mathbf{w})_k, the last step by symmetry. For a constant vector \mathbf{c}, \nabla_{\mathbf{w}}(\mathbf{c}^\top\mathbf{w}) = \mathbf{c}. With \mathbf{A} = \mathbf{X}^\top\mathbf{X} and \mathbf{c} = \mathbf{X}^\top\mathbf{y},

The gradient is the residual vector \mathbf{X}\mathbf{w} - \mathbf{y} mapped back to parameter space by \mathbf{X}^\top. Differentiating once more gives the Hessian \frac{2}{N}\mathbf{X}^\top\mathbf{X}, the same at every \mathbf{w}. It is positive semidefinite, because \mathbf{v}^\top\mathbf{X}^\top\mathbf{X}\mathbf{v} = \|\mathbf{X}\mathbf{v}\|^2 \ge 0 for every \mathbf{v}. So \mathcal{L} is a convex quadratic, and every point where its gradient vanishes is a global minimum.

The normal equations

Setting the gradient to zero gives the normal equations

If \mathbf{X} has full column rank, its d + 1 columns linearly independent, then \|\mathbf{X}\mathbf{v}\|^2 > 0 for every \mathbf{v} \ne \mathbf{0}, so \mathbf{X}^\top\mathbf{X} is positive definite and invertible, and the minimiser is unique:

Full column rank needs at least d + 1 linearly independent rows: no fewer examples than parameters. The formula is for reading; how to compute \mathbf{w}^\ast comes below.

The geometry: an orthogonal projection

As \mathbf{w} varies, \mathbf{X}\mathbf{w}, a linear combination of the columns of \mathbf{X}, sweeps out the column space of \mathbf{X}, a subspace of dimension d + 1 inside \R^N. Least squares picks the point of that subspace closest to \mathbf{y}. With the residual \mathbf{r} = \mathbf{y} - \mathbf{X}\mathbf{w}^\ast, the normal equations read

each column of \mathbf{X} has zero inner product with the residual. The residual is perpendicular, or normal, to the column space, which is where the equations get their name, and \hat{\mathbf{y}} = \mathbf{X}\mathbf{w}^\ast is the orthogonal projection of \mathbf{y} onto it: the foot of the perpendicular from \mathbf{y}, which is the closest point (Figure 1.1).

The projection is linear: \hat{\mathbf{y}} = \mathbf{P}\mathbf{y} with

the hat matrix, so called because it puts the hat on \mathbf{y}. (Statistics texts write it \mathbf{H}; this module keeps \mathbf{H} for the Hessian.) It is symmetric and idempotent, \mathbf{P}^2 = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}(\mathbf{X}^\top\mathbf{X})(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top = \mathbf{P}: projecting twice changes nothing. Its trace counts the dimensions it projects onto. Using \operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{B}\mathbf{A}), \operatorname{tr}\mathbf{P} = \operatorname{tr}\big((\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{X}\big) = \operatorname{tr}\mathbf{I}_{d+1} = d + 1, the number of fitted parameters. Section 9 turns this count into the effective degrees of freedom of ridge regression.

Two consequences follow when \mathbf{X} contains the column of ones. That column’s row of \mathbf{X}^\top\mathbf{r} = \mathbf{0} reads \sum_i r_i = 0: the residuals sum to zero, and \hat{\mathbf{y}} has the same mean \bar y as \mathbf{y}. And because \bar y\mathbf{1} also lies in the column space, \mathbf{y} - \bar y\mathbf{1} = (\hat{\mathbf{y}} - \bar y\mathbf{1}) + \mathbf{r} splits the centred data into a part inside the column space and a part perpendicular to it. Pythagoras gives

the total sum of squares is the explained plus the residual sum of squares. It defines the coefficient of determination

the fraction of the variance of \mathbf{y} about its mean that the fit explains; on the training data it lies between 0 and 1. Predicting the mean everywhere is the model with only an intercept, so R^2 compares the fit with that baseline. Section 7 uses R^2 on held-out data, where the identity no longer holds and R^2 can be negative.

When the columns are dependent

If one column of \mathbf{X} is a linear combination of the others, there is a \mathbf{v} \ne \mathbf{0} with \mathbf{X}\mathbf{v} = \mathbf{0}; then \mathbf{X}^\top\mathbf{X}\mathbf{v} = \mathbf{0}, and \mathbf{X}^\top\mathbf{X} is singular. Two ways it happens in practice: a temperature in °C next to the same temperature in °F, with an intercept, since T_{\text{F}} = 1.8\,T_{\text{C}} + 32 makes the °F column a combination of the °C column and the column of ones; and all K one-hot columns of a category next to an intercept, since the K columns sum to the column of ones. Every \mathbf{w}^\ast + t\mathbf{v} then fits equally well. The projection \hat{\mathbf{y}} is still unique, and so are the predictions at the training inputs and at any new input that obeys the same dependency; the weights are not. The pseudo-inverse solution \mathbf{X}^{+}\mathbf{y} picks, among all the minimisers, the one of smallest norm, and it is what np.linalg.lstsq returns. The ridge penalty of Section 9 also removes the ambiguity. Neither makes the individual weights of dependent columns meaningful.

Computing it

Never form (\mathbf{X}^\top\mathbf{X})^{-1}; solve. The condition number \kappa(\mathbf{X}), the ratio of the largest to the smallest singular value of \mathbf{X}, bounds how much rounding errors can be amplified: of the 16 significant digits of double precision, a solve can lose about \log_{10}\kappa. Three methods, in increasing order of robustness:

• Cholesky of \mathbf{X}^\top\mathbf{X}. Factor \mathbf{X}^\top\mathbf{X} = \mathbf{L}\mathbf{L}^\top with \mathbf{L} lower-triangular and solve two triangular systems. It is the fastest, but forming \mathbf{X}^\top\mathbf{X} squares the condition number, \kappa(\mathbf{X}^\top\mathbf{X}) = \kappa(\mathbf{X})^2, because the eigenvalues of \mathbf{X}^\top\mathbf{X} are the squared singular values of \mathbf{X}. At \kappa(\mathbf{X}) = 10^6 that is up to 12 lost digits instead of 6.
• QR of \mathbf{X}. Write \mathbf{X} = \mathbf{Q}\mathbf{R}, with orthonormal columns in \mathbf{Q} and \mathbf{R} upper-triangular. Since \mathbf{Q}^\top\mathbf{Q} = \mathbf{I}, the normal equations become \mathbf{R}^\top\mathbf{R}\mathbf{w} = \mathbf{R}^\top\mathbf{Q}^\top\mathbf{y}, that is \mathbf{R}\mathbf{w} = \mathbf{Q}^\top\mathbf{y}, solved by back-substitution. The method works with \mathbf{X} itself and never squares \kappa.
• SVD. \mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top gives \mathbf{w} = \mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{U}^\top\mathbf{y}. The singular values on the diagonal of \boldsymbol{\Sigma} display the conditioning, and inverting only the non-zero ones gives the minimum-norm solution of the rank-deficient case.

np.linalg.lstsq uses an SVD-based LAPACK driver and is the default to reach for. Each method costs about Nd^2 operations to form \mathbf{X}^\top\mathbf{X} or to factor \mathbf{X}, plus about d^3 for the small solve that follows.

The code, and the irreducible error

Here is all of it in NumPy, on the data of Lab 1: 200 inputs with three standard-normal features, true weights (1.5, -2.0, 0.5), intercept 0.7 and Gaussian noise of standard deviation 0.1. Gradient descent, the subject of the next section, is included for comparison.

import numpy as np

rng = np.random.default_rng(0)
N, d = 200, 3
X = rng.normal(size=(N, d))
w_true = np.array([1.5, -2.0, 0.5]); b_true = 0.7
y = X @ w_true + b_true + 0.1 * rng.normal(size=N)   # noise: the part no model recovers

Xb = np.hstack([X, np.ones((N, 1))])                 # absorb the bias

# closed form
w_closed = np.linalg.solve(Xb.T @ Xb, Xb.T @ y)

# gradient descent
w = np.zeros(d + 1); eta = 0.1
for step in range(500):
    grad = (2 / N) * Xb.T @ (Xb @ w - y)
    w -= eta * grad
print(np.round(w_closed, 3), np.round(w, 3))         # both near [1.5, -2.0, 0.5, 0.7]

# the robust default: an SVD-based solver that never forms Xb.T @ Xb
w_lstsq = np.linalg.lstsq(Xb, y, rcond=None)[0]
rss = np.sum((y - Xb @ w_lstsq) ** 2)               # residual sum of squares
print(np.round(w_lstsq, 3))
print(f"residual RMS {np.sqrt(rss / N):.4f}")
print(f"sqrt(RSS / (N - 4)) {np.sqrt(rss / (N - 4)):.4f}")

[ 1.491 -2.01   0.505  0.696] [ 1.491 -2.01   0.505  0.696]
[ 1.491 -2.01   0.505  0.696]
residual RMS 0.1000
sqrt(RSS / (N - 4)) 0.1010

The term 0.1\cdot\mathcal{N}(0, 1) in the data is the irreducible error: no function of \mathbf{x} predicts it, so no model can bring the expected squared error on new data below its variance, 0.01. A residual RMS equal to the noise level says that nothing recoverable is left. A model whose training residuals fall well below the noise level is fitting the noise, which new data will not share. That is the first appearance of the central problem of Section 8.

The normal equations solve least squares in one step. Gradient descent solves it by iteration,

with a step size \eta > 0 called the learning rate. There are three reasons to iterate when a closed form exists. Cost: the solve takes about Nd^2 + d^3 operations, a gradient step about Nd, two matrix–vector products. Memory: the solve needs all the data at once, while the stochastic steps of Section 4 need only a few rows. Generality: no model from Module 02 on has a closed form, and gradient descent is what trains them all. Least squares is the one case in which its behaviour can be derived exactly, which is why it is derived here. The rule of thumb: for small d, call lstsq; for large N or d, or for data that arrive as a stream, use the stochastic gradient descent of Section 4; for an ill-conditioned problem, rescale the features first (below).

The quadratic, exactly

Write \mathbf{H} = \frac{2}{N}\mathbf{X}^\top\mathbf{X} for the Hessian of Section 2. Put \mathbf{w} = \mathbf{w}^\ast + \mathbf{e} in the loss and expand:

The middle term vanishes by the normal equations. Dividing by N,

and its gradient is \nabla\mathcal{L}(\mathbf{w}) = \mathbf{H}(\mathbf{w} - \mathbf{w}^\ast). The loss is its minimum plus a bowl whose shape is \mathbf{H}. Follow the error \mathbf{e}_t = \mathbf{w}_t - \mathbf{w}^\ast: subtracting \mathbf{w}^\ast from both sides of the update gives

Because \mathbf{H} is symmetric it has an orthonormal basis of eigenvectors, \mathbf{H} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top, with eigenvalues \lambda_{\min} = \lambda_1 \le \dots \le \lambda_{d+1} = \lambda_{\max}, all non-negative. In the rotated coordinates \tilde{\mathbf{e}} = \mathbf{Q}^\top\mathbf{e} the matrix \mathbf{I} - \eta\mathbf{H} becomes the diagonal \mathbf{I} - \eta\boldsymbol{\Lambda}, and each component evolves alone:

That one line contains everything there is to know about the learning rate on a quadratic.

Stability

Every component shrinks, from every start, only if |1 - \eta\lambda_i| < 1 for every i, that is 0 < \eta\lambda_i < 2. The steepest direction binds first:

Follow the factor 1 - \eta\lambda_{\max} of the steepest direction as \eta grows (Figure 1.2 shows the first, third and fifth regimes on one problem):

• \eta < 1/\lambda_{\max}: every factor lies between 0 and 1, and every component shrinks monotonically, without changing sign.
• \eta = 1/\lambda_{\max}: the steepest factor is 0, and that component is solved in one step.
• 1/\lambda_{\max} < \eta < 2/\lambda_{\max}: the steepest factor is negative, between −1 and 0. That component changes sign at every step, so the iterate zig-zags across the valley while it converges.
• \eta = 2/\lambda_{\max}: the factor is −1, and that component flips sign for ever, neither growing nor decaying.
• \eta > 2/\lambda_{\max}: the factor is below −1, and that component grows geometrically; the loss explodes.

The boundary is sharp. Lab 1 runs 200 steps at 0.99 and at 1.01 times 2/\lambda_{\max} and ends with a training MSE of 0.0104 in the first case and 3.8\times10^3 in the second.

Speed

Stability is set by the steepest direction; speed is set by the flattest. The error norm shrinks per step by at most \rho(\eta) = \max_i|1 - \eta\lambda_i|, the largest factor in absolute value (the rotation \mathbf{Q} preserves norms). Since |1 - \eta\lambda| is a V in \lambda, its largest value over the eigenvalues is at one of the two ends: \rho = \max(|1 - \eta\lambda_{\min}|, |1 - \eta\lambda_{\max}|). Raising \eta shrinks the first term and, beyond 1/\lambda_{\max}, grows the second. The best fixed step makes them equal:

where \kappa = \lambda_{\max}/\lambda_{\min} is the condition number of \mathbf{H} (for \mathbf{H} \propto \mathbf{X}^\top\mathbf{X} it is \kappa(\mathbf{X})^2 in the notation of Section 2). Reducing the error by a factor \epsilon takes \rho^t \le \epsilon, that is

using -\ln\rho^\ast = \ln(1 + 1/\kappa) - \ln(1 - 1/\kappa) \approx 2/\kappa. By (3.1) the excess loss is quadratic in the error, so it shrinks as \rho^{2t}. Even at the best fixed step, the number of steps grows in proportion to \kappa: that is the cost of ill-conditioning. Module 02, Section 7 shows how momentum reduces it to about \sqrt{\kappa}.

Where the conditioning comes from

The condition number is a property of the features, and it can be read off them.

Units. If the features are centred and uncorrelated, with standard deviations s_j, then \frac{1}{N}\mathbf{X}^\top\mathbf{X} is diagonal: s_j^2 for each feature and 1 for the intercept, whose column of ones is orthogonal to every centred column. So \kappa = (s_{\max}/s_{\min})^2 (provided the intercept’s 1 lies between the extremes). Take a span length in millimetres with a standard deviation of 1000 mm next to an elastic modulus in gigapascals with a standard deviation of 0.001 GPa, which is 1 MPa of batch-to-batch scatter. The standard deviations differ by a factor of 10^6, so \kappa = 10^{12}, and even at the best fixed step gradient descent needs about (\kappa/2)\ln(10^6) \approx 6.9\times10^{12} iterations to gain six digits. After standardisation \kappa \approx 1.

Offsets. A feature that is not centred couples to the intercept. One feature with mean 100 and standard deviation 1, next to the column of ones, gives

because the off-diagonal entry is the feature’s mean and the corner is its mean square, 100^2 + 1^2. The trace is 10002 and the determinant 10001 - 10000 = 1, so the eigenvalues are about 10002 and 1/10002 = 1.0\times10^{-4}, and \kappa \approx 1.0\times10^{8}, from a feature of unremarkable scale. The loss hardly changes when the weight and the intercept move together, one up and the other down by 100 times as much, which is a long, flat valley. Centring alone fixes it: subtract the mean and the matrix becomes the identity.

Correlation. Two standardised features with correlation r give the correlation matrix \begin{bmatrix} 1 & r \\ r & 1 \end{bmatrix}, whose eigenvalues are 1 + r along (1, 1) and 1 - r along (1, -1), so \kappa = (1 + r)/(1 - r): 19 at r = 0.9 and 199 at r = 0.99. Standardisation fixes units; it does not fix correlation.

The remedy for units and offsets is to standardise each feature, x'_j = (x_j - \mu_j)/s_j, with the mean \mu_j and standard deviation s_j computed on the training set only and applied unchanged to validation data, test data and inputs in service. Computing them on all the data lets the held-out rows shape the fit: leakage in miniature, which Section 10 treats in general. The standardised model is the same model in new coordinates, so its weights map back exactly:

Beyond quadratics

Near any smooth minimum a loss is approximately quadratic, with \mathbf{H} its Hessian at the minimum, so \eta < 2/\lambda_{\max}(\mathbf{H}) also governs local stability for the networks of Module 02. For logistic regression (Section 6) the curvature changes with \mathbf{w}, and a bound computed from its worst case is sufficient but not necessary: in Lab 2 the loss still falls at more than eight times that bound.

The gradient of the empirical risk is an average over all N examples, so computing it exactly costs a pass over the data per step. Stochastic gradient descent (SGD) replaces it with the average over a random mini-batch \mathcal{B} of B examples:

where \ell_i(\theta) = \ell(f_\theta(\mathbf{x}_i), y_i). B = N is the full-batch gradient descent of Section 3; B = 1 is SGD in its original form.

Unbiased, with variance falling as 1/B

Draw the B indices independently and uniformly from 1, \dots, N (with replacement). Each drawn gradient \nabla\ell_{i_k} is then a random vector with mean \frac{1}{N}\sum_i\nabla\ell_i = \nabla\mathcal{L} and covariance

the spread of the per-example gradients. Linearity of expectation gives \E[\mathbf{g}_{\mathcal{B}}] = \nabla\mathcal{L}: the mini-batch gradient is unbiased. Independence makes the covariances of the B terms add, while the factor 1/B enters squared:

A batch four times larger halves the noise’s standard deviation at four times the cost. Sampling without replacement, as a shuffled pass through the data does, multiplies the covariance by (N - B)/(N - 1): slightly less, and exactly zero at B = N, where the batch is the whole dataset.

One pass over the data, an epoch, is \lceil N/B\rceil updates instead of one, each costing about Bd operations for a linear model instead of Nd. Far from the minimum the per-example gradients largely agree, so a mini-batch step makes nearly the progress of a full step at a fraction of the cost; that is why SGD wins on large data. In practice: reshuffle every epoch, because data stored in order (by class, by time, by machine) would otherwise give batches that are not samples of the whole; fix the random seeds so that runs can be repeated; and choose B between 32 and a few thousand. Language models train on batches of millions of tokens, for reasons given in Module 08, Section 6.

The loop that does it, continuing the code of Section 2:

rng = np.random.default_rng(1)
w = np.zeros(d + 1); eta, B = 0.1, 32
for epoch in range(50):
    order = rng.permutation(N)                       # reshuffle every epoch
    for k in range(0, N, B):
        idx = order[k:k + B]                         # the last batch holds 8 rows
        grad = (2 / len(idx)) * Xb[idx].T @ (Xb[idx] @ w - y[idx])
        w -= eta * grad
excess = np.mean((Xb @ w - y) ** 2) - np.mean((Xb @ w_lstsq - y) ** 2)
print(np.round(w, 3), f"excess loss {excess:.1e}")

[ 1.489 -2.016  0.507  0.688] excess loss 1.3e-04

After 350 updates the weights are close to the least-squares solution but not on it, and more epochs at the same step would not bring them closer. The reason is the subject of the rest of this section.

The noise floor

At the minimum the full gradient is zero, but the per-example gradients are not, so the mini-batch gradient is not zero either and the iterate keeps moving. Take a one-dimensional quadratic, \mathcal{L}(\theta) = \mathcal{L}(\theta^\ast) + \frac{1}{2}\lambda(\theta - \theta^\ast)^2, and model the mini-batch gradient as the true gradient \lambda(\theta - \theta^\ast) plus noise \xi_t, drawn afresh at each step, with mean zero and variance s^2/B, where s^2 is the variance of the per-example gradients. The update gives

Square and take expectations. The cross term vanishes because \xi_t has mean zero and is independent of \theta_t, which earlier batches decided. With V_t = \E[(\theta_t - \theta^\ast)^2],

For a stable step, |1 - \eta\lambda| < 1, the distance V_t - V to the fixed point shrinks by (1 - \eta\lambda)^2 per step, so V_t converges to the V that solves V = (1 - \eta\lambda)^2V + \eta^2s^2/B. Using 1 - (1 - \eta\lambda)^2 = \eta\lambda(2 - \eta\lambda):

A constant step does not converge. The iterate hovers around \theta^\ast with variance V, and the expected excess loss is \frac{1}{2}\lambda V = \eta s^2/\big(2B(2 - \eta\lambda)\big): a noise floor proportional to \eta/B.

Getting below the floor

There are two ways down. Enlarge the batch, which divides the floor by B at a proportional cost per update; or let the step decay. Robbins and Monro (1951) gave the conditions on a schedule \eta_t under which the iterate converges:

The first lets the steps add up to any distance, so the iterate can reach the minimum from any start. The second bounds the noise the steps inject, whose variance enters as \eta_t^2s^2/B in the recursion above. The schedule \eta_t = \eta_0/(1 + t/\tau) meets both: it stays near \eta_0 for the first \tau or so updates and then behaves like \eta_0\tau/t, whose sum diverges like \ln t while the sum of its squares converges.

Return to the widget of Section 3 and raise its gradient-noise slider \sigma_g: the loss falls as before, then flattens at the dashed line of the predicted floor. With \kappa = 20 and \sigma_g = 0.3 the floor is 0.0046 at \eta = 0.1, and halving \eta to 0.05 halves it, to 0.0023. At large steps the factor 2 - \eta\lambda also matters: the floor is 0.0264 at \eta = 0.5 and 0.0681 at \eta = 1.0.

Choosing the learning rate

When the curvature is known, as for least squares, start from 2/\lambda_{\max} and stay a safe factor below it. Otherwise measure: run a short sweep of learning rates spaced by a factor of about 3 on a log scale, from 10^{-4} to 1 (0.0001, 0.0003, 0.001, ..., 0.3, 1), for a few hundred updates each, take the largest \eta whose loss falls smoothly, and decay it over the run. Then read the loss curve on a log scale:

• rising, or NaN within a few steps: \eta is beyond the stability limit; divide it by 3 to 10;
• a straight but shallow descent: \eta is too small, or the problem is ill-conditioned; compute \kappa before raising \eta;
• a fall that turns into a noisy plateau: the noise floor; decay the step or enlarge the batch, and do not read the plateau as convergence.

Momentum and Adam are in Module 02, Sections 7 and 8, and the learning-rate range test and schedules in Section 9 of that module.

Noise, flat minima and non-convex losses

The noise has a second effect, which may be beneficial: small batches tend to find flatter minima, and flatter minima tend to generalise better. Keskar et al. (2017) observed that large-batch training of networks converges to sharper minima that generalise worse. This is an empirical observation with a partial theory, not a law, and the link between flatness and generalisation is still debated.

For non-convex models, which means every model from Module 02 on, gradient descent finds a local minimum or a saddle point, with no guarantee that it is the best one. In practice, for large networks, the minima found are usually good enough, and the difficulty lies elsewhere: in conditioning, initialisation and the learning rate, which Module 02 takes up.

Sections 2–4 minimised squared error without asking why the errors should be squared rather than, say, taken in absolute value. This section answers the question, and the answer turns a choice that looks arbitrary into an assumption you can check: a loss is the negative log-likelihood of the data under a model of the noise. Every loss in this module, and the cross-entropy that trains every language model in the series, comes from that one principle. It needs a little probability first.

A probability refresher

A random variable is a quantity whose value is uncertain until it is observed: the next reading of a strain gauge, the label of the next weld inspected. A discrete random variable has a probability mass function p(y) = P(Y = y), whose values lie between 0 and 1 and sum to

1. A continuous one has a probability density p(y), and probabilities are areas under it, P(a \le Y \le b) = \int_a^b p(y)\,dy. A density is not a probability. It can exceed 1 wherever the distribution is concentrated; only its integral must equal 1.

Four distributions do most of the work in this series:

• Gaussian: \mathcal{N}(y;\mu,\sigma^2) = (2\pi\sigma^2)^{-1/2}\exp\big(-(y-\mu)^2/(2\sigma^2)\big), with mean \mu and variance \sigma^2. Its peak density 1/(\sqrt{2\pi}\,\sigma) exceeds 1 whenever \sigma < 0.399.
• Laplace: (2b)^{-1}\exp(-\lvert y-\mu\rvert/b), with mean \mu and variance 2b^2. Against a Gaussian of the same variance it is more sharply peaked and has heavier tails: at y = \mu its density is 1/\sqrt2 = 0.707 against the Gaussian’s 0.399 at unit variance.
• Bernoulli, for y \in \{0, 1\}: p^y(1-p)^{1-y}, which is p when y = 1 and 1 - p when y = 0.
• Categorical, for y \in \{1, \dots, K\}: \prod_k p_k^{[y=k]}, where [y = k] is 1 when y = k and 0 otherwise, so the product picks out p_y.

The expectation \E[X] is the probability-weighted average of X (a sum over a mass function, an integral over a density), and the variance is \operatorname{Var}(X) = \E[(X - \E[X])^2]. Three rules follow from the definitions. For constants a and c, \E[aX + c] = a\E[X] + c and \operatorname{Var}(aX + c) = a^2\operatorname{Var}(X); for independent X and Y, \operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y). Section 8 uses them to compute the variance of a fitted curve, and Section 12 the uncertainty of a fitted material constant.

Two random variables are independent when their joint probability factorises, p(y_1, y_2) = p(y_1)\,p(y_2). The standard assumption about a dataset is that its examples are independent and identically distributed (i.i.d.): each is drawn separately from the same P. The probability of the whole dataset is then a product over its examples. Engineering data break the assumption often. Ten measurements of one specimen share that specimen’s quirks; consecutive readings of one sensor share its drift. Treating them as independent overstates how much evidence they hold, and splitting them at random between training and test leaks information across the split (group and time leakage, Section 10).

Likelihood, and why we take its logarithm

Let a model assign a probability, or a density, p(y \mid \mathbf{x};\theta) to every possible output. For i.i.d. data the likelihood is

read as a function of \theta with the data held fixed. It is not a probability distribution over \theta: it need not integrate to 1 over \theta, and it says how well each \theta explains the data, not how probable \theta is. Maximum likelihood estimation chooses the \theta under which the observed data are most probable.

In practice one maximises the log-likelihood \sum_i \ln p(y_i \mid \mathbf{x}_i;\theta) instead, for two reasons. The logarithm is increasing, so it has the same maximiser. And products of many probabilities underflow: a thousand probabilities of 0.1 multiply to 10^{-1000}, far below the smallest normal double-precision number (about 2.2\times10^{-308}), so the product is stored as 0.0 and every \theta looks equally bad. The sum of the logarithms, 1000\ln 0.1 = -2302.6, is an ordinary number.

import numpy as np

probs = np.full(1000, 0.1)                 # a thousand events, each of probability 0.1
print(np.prod(probs))                      # the likelihood itself underflows
print(f"{np.log(probs).sum():.1f}")        # the log-likelihood is an ordinary number

0.0
-2302.6

Gaussian noise gives squared error

Suppose the data were generated as y = f_\theta(\mathbf{x}) + \varepsilon, with \varepsilon \sim \mathcal{N}(0, \sigma^2) independent across examples. Then y given \mathbf{x} is Gaussian with mean f_\theta(\mathbf{x}), and the log-likelihood is

The first term does not depend on \theta. The second is the sum of squared errors times the negative constant -1/(2\sigma^2). Maximising the log-likelihood over \theta is therefore minimising \sum_i (y_i - f_\theta(\mathbf{x}_i))^2, which is N times the mean squared error of Section 2. The value of \sigma scales and shifts the objective but does not move its maximiser: least squares is the maximum-likelihood fit under Gaussian noise of any width.

The noise level itself can be estimated by maximum likelihood. Write v = \sigma^2 and r_i = y_i - f_{\hat\theta}(\mathbf{x}_i) for the residuals of the fitted model, and set the derivative with respect to v to zero:

This estimate is biased low, for the same reason as the training loss of Section 1: the residuals were made small by the fit. For a linear model with p parameters, \E[\sum_i r_i^2] = (N - p)\sigma^2, a standard result (Hastie et al. 2009, Section 3.2) with a one-line proof from Section 2’s hat matrix \mathbf{P}. The residual is \mathbf{r} = (\mathbf{I} - \mathbf{P})\mathbf{y} = (\mathbf{I} - \mathbf{P})\boldsymbol{\varepsilon}, because (\mathbf{I} - \mathbf{P})\mathbf{X} = \mathbf{0}, so

(The middle step uses that \mathbf{I} - \mathbf{P} is symmetric and idempotent, so \lVert(\mathbf{I} - \mathbf{P})\boldsymbol{\varepsilon}\rVert^2 = \boldsymbol{\varepsilon}^\top(\mathbf{I} - \mathbf{P})\boldsymbol{\varepsilon}, and that \E[\boldsymbol{\varepsilon}^\top\mathbf{A}\boldsymbol{\varepsilon}] = \sigma^2\operatorname{tr}\mathbf{A} for independent noise of variance \sigma^2.)

The residual lives in the N - p dimensions the model cannot reach, and only noise in those dimensions shows in it. The unbiased estimate therefore divides by N - p. That is the difference between Section 2’s residual RMS of 0.1000, which is the maximum-likelihood \hat\sigma, and its \sqrt{\text{RSS}/(N-4)} = 0.1010.

Laplace noise gives absolute error

If instead \varepsilon has the Laplace density, the negative log-likelihood of one example is \ln(2b) + \lvert y - f_\theta(\mathbf{x})\rvert/b, and maximum likelihood minimises \sum_i \lvert y_i - f_\theta(\mathbf{x}_i)\rvert, the absolute error. The two losses differ most clearly for the simplest model, a constant prediction c. Under squared error the derivative of \sum_i (y_i - c)^2 is -2\sum_i (y_i - c), which is zero at the mean. Under absolute error the derivative of \sum_i \lvert y_i - c\rvert is -\sum_i \operatorname{sign}(y_i - c), which is zero when as many points lie above c as below it: at the median. The mean responds to every point in proportion to its distance; the median only counts how many points lie on each side. That is why an absolute-error fit is robust to outliers and a squared-error fit is not.

Huber’s loss, quadratic for small residuals and linear beyond a threshold, is the usual compromise between the two; Module 02, Section 12 uses it.

Different noise per point: weighted least squares

Measurements often come with their own uncertainties: a load cell is less precise at the bottom of its range, a reading taken during a disturbance is flagged as poor. If point i has known noise standard deviation \sigma_i, the Gaussian negative log-likelihood is

The first sum does not involve \theta, so maximum likelihood minimises \sum_i w_i (y_i - f_\theta(\mathbf{x}_i))^2 with weights w_i = 1/\sigma_i^2. This is weighted least squares: a precise point pulls hard on the fit, an uncertain one gently. It is the likelihood of the hydrogel fit in Section 12, where every stress reading carries its own uncertainty.

The recipe

Read a row across: assuming Gaussian noise, fitting by squared error and summarising a sample by its mean are one decision, not three. The last two rows are the classifiers of Section 6; their best constant is found the same way, by setting a derivative to zero.

The average negative log-likelihood has a name that recurs throughout the series. For a data distribution p_\text{data} and a model distribution p_\theta, adding and subtracting \ln p_\text{data} inside the expectation gives

The left side is the cross-entropy H(p_\text{data}, p_\theta). It splits into the entropy of the data, which no choice of \theta can change, and the Kullback–Leibler divergence, which is never negative and is zero only when p_\theta = p_\text{data}. Maximising likelihood therefore moves the model’s distribution towards the data’s. Module 07, Section 1 reports a language model’s quality as its perplexity, which is \exp of exactly this cross-entropy, per token, on held-out text.

Adding a prior distribution over \theta turns maximum likelihood into maximum a posteriori (MAP) estimation, and the logarithm of the prior becomes a penalty on the parameters: Section 9 derives ridge and lasso regression this way.

Regression predicts a number; classification predicts one of K labels: defective or acceptable, benign or malignant, the next token of a text. The Bernoulli and categorical rows of Section 5’s table already give the loss. What is missing is a model whose output is a probability. This section builds the linear one, derives its gradient and curvature, and shows how it fails.

Probabilities from a linear score

For y \in \{0, 1\}, logistic regression models

with the intercept absorbed into \mathbf{w} as in Section 2. Here \sigma is the logistic (sigmoid) function, not a standard deviation. It maps any real score z, the logit, into (0, 1). Two properties are used below. First, \sigma(-z) = 1 - \sigma(z). Second, its derivative is

Solving p = \sigma(z) for z gives z = \ln\big(p/(1-p)\big), the logarithm of the odds. So logistic regression is linear in the log-odds: \ln(p/(1-p)) = \mathbf{w}^\top\mathbf{x}. Raising feature j by one unit adds w_j to the log-odds, which multiplies the odds by e^{w_j}. A coefficient of 0.7 doubles the odds per unit (e^{0.7} = 2.01). That is how a fitted coefficient is read.

The loss and its gradient

The Bernoulli negative log-likelihood of one example, with \hat p = \sigma(z), is the binary cross-entropy

Its derivative with respect to the logit follows from the chain rule in two factors:

The sigmoid’s derivative cancels both denominators, and what is left is predicted minus observed. Since z_i = \mathbf{w}^\top\mathbf{x}_i gives \partial z_i/\partial\mathbf{w} = \mathbf{x}_i, the gradient of the mean loss is

the same shape as the least-squares gradient \frac{2}{N}\mathbf{X}^\top(\mathbf{X}\mathbf{w} - \mathbf{y}), with \hat{\mathbf{p}} in place of \hat{\mathbf{y}} (the factor 2 belonged to the square).

Curvature, convexity and the step size

Differentiating once more, \partial\hat p_i/\partial\mathbf{w} = \hat p_i(1-\hat p_i)\mathbf{x}_i, so the Hessian is

For any \mathbf{v}, \mathbf{v}^\top\mathbf{H}\mathbf{v} = \frac{1}{N}\sum_i \hat p_i(1-\hat p_i)(\mathbf{x}_i^\top\mathbf{v})^2 \ge 0, so \mathbf{H} is positive semidefinite and the loss is convex: every minimum is global. But the gradient is nonlinear in \mathbf{w}, so setting it to zero has no closed-form solution. The options are gradient descent; Newton’s method, which here is called iteratively reweighted least squares because each Newton step solves a least-squares problem weighted by \mathbf{S}; or a quasi-Newton method. scikit-learn’s default solver is L-BFGS, a quasi-Newton method.

Section 3’s stability analysis carries over with one change. Since \hat p(1-\hat p) \le \tfrac14, with equality at \hat p = \tfrac12, the largest curvature satisfies \lambda_{\max}(\mathbf{H}) \le \lambda_{\max}(\mathbf{X}^\top\mathbf{X}/N)/4, and gradient descent is guaranteed stable for

This bound is sufficient, not necessary. The curvature reaches its worst case only where every \hat p_i is near \tfrac12; as the fit improves most \hat p_i move towards 0 or 1 and the curvature falls far below the bound. Lab 2 computes the bound for its data as \eta < 0.59, and in its extension the loss still falls at \eta = 5.

Why not squared error

Squared error through a sigmoid, \tfrac12(\hat p - y)^2, has gradient

with respect to the logit. The extra factor \hat p(1-\hat p) vanishes as \hat p approaches 0 or 1, and that includes the case where the model is confidently wrong, \hat p near 1 for y = 0. There, cross-entropy’s gradient \hat p - y is close to its largest value, while the squared-error gradient is close to zero: the model learns slowest exactly where it is most wrong (Figure 1.4, right). Squared error through a sigmoid is also not convex in \mathbf{w}: for y = 0 it rises and then flattens at 0.5, and gradient descent can stall on the plateau. A confidently wrong model must be able to learn fast. Exercise 4 puts numbers on the difference.

Computing the loss stably

The loss is better written in terms of the logit. Using \ln\sigma(z) = -\ln(1 + e^{-z}), \ln(1 - \sigma(z)) = \ln\sigma(-z) = -\ln(1 + e^{z}) and \ln(1 + e^{-z}) = \ln(1 + e^{z}) - z,

Evaluated naively, by computing \hat p = \sigma(z) first and then its logarithms, the loss fails at large \lvert z\rvert. The stable route computes \operatorname{softplus}(z) = \max(z, 0) + \ln(1 + e^{-\lvert z\rvert}), whose exponential never overflows; NumPy provides it as np.logaddexp(0, z). Module 02, Section 12 generalises this to the log-sum-exp of many logits.

import numpy as np

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

def bce_naive(z, y):
    p = sigmoid(z)
    return -(y * np.log(p) + (1 - y) * np.log(1 - p))

def bce_stable(z, y):
    # softplus(z) - y*z, with softplus(z) = ln(1 + e^z) computed by logaddexp
    return np.logaddexp(0.0, z) - y * z

with np.errstate(divide="ignore", invalid="ignore"):
    for z, y in [(40.0, 1), (-40.0, 1), (40.0, 0)]:
        print(f"z = {z:5.1f}, y = {y}:  naive {bce_naive(z, y):5.1f}"
              f"  stable {bce_stable(z, y):5.1f}")

z =  40.0, y = 1:  naive   nan  stable   0.0
z = -40.0, y = 1:  naive  40.0  stable  40.0
z =  40.0, y = 0:  naive   inf  stable  40.0

The decision boundary

Write the intercept separately again, z = \mathbf{w}^\top\mathbf{x} + b. A probability becomes a decision through a threshold t: predict 1 when \hat p \ge t. Because \sigma is increasing, \hat p \ge t is the same as z \ge \ln\big(t/(1-t)\big). At t = 0.5 the condition is \mathbf{w}^\top\mathbf{x} + b \ge 0, and the decision boundary \mathbf{w}^\top\mathbf{x} + b = 0 is a hyperplane with normal vector \mathbf{w}. The signed distance of a point from it is (\mathbf{w}^\top\mathbf{x} + b)/\lVert\mathbf{w}\rVert, so the logit is \lVert\mathbf{w}\rVert times that distance: \lVert\mathbf{w}\rVert sets how fast the probability changes as a point moves away from the boundary. Another threshold moves the boundary parallel to itself, to \mathbf{w}^\top\mathbf{x} + b = \ln\big(t/(1-t)\big).

The boundary is flat in whatever features the model is given. A curved boundary needs curved features: adding x_1^2, x_2^2 and x_1x_2 lets it be any conic, an ellipse among them. The concentric circles of Module 02, Section 1 need such features, or a network that learns them.

Separable data

Suppose some hyperplane separates the classes perfectly: z_i > 0 for every positive example and z_i < 0 for every negative one. Multiply \mathbf{w} and b by any c > 1. Every logit grows in size with its sign unchanged, every \hat p_i moves towards its label, and every example’s loss falls. The loss can always be lowered further, its infimum of 0 is approached only as \lVert\mathbf{w}\rVert \to \infty, and the maximum-likelihood estimate does not exist. Gradient descent obeys: \lVert\mathbf{w}\rVert grows without limit and the probabilities saturate at 0 and 1, claiming a certainty the data cannot support. A penalty on \lVert\mathbf{w}\rVert (Section 9) or early stopping is required. scikit-learn’s LogisticRegression includes an L2 penalty by default (C=1.0), which is why it never shows the problem; Lab 2 removes the penalty and watches \lVert\mathbf{w}\rVert climb from 3.2 after 100 steps to 51.2 after 100,000.

Softmax regression

For K classes the model computes K logits \mathbf{z} = \mathbf{W}\mathbf{x}, one row \mathbf{w}_k of \mathbf{W} per class, and turns them into probabilities with the softmax

Differentiate the loss with respect to one logit z_k. The second term gives e^{z_k}/\sum_j e^{z_j} = \hat p_k. The first gives -1 if k is the true class y and 0 otherwise. So

It is the binary pattern again, predicted minus observed, with the observed label one-hot encoded. Stacking the predicted probabilities as the rows of \hat{\mathbf{P}} (N\times K) and the one-hot labels as the rows of \mathbf{Y} gives the gradient for the whole weight matrix,

a K\times(d+1) matrix, the shape of \mathbf{W}.

Adding the same constant to every logit changes nothing, because e^{z_k + c}/\sum_j e^{z_j + c} = e^{z_k}/\sum_j e^{z_j}. One class’s weights can therefore be fixed at zero without loss. With K = 2,

which is logistic regression with weights \mathbf{w}_1 - \mathbf{w}_0. The predicted class is the one with the largest logit, and classes j and k tie where (\mathbf{w}_j - \mathbf{w}_k)^\top\mathbf{x} + (b_j - b_k) = 0, a hyperplane. Each decision region is an intersection of half-spaces, so the boundaries are piecewise linear (Figure 1.5, right).

One output layer for the whole series

Logits, a softmax and a cross-entropy loss form the output layer of every classifier in this series. A language model is such a classifier, whose classes are the entries of its vocabulary (tens to hundreds of thousands of them), applied at every position of a text to predict the next token (Module 07). The gradient \hat p - y, the overconfidence it produces on separable data and the need to compute logarithms stably all carry over unchanged.

The loss is what the optimiser can minimise: smooth, differentiable, averaged over examples. The metric is what the decision depends on: how many defective welds are missed, how many good ones are scrapped, how far a predicted sag is from the real one in millimetres. They are rarely the same quantity. A classifier can minimise cross-entropy well and serve its decision badly, when the decision needs high recall at a fixed false-positive rate. Choose the metric from the decision first; then choose the threshold, and if necessary the model, to serve it.

Every reported number carries its baseline. An RMSE of 0.3 mm on a sag prediction means nothing until you know that predicting the mean sag gives 1.1 mm and the closed-form solver gives 0.05 mm: the model is far better than knowing nothing and six times worse than the physics.

Regression metrics

For predictions \hat y_i of targets y_i on n held-out cases, with \bar y their mean:

RMSE is in the units of y and dominated by the largest errors. MAE is in the same units and, like the median of Section 5, robust to a few outliers. R^2 comes from Section 2’s Pythagoras: the fraction of the variance of y that the model explains, measured against the baseline of predicting the mean. On training data with an intercept it lies between 0 and 1. On held-out data nothing bounds it below: a model worse than predicting the mean has \text{SS}_\text{res} > \text{SS}_\text{tot} and a negative R^2. The sag model above has R^2 = 1 - (0.3/1.1)^2 = 0.93, which sounds excellent and hides the factor of six. Percentage errors such as \lvert y - \hat y\rvert/\lvert y\rvert are popular and break when y can be near zero, where a tiny error becomes an enormous percentage.

The confusion matrix

A binary classifier with a fixed threshold produces four kinds of outcome: true positives (TP), false positives (FP), false negatives (FN) and true negatives (TN). Every rate is a ratio of these counts:

• accuracy = (\text{TP} + \text{TN})/n;
• precision = \text{TP}/(\text{TP} + \text{FP}): of the cases flagged, the fraction that are positive;
• recall = \text{TP}/(\text{TP} + \text{FN}), also called sensitivity or the true-positive rate (TPR): of the positives, the fraction flagged;
• specificity = \text{TN}/(\text{TN} + \text{FP}), and the false-positive rate \text{FPR} = \text{FP}/(\text{FP} + \text{TN}) = 1 - \text{specificity};
• F1 = 2\cdot\text{precision}\cdot\text{recall}/(\text{precision} + \text{recall}), the harmonic mean of precision and recall.

A single rate hides which mistakes are being made, so state the confusion matrix whenever you state a rate. F1 is a harmonic rather than an arithmetic mean because the harmonic mean is dominated by the smaller value: a classifier cannot score well by buying one of precision or recall with the other.

Accuracy misleads when the classes are imbalanced. If 1% of cases are positive, the rule “always negative” scores 99% accuracy and finds nothing.

The base rate

Precision depends on how common positives are, and Bayes’ rule says exactly how. Let \pi be the prevalence, the fraction of cases that are positive. Then

The denominator is the probability of a flag, from a positive case or from a negative one. When \pi is small, the false flags \text{FPR}\cdot(1-\pi) dominate it even for a small FPR. A good detector of a rare event floods its users with false alarms.

Ranking: ROC and precision–recall curves

A classifier that outputs a score can be used at any threshold. Sweeping the threshold from high to low traces the ROC curve, TPR against FPR, from (0, 0), where nothing is flagged, to (1, 1), where everything is. The area under it, the AUC, has a direct meaning: it is the probability that a randomly chosen positive is scored above a randomly chosen negative (the Mann–Whitney statistic, with ties counted as one half). So it can be computed by counting pairs. TPR is computed among the positives and FPR among the negatives, so neither changes when the class proportions change, and neither does the AUC. That makes the AUC a good measure of how well a model ranks, and a misleading one when positives are rare: a 1% false-positive rate looks small on an ROC plot and can be most of the alarms.

When positives are rare, look at the precision–recall curve, precision against recall as the threshold sweeps. It shows what the person reading the alarms will experience. Its summary, the average precision, averages the precision over the recall levels at which positives are found. Its baseline is not one half: a classifier that scores at random has precision equal to the prevalence at every recall.

Choosing the threshold

The threshold is a decision, not a property of the model. Suppose the predicted probabilities are calibrated (defined below), a false positive costs C_\text{FP} and a false negative C_\text{FN}. For a case with probability p of being positive, flagging it costs (1 - p)\,C_\text{FP} in expectation, the cost incurred if it is negative, and passing it costs p\,C_\text{FN}. Flag it when p\,C_\text{FN} > (1 - p)\,C_\text{FP}, which rearranges to

The default 0.5 is right only when the two errors cost the same. The alternative is a constraint, such as recall at least 0.99 at the lowest false-positive rate that achieves it, read off the ROC curve on validation data. A model trained by minimising cross-entropy still needs one of these steps before it serves a decision that needs recall at a fixed false-positive rate.

Calibration

This section defines calibration for the whole series. A model is calibrated when its probabilities mean what they say: among the cases given \hat p = 0.8, 80% are positive. Formally, P(y = 1 \mid \hat p = p) = p for every p. Calibration matters whenever a probability is used as a probability: in the cost threshold above, in a risk estimate, or to decide when a model should abstain.

A reliability diagram groups the predictions into bins by \hat p and plots, for each bin, the observed fraction of positives against the mean predicted probability; a calibrated model lies on the diagonal (Figure 1.6 draws Lab 2’s). The expected calibration error summarises the diagram. With equal-width bins, bin b holding n_b of the N predictions, mean predicted probability \text{conf}_b and fraction of positives \text{freq}_b,

For K classes, or for a model’s answers to questions, the usual form bins the confidence of the predicted class, its largest probability, and compares it with the accuracy in the bin:

This is the top-label form of Guo et al. (2017); Modules 02, 07 and 09 use it and refer back here.

Three cautions. First, state the binning, 10 or 15 equal-width bins as usual or 5 for a test set of about a hundred, because ECE changes with it: Lab 2’s predictions give 0.023 with 5 bins and 0.039 with 10. Second, ECE is biased upward on small samples. Even a perfectly calibrated model’s bin frequencies scatter around their confidences, and the absolute value turns that scatter into a positive error. Third, ECE measures calibration, not discrimination: a model that predicts the base rate for every case is calibrated and useless. Report ECE together with the AUC and the Brier score \frac1N\sum_i(\hat p_i - y_i)^2, the mean squared error of the probabilities. The Brier score is a proper scoring rule: its expectation is smallest when the predicted probabilities are the true ones, so it rewards calibration and discrimination together.

import numpy as np

def ece_binary(p, y, n_bins=10):
    """Binary ECE of probabilities p against labels y in {0, 1}, equal-width bins."""
    edges = np.linspace(0.0, 1.0, n_bins + 1)
    bins = np.clip(np.digitize(p, edges[1:-1]), 0, n_bins - 1)
    ece = 0.0
    for b in range(n_bins):
        in_bin = bins == b
        if in_bin.any():                       # empty bins contribute nothing
            gap = abs(y[in_bin].mean() - p[in_bin].mean())
            ece += in_bin.mean() * gap         # weight n_b / N
    return ece

# the three-bin example below: 50 cases at 0.2 (5 positive), 30 at 0.5 (15), 20 at 0.9 (14)
p = np.repeat([0.2, 0.5, 0.9], [50, 30, 20])
y = np.concatenate([np.r_[np.ones(5), np.zeros(45)],
                    np.r_[np.ones(15), np.zeros(15)],
                    np.r_[np.ones(14), np.zeros(6)]])
print(f"ECE   {ece_binary(p, y):.4f}")
print(f"Brier {np.mean((p - y) ** 2):.4f}")

ECE   0.0900
Brier 0.1750

Models trained for long with cross-entropy tend to be overconfident, because the loss keeps rewarding larger logits on training examples that are already classified correctly. The remedy is to recalibrate on held-out data. Platt scaling fits a one-feature logistic regression from the model’s score to the label; isotonic regression fits a non-decreasing step function, which is more flexible and needs more data. For networks the standard method is temperature scaling, which divides all logits by one fitted constant (Module 02, Section 11, Module 07, Section 10). The direction of a miscalibration should be measured, not assumed: Lab 2’s L2-regularised model, with a mean confidence in its predicted class of 0.956 against an accuracy of 0.986, errs slightly the other way.

The training loss is biased low (Section 1), so it cannot say how a model will do on new inputs. This section makes the gap measurable: it sets up the data so that the gap can be seen, derives where it comes from, and gives the two tools that manage it in practice, validation curves and cross-validation.

Three sets

Split the data before doing anything else:

• the training set is what the optimiser sees;
• the validation set is used to choose between models and settings (a degree, a \lambda, a k);
• the test set is touched once, at the end, to report a number.

Every look at the test set that changes something, whether a feature, a setting or the choice of model, turns it into a validation set, and its number stops meaning what you will claim it means. Section 10 develops this discipline.

The capacity curve

The experiment of Lab 3: y = \sin(2\pi x) + \varepsilon with \varepsilon \sim \mathcal{N}(0, 0.3^2), N = 20 training points on the evenly spaced grid x_i = i/19, and polynomials of degree 0 to 15 fitted by least squares. The polynomials are written in the Legendre basis P_k(2x - 1) rather than as powers x^k. The fitted functions are the same, but the Legendre design matrix at degree 15 has condition number 47 against 6.7\times10^{5} for powers of 2x - 1, so the solve is accurate (Section 2). A validation set of 1,000 points at random x measures the error on new inputs.

Training RMSE never increases with the degree: each family of polynomials contains the previous one, so least squares can always do at least as well. Validation RMSE falls, then rises (Figure 1.7). On the left the model underfits: it cannot represent the signal, and both errors are high. On the right it overfits: it fits noise that the validation set does not share, and the gap between the two errors widens. The degree is one capacity knob. Others are the number of parameters, the depth of a tree, 1/k in k-nearest neighbours and 1/\lambda in ridge regression.

The bias–variance decomposition

Why validation error is U-shaped follows from an identity. Fix an input \mathbf{x}. A new observation there is y = f^\ast(\mathbf{x}) + \varepsilon, with \E[\varepsilon] = 0, \operatorname{Var}(\varepsilon) = \sigma^2, and \varepsilon independent of the training set \mathcal{D}. The fitted model depends on \mathcal{D}, which is random. Write \hat f = \hat f_{\mathcal{D}}(\mathbf{x}) for its prediction at \mathbf{x} and \bar f = \E_{\mathcal{D}}[\hat f] for the average prediction over training sets. Add and subtract \bar f and f^\ast:

Square it:

Take the expectation over \mathcal{D} and \varepsilon, term by term. \bar f - f^\ast is a constant, so \E[(\hat f - \bar f)(\bar f - f^\ast)] = (\bar f - f^\ast)\,\E_{\mathcal{D}}[\hat f - \bar f] = 0, because \E_{\mathcal{D}}[\hat f] = \bar f. Since \varepsilon is independent of \mathcal{D}, \E[(\hat f - \bar f)\varepsilon] = \E_{\mathcal{D}}[\hat f - \bar f]\,\E[\varepsilon] = 0. And \E[(\bar f - f^\ast)\varepsilon] = (\bar f - f^\ast)\E[\varepsilon] = 0. What remains is

Averaging over \mathbf{x} drawn from P gives the expected test error.

Bias is what the model family cannot represent, even on average over training sets. Variance is how much the fit changes from one training sample to another. Noise is irreducible: no model predicts \varepsilon. Simple models have high bias and low variance; flexible ones the reverse, and the validation curve is their sum. The decomposition is an exact identity for squared error. For the 0–1 loss of classification there is no clean additive version, only the same picture used loosely.

Learning curves

A learning curve plots training and validation error against the number of training points N for a fixed model. For least squares with p parameters, Section 5’s \E[\text{RSS}] = (N - p)\sigma^2 says the expected training MSE is bias² on the training points plus \sigma^2(1 - p/N), below the noise level, while the expected validation MSE is bias² + variance + \sigma^2, above it. As N grows the variance shrinks, roughly as \sigma^2 p/N, and both curves approach bias² + \sigma^2 from opposite sides.

The shape is a diagnosis. A large gap that persists means variance: more data will help. Both curves high and close together means bias: more data will not help, and more capacity or better features will. In Lab 3’s setting (Figure 1.8), at N = 20 degree 9 has expected training and validation MSE 0.045 and 0.132, a gap of 0.087, against 0.078 and 0.111 for degree 3. By N = 1{,}000 both pairs are within 0.001 of their limits, \sigma^2 + bias² = 0.0946 for degree 3 and 0.0900 for degree 9, and from about N = 115 on degree 9 is the better model. The best degree moves to the right as data accumulate.

k-fold cross-validation

A single validation set wastes data when data are scarce, and its verdict depends on which points it happened to get. k-fold cross-validation splits the data into k folds of nearly equal size, trains on k - 1 of them, validates on the remaining one, rotates through all k choices and averages the k validation errors. Every point is used for validation once and for training k - 1 times. The usual choices: k = 5 or 10, so each fit sees 80–90% of the data at the cost of k fits; stratified folds for classification, so each fold keeps the class proportions; and leave-one-out (k = N) for very small datasets.

The k fold errors also give a spread. Their mean comes with a standard error \text{sd}/\sqrt{k}, where sd is the standard deviation of the fold errors. It is approximate and optimistic: the training sets of any two folds overlap (with k = 5 they share three quarters of their points), so the fold errors are positively correlated and their mean varies more than \text{sd}/\sqrt{k} suggests. It is still the only honest way to say that one setting beats another by less than the noise.

Double descent

The U-shaped curve is the classical expectation, and for the models of this module it is the right one. It is not a law. When models are fitted all the way to interpolation, and the fitting rule picks the minimum-norm solution among the many that fit, test error can peak where the number of parameters p reaches the number of training points N, the interpolation threshold, and fall again beyond it (Figure 1.9). Belkin et al. (2019) named this double descent; Nakkiran et al. (2020) found it in deep networks, along the axes of model size and training time.

It does not contradict the decomposition, which is an identity. It changes how the variance behaves with size. At p = N exactly one function fits the data, and it must chase every noise value. Past the threshold infinitely many functions fit, and a rule that prefers small norms picks a smooth one, whose variance falls as p grows. Lab 3’s design shows it. The exact formulas, with np.linalg.lstsq supplying the minimum-norm fit beyond degree 19, give an expected test MSE of 0.96 at degree 15, about 16,000 at degree 19 (p = N = 20), 0.19 at degree 25 and 0.14 at degree 100: a second descent, but not below the classical minimum of 0.111 at degree 3. The practical lesson is to measure: do not assume a single U, and do not assume the second descent beats it. The guided reading for this session is the Belkin paper.

Section 8 traced the failure of flexible models to their variance. Regularisation trades some of that variance for bias on purpose, usually by adding a penalty to the loss: \mathcal{L}_\lambda(\theta) = \mathcal{L}(\theta) + \lambda\,\Omega(\theta), with \Omega = \lVert\mathbf{w}\rVert_2^2 for ridge regression and \lVert\mathbf{w}\rVert_1 for the lasso. The penalty is not an arbitrary add-on: it states what you believed about the parameters before seeing the data. In this section \lambda is the regularisation strength, and the eigenvalues of \mathbf{X}^\top\mathbf{X} are written s_i^2.

The prior is the regulariser

Bayes’ rule applied to the parameters reads

The prior p(\theta) says which parameters are plausible before the data, the likelihood of Section 5 how well each explains the data. The maximum a posteriori (MAP) estimate maximises the posterior’s logarithm:

Maximum likelihood is MAP with a flat prior. Any other prior adds -\ln p(\theta) to the negative log-likelihood: the prior is the regulariser.

Ridge as MAP, constants kept

Take Gaussian noise, y_i = \mathbf{w}^\top\mathbf{x}_i + \varepsilon_i with \varepsilon_i \sim \mathcal{N}(0, \sigma^2), and the prior \mathbf{w} \sim \mathcal{N}(\mathbf{0}, \tau^2\mathbf{I}): each weight is believed to be of size about \tau. The negative log-posterior is

Multiplying by the positive constant 2\sigma^2/N leaves the minimum where it is and gives this module’s mean loss plus a penalty:

A tighter prior (smaller \tau) or noisier data enlarge the penalty. More data shrink it: the evidence outweighs the prior, and with enough data MAP and maximum likelihood agree.

The closed form always exists

Set the gradient of (9.1) to zero, as in Section 2, and multiply by N/2:

\mathbf{X}^\top\mathbf{X} has eigenvalues s_i^2 \ge 0, the squared singular values of \mathbf{X}. Adding \lambda N\mathbf{I} adds \lambda N to each, so every eigenvalue is at least \lambda N > 0 and the system is solvable with a duplicated column, with more columns than rows, with any data.

Each direction shrunk by its own factor

With the thin singular value decomposition \mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top (orthonormal columns \mathbf{u}_i and \mathbf{v}_i, singular values s_i), \mathbf{X}^\top\mathbf{X} = \mathbf{V}\boldsymbol{\Sigma}^2\mathbf{V}^\top, the solution is \mathbf{w}_\lambda = \mathbf{V}(\boldsymbol{\Sigma}^2 + \lambda N\mathbf{I})^{-1}\boldsymbol{\Sigma}\mathbf{U}^\top\mathbf{y}, and the fitted values are

Least squares keeps every component \mathbf{u}_i^\top\mathbf{y} whole: Section 2’s projection. Ridge multiplies each by its own shrinkage factor s_i^2/(s_i^2 + \lambda N), near 1 where the data vary strongly and near 0 where they barely vary.

The variance says why. Along \mathbf{v}_i the least-squares coefficient is \mathbf{u}_i^\top\mathbf{y}/s_i; its noise part has variance \sigma^2/s_i^2, so a direction the data barely constrain is set by noise divided by a small number. The ridge coefficient, s_i\,\mathbf{u}_i^\top\mathbf{y}/(s_i^2 + \lambda N), has variance \sigma^2 s_i^2/(s_i^2 + \lambda N)^2, which peaks at s_i^2 = \lambda N and never exceeds \sigma^2/(4\lambda N). Ridge caps the variance in every direction and pays with bias only where the data had little to say.

The sum of the factors, \operatorname{df}(\lambda) = \sum_i s_i^2/(s_i^2 + \lambda N), is the effective degrees of freedom. It equals the number of columns at \lambda = 0 when they are independent (the trace of the hat matrix, Section 2) and falls smoothly to 0 as \lambda grows: a continuous capacity knob.

Weight decay

A gradient step on the ridge loss is

Every weight is first multiplied by 1 - 2\eta\lambda, slightly below 1, which is why the L2 penalty is called weight decay. AdamW’s decoupled version is in Module 02, Section 8.

Lasso: a Laplace prior, and exact zeros

A Laplace prior, p(w_j) = \frac{1}{2b}\exp(-\lvert w_j\rvert/b), peaked at zero with heavier tails, has negative logarithm \lVert\mathbf{w}\rVert_1/b plus a constant. The scaling of (9.1) gives the lasso, \frac{1}{N}\lVert\mathbf{X}\mathbf{w} - \mathbf{y}\rVert^2 + \lambda\lVert\mathbf{w}\rVert_1 with \lambda = 2\sigma^2/(Nb). It drives coefficients exactly to zero, so it selects features. Geometrically (Figure 1.10), penalising is equivalent to minimising the loss inside a ball \Omega(\mathbf{w}) \le t: grow the elliptical loss contours until they touch the ball. The L2 ball is round and is touched where every coordinate is non-zero. The L1 ball is a diamond with corners on the axes, and an ellipse coming from a general direction tends to meet a corner first.

The algebra is exact for an orthonormal design, \mathbf{X}^\top\mathbf{X} = N\mathbf{I}. The least-squares solution is then \hat{\mathbf{w}} = \mathbf{X}^\top\mathbf{y}/N, and expanding the square gives

which separates into one problem per coordinate.

• Ridge: minimise (w_j - \hat w_j)^2 + \lambda w_j^2; setting 2(w_j - \hat w_j) + 2\lambda w_j = 0 gives w_j = \hat w_j/(1 + \lambda).
• Lasso: minimise h(w_j) = (w_j - \hat w_j)^2 + \lambda\lvert w_j\rvert. For w_j > 0, 2(w_j - \hat w_j) + \lambda = 0 gives w_j = \hat w_j - \lambda/2, consistent only if \hat w_j > \lambda/2; symmetrically w_j = \hat w_j + \lambda/2 if \hat w_j < -\lambda/2. Otherwise the minimum is at the kink: at w_j = 0 the slope is -2\hat w_j + \lambda \ge 0 to the right and -2\hat w_j - \lambda \le 0 to the left. Together, soft-thresholding:

Ridge scales; the lasso subtracts and clips (Figure 1.11 shows both on real data). In general the lasso has no closed form, and coordinate descent, soft-thresholding one coordinate at a time, solves it. Among strongly correlated features the lasso keeps one more or less arbitrarily; the elastic net, which adds both penalties, keeps the group.

In practice

• Standardise first. The penalty treats all coefficients alike, but a coefficient’s size depends on its feature’s units: a length in millimetres needs a coefficient a thousand times smaller than in metres, and ridge would penalise it a million times less.
• Do not penalise the intercept. Centre \mathbf{y} and the features, or exclude b, as Lab 3 does by replacing \mathbf{I} with \mathbf{I}', whose first diagonal entry is zero.
• Choose \lambda on a logarithmic grid by validation or cross-validation. The one-standard-error rule (Hastie et al., 2009) takes the largest \lambda whose cross-validated error is within one standard error of the best.
• Map \lambda to the library, by dividing each library objective by a constant until it matches (9.1):

Other regularisers

Early stopping regularises through the optimiser. Gradient descent from \mathbf{w}_0 = \mathbf{0} shrinks the error along each eigenvector of \mathbf{H}, eigenvalue h_i here, by 1 - \eta h_i per step (Section 3), so after t steps that component of \mathbf{w}_t is the least-squares value times 1 - (1 - \eta h_i)^t. Ridge’s factor in the same basis is h_i/(h_i + 2\lambda), since h_i = 2s_i^2/N. Both are near 1 for large h_i; for small h_i they are about \eta t\,h_i and h_i/(2\lambda), which match when \lambda \approx 1/(2\eta t). Stopping early behaves roughly like ridge with a penalty that weakens as training continues. Data augmentation (Module 03, Section 10), and dropout and noise injection (Module 02, Section 11), regularise networks.

A test score is a fair estimate of the expected risk of Section 1 only if the test cases played no part in building the model, resemble the cases met in service, and the number is read with its sampling noise. This section is the series’ reference for standard errors, the bootstrap and McNemar’s test; Modules 07 and 09 apply them to language-model benchmarks.

How noisy a test number is

Score n independent cases with a model of true accuracy p. The number correct is binomial, so the observed accuracy \hat p has variance p(1 - p)/n and standard error

and a 95% interval is about \hat p \pm 1.96\,\operatorname{SE} by the central limit theorem. The formula holds for any proportion, with n the cases it is computed over: for recall, the positives. The normal approximation fails with fewer than about ten errors or ten successes; use the bootstrap or an exact binomial interval there. To size a test set for a half-width h, invert the formula:

Test-set reuse and selection

Test-set reuse is looking at the test number, changing the model and looking again; each change fits the model a little to the test set. The remedy is procedural: open the test set once and report the number with the date it was opened.

Selection does the same to a validation set in one step. Each configuration’s score is its true performance plus noise, and the largest of many scores belongs to the configuration whose noise was most favourable, so the winner is biased upwards even when no configuration is better.

The remedy is nested cross-validation (Figure 1.12): an inner loop inside each outer training portion does the selection, and the outer fold scores the model it chose, so the outer score measures the whole procedure. It costs k_\text{outer}\times k_\text{inner}\times the number of configurations fits, 5\times3\times20 = 300 for twenty values of \lambda. The deployed model comes from running the inner selection once on all the data.

Leakage

Leakage is any path by which information about the held-out labels reaches the model during training. Its forms, each with a fix:

1. Preprocessing fitted on all the data: scaling, imputation, feature selection, target encoding. Fit them inside each fold with a scikit-learn Pipeline. Scaling usually leaks little; selection and target encoding can leak a great deal.
2. Duplicates and near-duplicates across the split, such as a weld image saved twice. Deduplicate before splitting.
3. Groups: repeated measurements of one part, specimen, patient or machine, split by row. The model learns to recognise units. Split by group (GroupKFold).
4. Time: a shuffled time series lets the model train on the future. Use forward-chaining splits with a gap before each validation block (TimeSeriesSplit with gap).
5. Target leakage: a feature recorded after the label was known, such as the repair cost in a failure predictor. Check when each feature is available in service.
6. Test-set reuse, above.

The rule behind all six: split by the unit that will be new at deployment (a new part, a new patient, a new day), never by row (Figure 1.13).

Distribution shift

The test set is drawn from the same P as training; deployment is not. Under covariate shift the inputs change and the input–output relation does not: a new sensor supplier with a different offset. Under label shift the class proportions change: a process improvement halves the defect rate, and precision falls with the prevalence (Section 7). Under concept drift the relation itself changes: wear alters how vibration relates to remaining life. Evaluate on data from the conditions the model will meet, even if there is less of it, and in production monitor the inputs and outputs, not only the accuracy.

Bootstrap intervals

The bootstrap estimates sampling noise by resampling: draw n test cases with replacement, recompute the metric, repeat R = 2{,}000 times, and take the 2.5th and 97.5th percentiles as a 95% interval. It needs no formula, so it serves the AUC, F1 or ECE as well as accuracy.

To compare two models on one test set, use the paired bootstrap: draw the indices once per replicate, score both models on them and take the interval of the difference. Since \operatorname{Var}(\hat a - \hat b) = \operatorname{Var}\hat a + \operatorname{Var}\hat b - 2\operatorname{Cov}(\hat a, \hat b) and models find the same cases easy and the same ones hard, the paired interval is much narrower than two separate intervals suggest. With grouped data, resample groups (the cluster bootstrap).

Two caveats. The bootstrap over a test set measures test-sampling noise, not the variability from retraining. And it cannot invent evidence: when a difference rests on a handful of cases where the models disagree, every replicate is built from those few cases and the percentile interval is too narrow. Use the exact test below.

McNemar’s test

Score classifiers A and B on the same n cases. Cases both get right, or both get wrong, say nothing about which is better; only the discordant ones do: n_{10} that A alone gets right, n_{01} that B alone gets right. If the models are equally good, each discordant case is a fair coin, so n_{10} \sim \operatorname{Binomial}(n_{10} + n_{01}, \tfrac12), and the exact two-sided p-value is \min\big(1,\; 2\,P(X \le \min(n_{10}, n_{01}))\big) (scipy.stats.binomtest). This is McNemar’s test (McNemar, 1947), recommended for comparing classifiers by Dietterich (1998). The textbook statistic

compared with 3.84 (5% for \chi^2 with one degree of freedom), is a large-count approximation that rejects too easily at moderate counts; its continuity-corrected form (\lvert n_{10} - n_{01}\rvert - 1)^2/(n_{10} + n_{01}) tracks the exact test. The exact test costs nothing, so use it, and report n_{10} and n_{01} with it.

Applications built on language models raise evaluation problems of their own; the AI Agents series has a module on them.

Linear models are where every idea of this module is easiest to see, and often a good baseline. When the boundary bends or the features interact, other methods do better. This section is a tour with enough mechanism to know when each is the right baseline; the evaluation of Sections 7 to 10 is how you find out.

k-nearest neighbours and the curse of dimensionality

k-nearest neighbours (k-NN) predicts by averaging the targets, or taking a vote of the labels, of the k training points closest to the query. There is no training; a prediction costs about Nd operations to find the neighbours. k is the capacity knob: k = 1 interpolates the training data, and a large k averages so widely that it underfits. Distances mix features, so scaling is essential: unscaled, the feature with the largest units decides who is a neighbour. In a few dimensions k-NN is a strong baseline.

In many dimensions it fails, because neighbourhoods stop being local. For data spread uniformly in the unit cube, a sub-cube holding a fraction r of the data has edge r^{1/d}.

Trees, forests and boosting

A decision tree splits the input space recursively by axis-aligned thresholds, x_j \le t, choosing at each node the feature and threshold that most reduce the impurity of the two children: the variance of the targets for regression; for classification the Gini impurity 1 - \sum_k p_k^2 or the entropy, with p_k the class fractions in the node. Depth is capacity. An increasing rescaling of a feature moves a threshold on it without changing which points fall on each side, so trees need no scaling. They have high variance: a small change in the data can change the first split and everything below it.

A random forest grows M deep trees, each on a bootstrap sample of the rows and considering a random subset of the features at each split, and averages them. If each tree’s prediction has variance \sigma^2 and any two have correlation \rho, the average has variance

the M variances plus the M(M - 1) covariances. More trees remove only the second term; the first falls only if the trees are decorrelated, which is what the random feature subsets are for. Each bootstrap sample leaves out about a third of the rows, so scoring every row with the trees that did not see it gives the out-of-bag error, a validation estimate for free.

Gradient boosting builds an additive model F_m(\mathbf{x}) = F_{m-1}(\mathbf{x}) + \nu\,h_m(\mathbf{x}) from small trees h_m, each fitted to the negative gradient of the loss with respect to the prediction, evaluated at F_{m-1}. For squared error \tfrac12(y - F)^2 that is y - F, so each tree is fitted to the current residuals y_i - F_{m-1}(\mathbf{x}_i). It is gradient descent in function space: the learning rate \nu is traded against the number of trees M, which early stopping on validation data chooses. XGBoost, LightGBM and scikit-learn’s HistGradientBoosting estimators are the common implementations (as of 2026). Boosting needs no scaling and handles missing values natively. On tabular data of up to a few hundred thousand rows it is very often the best method and the first thing to try; Grinsztajn et al. (2022) give benchmark evidence that tree ensembles still beat deep networks on typical tabular data. “Very often” is not “always”: Exercise 13 tests the claim on two datasets, and on one of them ridge regression wins.

Kernels and support vector machines

Replace \mathbf{x} by a feature map \phi(\mathbf{x}), possibly of very high dimension, and stack the mapped inputs as the rows of \boldsymbol{\Phi}. Ridge regression gives \mathbf{w} = (\boldsymbol{\Phi}^\top\boldsymbol{\Phi} + \lambda N\mathbf{I})^{-1}\boldsymbol{\Phi}^\top\mathbf{y}. Multiplying out shows the push-through identity,

and multiplying on the left by (\boldsymbol{\Phi}^\top\boldsymbol{\Phi} + \lambda N\mathbf{I})^{-1} and on the right by (\boldsymbol{\Phi}\boldsymbol{\Phi}^\top + \lambda N\mathbf{I})^{-1} turns it into (\boldsymbol{\Phi}^\top\boldsymbol{\Phi} + \lambda N\mathbf{I})^{-1}\boldsymbol{\Phi}^\top = \boldsymbol{\Phi}^\top(\boldsymbol{\Phi}\boldsymbol{\Phi}^\top + \lambda N\mathbf{I})^{-1}. Hence

and a prediction is

Only inner products appear, so \phi never has to be formed: a kernel k computes them directly. This is the kernel trick, and the method is kernel ridge regression. The RBF kernel k(\mathbf{x}, \mathbf{x}') = \exp(-\gamma\lVert\mathbf{x} - \mathbf{x}'\rVert^2) corresponds to an infinite-dimensional \phi: in one dimension, expanding \exp(2\gamma xx') as a power series gives one feature e^{-\gamma x^2}\sqrt{(2\gamma)^j/j!}\,x^j for every j \ge 0. The price is the N\times N matrix \mathbf{K}: O(N^3) to fit, and N kernel evaluations per prediction.

The support vector machine (SVM), for labels y_i \in \{-1, +1\}, keeps the L2 penalty and replaces squared error by the hinge loss \max(0, 1 - y_i f(\mathbf{x}_i)). A point classified with margin y_i f(\mathbf{x}_i) \ge 1 contributes no loss and no gradient, so the solution depends only on the support vectors, the points on or inside the margin; all others have \alpha_i = 0, and prediction needs only the support vectors. The boundary is the one with the largest margin the penalty allows. With an RBF kernel, SVMs were the state of the art on medium-sized data before deep learning.

Gaussian processes and neural networks

A Gaussian process puts a prior distribution over functions, specified by a kernel, and returns a posterior mean and variance with every prediction. The posterior mean is kernel ridge regression with \lambda N equal to the noise variance. The cost is cubic in N, and the uncertainty is the point: Gaussian processes are the standard surrogate models for expensive engineering simulations and the engine of Bayesian optimisation.

Neural networks learn the features as well as the final linear map. Their advantage is on images, text and signals, where good features are not known and data are plentiful; they are the rest of this series.

Unsupervised workhorses: k-means and PCA

k-means partitions data into k clusters by minimising J = \sum_i\lVert\mathbf{x}_i - \boldsymbol{\mu}_{c(i)}\rVert^2. Lloyd’s algorithm alternates assigning each point to its nearest centre, which minimises J over the assignments, and moving each centre to its cluster’s mean, which minimises J over the centres (Section 5’s best constant under squared error). Neither step increases J, so it converges, but to a local minimum that depends on the start: use k-means++ initialisation and several restarts. Choose k by the downstream use, the elbow of J against k, or the silhouette score; scale first.

Principal component analysis (PCA) finds the directions of largest variance. The projection \mathbf{u}^\top\mathbf{x} with \lVert\mathbf{u}\rVert = 1 has variance \mathbf{u}^\top\mathbf{C}\mathbf{u}, for the sample covariance \mathbf{C}; maximising it with a Lagrange multiplier gives \mathbf{C}\mathbf{u} = \lambda\mathbf{u}, with variance \lambda (here \lambda is an eigenvalue again, as in Section 3, not the ridge strength of the kernel subsection above). The principal directions are the eigenvectors of \mathbf{C} in order of eigenvalue, equivalently the right singular vectors of the centred data. Component i explains the fraction \lambda_i/\sum_j\lambda_j of the variance, and the projection onto the first k is \mathbf{z} = \mathbf{V}_k^\top(\mathbf{x} - \boldsymbol{\mu}). Uses: visualisation, compression, denoising, and decorrelating features, which improves the conditioning of Section 3. PCA is linear and sensitive to scaling; its neural relative, the autoencoder, is in Module 05, Section 2.

No free lunch, and a comparison

Averaged over all possible problems, no learning method beats any other (Wolpert, 1996). A method wins on a problem because its assumptions match that problem, which is why the evaluation of Sections 7 to 10 matters more than the choice of method.