Statistical learning, from least squares to the exponential family

Why this subject matters now

It is tempting to treat classical statistical learning as a historical layer beneath deep learning, something to be skimmed on the way to transformers. The opposite is closer to the truth. The loss functions that train modern networks are exponential-family negative log-likelihoods. Softmax cross-entropy is the categorical GLM objective, mean squared error is the Gaussian one, and the reason both have the gradient "prediction minus target" is a single theorem about canonical links, derived below in four lines. The optimizer folklore of deep learning, why second-order methods converge in a handful of steps when the curvature fits in memory, why weight decay is a prior, why the Fisher information is the expected Hessian, is generalized-linear-model theory transplanted. Kernel methods, which looked superseded a decade ago, returned to the center of the field when Jacot, Gabriel, and Hongler at EPFL showed that infinitely wide networks train as kernel regression under the neural tangent kernel, and the sharpest available analyses of why overparameterized models generalize, benign overfitting from Bartlett and collaborators at Berkeley and the double-descent curve from Belkin and coauthors, are statements about ridge regression in high dimension.

Meanwhile the deployed world still runs on the classical stack directly. Insurance pricing, epidemiology, and demand forecasting run on Poisson and gamma GLMs because the coefficients are rate ratios that an actuary or a regulator can read. Advertising and recommendation systems run regularized logistic regression at billions of rows because it trains in one pass and calibrates. And on tabular data with heterogeneous columns, gradient-boosted trees still beat neural networks. Grinsztajn, Oyallon, and Varoquaux at INRIA benchmarked this carefully in 2022 across 45 datasets and found the gap persists even after substantial tuning of the neural models, a result that has held up.

What changed in the last five years is not that this material became obsolete but that the bar moved. Knowing the formulas is table stakes. What distinguishes a strong practitioner is being able to derive them, because the derivations are what transfer when the model in front of them is not in a textbook. This page is a single unified treatment rather than a tour. Shorter standalone pages exist here for linear regression, logistic regression, support vector machines, Gaussian mixtures, PCA, decision trees, gradient boosting, Naive Bayes, and k-means. This one derives the common framework those pages sit inside, and most of them collapse into special cases of two or three theorems.

The setup, from risk to the Bayes predictor and where error comes from

Risk and empirical risk

Fix a joint distribution \(\P\) over pairs \((x, y)\) with \(x \in \mathcal{X}\) and \(y \in \mathcal{Y}\). A predictor is a function \(h: \mathcal{X} \to \mathcal{Y}'\), a loss is a function \(\ell(\hat y, y) \ge 0\), and the object every learning algorithm is trying to make small is the risk, the expected loss on a fresh draw.

The distribution is unknown. What is available is a sample \(\D = \{(x_i, y_i)\}_{i=1}^{n}\) drawn i.i.d. from \(\P\). Replacing the expectation by the sample average gives the empirical risk \(\hat R_n(h) = \frac{1}{n}\sum_{i=1}^{n} \ell(h(x_i), y_i)\), and the induction principle almost all of supervised learning uses is empirical risk minimization, which chooses \(\hat h = \argmin_{h \in \mathcal{H}} \hat R_n(h)\) over some hypothesis class \(\mathcal{H}\). Everything interesting comes from the two gaps this introduces, the gap between \(\mathcal{H}\) and all functions, and the gap between \(\hat R_n\) and \(R\).

The Bayes-optimal predictor, derived

Before restricting to a hypothesis class, ask what the best possible predictor is when \(\P\) is known and \(h\) may be any function. Because the loss decomposes over points, the minimization can be done pointwise, by conditioning on \(x\) and minimizing the conditional expected loss separately at each \(x\).

Squared loss. With \(\ell(\hat y, y) = (\hat y - y)^2\), at a fixed \(x\) we minimize \(q(c) = \E[(c - y)^2 \mid x]\) over the scalar \(c\). Expanding,

which vanishes at \(c = \E[y \mid x]\), and \(q''(c) = 2 > 0\) confirms a minimum. Therefore

the conditional mean is the Bayes predictor for squared loss. Two consequences follow immediately. First, the minimum achievable risk is \(\E_x[\Var(y \mid x)]\), which is the irreducible noise, the part no amount of data or model capacity removes. Second, and this is the sentence worth remembering, every regression method that minimizes squared error is an attempt to estimate a conditional expectation, whether it is a linear model, a random forest, or a neural network. The same computation with absolute loss gives the conditional median (the derivative of \(\E[|c - y| \mid x]\) is \(\P(y < c \mid x) - \P(y > c \mid x)\), zero exactly at the median), and with the pinball loss at level \(\tau\) it gives the conditional \(\tau\)-quantile, which is why quantile regression exists.

Zero-one loss. For classification with \(y \in \{1, \dots, K\}\) and \(\ell(\hat y, y) = \mathbb{1}[\hat y \ne y]\), the conditional risk of predicting class \(k\) at \(x\) is

which is minimized by taking the largest posterior probability,

This is the Bayes classifier, and its risk \(\E_x[1 - \max_k \P(y=k \mid x)]\) is the Bayes error. Note what this argument does not require. There is no assumption about the shape of \(\P(y \mid x)\), no independence between features, nothing. It also explains why classification and probability estimation are different problems with different difficulty. Getting the argmax right needs only the correct ordering of the posteriors, which is why a badly calibrated model can still have excellent accuracy, and why the calibration section at the end of this page is not optional for anyone who uses the probabilities downstream.

Approximation, estimation, and optimization

Let \(h^\ast\) be the Bayes predictor, \(h^\ast_{\mathcal{H}} = \argmin_{h \in \mathcal{H}} R(h)\) the best predictor in the class, \(\hat h_n = \argmin_{h \in \mathcal{H}} \hat R_n(h)\) the empirical minimizer, and \(\tilde h_n\) whatever the optimizer actually returned after a finite compute budget. Add and subtract.

The identity is trivial. The content is in what controls each term. Approximation error is a property of \(\mathcal{H}\) alone and shrinks as the class grows richer. Estimation error is the price of using \(n\) samples instead of \(\P\), and it grows with the capacity of \(\mathcal{H}\). The uniform-convergence machinery in Shalev-Shwartz and Ben-David bounds it by quantities like the Rademacher complexity or the VC dimension, typically at rate \(O(\sqrt{\text{capacity}/n})\). Optimization error is the gap the solver leaves behind, and it is the only one of the three that a practitioner can watch directly by looking at the training loss. Bottou and Bousquet's observation, that at fixed compute the three terms trade off against each other and the optimal choice is usually to under-optimize a larger model rather than fully optimize a smaller one, is the reason large-scale learning uses SGD rather than Newton's method even though Newton converges in a handful of iterations. The measurement below quantifies both sides of that tradeoff.

A convention holds throughout. \(X \in \R^{n \times d}\) stacks inputs as rows, \(y \in \R^n\) stacks targets, a leading 1 is folded into each \(x_i\) so the intercept is a coordinate of \(\theta\), and \(\hat y\) denotes fitted values. Where a factor of \(\tfrac12\) appears in a loss it is there only to cancel the 2 from differentiating a square.

Least squares

The problem and the calculus derivation

The training set is \(n\) pairs \((x_i, y_i)\) with \(x_i \in \R^d\) and \(y_i \in \R\). The hypothesis is linear, \(h_\theta(x) = \theta\T x\), and the fitting criterion is the sum of squared residuals,

To differentiate, expand the quadratic.

Two matrix-calculus facts do the work, and both reduce to coordinates. For the linear term, \(\partial (b\T \theta)/\partial \theta_k = b_k\), so \(\nabla_\theta\, b\T\theta = b\), and here \(b = X\T y\). For the quadratic term with symmetric \(A = X\T X\), write \(\theta\T A \theta = \sum_{j,k} A_{jk}\theta_j \theta_k\). Differentiating with respect to \(\theta_m\) picks up \(A_{mk}\theta_k\) from the \(j{=}m\) terms and \(A_{jm}\theta_j\) from the \(k{=}m\) terms, giving \(2\sum_k A_{mk}\theta_k\) by symmetry, that is \(\nabla_\theta\, \theta\T A\theta = 2A\theta\). Therefore

and setting the gradient to zero gives the normal equations

The Hessian is \(X\T X\), which is positive semidefinite because \(v\T X\T X v = \lVert Xv\rVert^2 \ge 0\), so \(J\) is convex and the stationary point is a global minimum. When \(X\) has linearly dependent columns, \(X\T X\) is singular, the minimizer is a whole affine set, and the Moore-Penrose pseudoinverse \(\theta = X^{+} y\) selects the minimum-norm member. That choice is what lstsq routines return and it matters again in the discussion of interpolation and double descent later.

The projection view

The same equations fall out of geometry with no calculus at all. The set of achievable prediction vectors \(\{X\theta : \theta \in \R^d\}\) is the column space of \(X\), a \(d\)-dimensional subspace of \(\R^n\) when the columns are independent. Minimizing \(\lVert X\theta - y\rVert\) means finding the point of that subspace closest to \(y\), which is the orthogonal projection of \(y\) onto it. The defining property of the projection is that the residual \(r = y - X\theta\) is orthogonal to the subspace, hence orthogonal to every column of \(X\).

Same normal equations, now with content. Least squares does not make residuals small in every direction, it makes them exactly orthogonal to everything the model can express. The fitted values are \(\hat y = X(X\T X)^{-1}X\T y = Hy\), where \(H\) is the hat matrix. \(H\) is symmetric and idempotent, since \(H^2 = X(X\T X)^{-1}X\T X(X\T X)^{-1}X\T = X(X\T X)^{-1}X\T = H\), and symmetric idempotent matrices are exactly orthogonal projections. Its diagonal entries \(h_{ii} = x_i\T(X\T X)^{-1}x_i\), the leverages, measure how much observation \(i\) pulls its own fit. They satisfy \(0 \le h_{ii} \le 1\) and sum to \(\tr(H) = d\). The leverages return with a starring role in cross-validation.

The probabilistic view, where Gaussian noise makes least squares maximum likelihood

So far squared error was a choice. A generative assumption turns it into a consequence. Suppose the targets are generated as

so \(p(y_i \mid x_i; \theta) = \frac{1}{\sqrt{2\pi}\,\sigma} \exp\!\big({-\frac{(y_i - \theta\T x_i)^2}{2\sigma^2}}\big)\). Independence makes the likelihood of the dataset a product, and the log turns it into a sum.

The first term does not involve \(\theta\), and the second is \(-J(\theta)/\sigma^2\). Maximizing \(\ell\) over \(\theta\) is therefore exactly minimizing the sum of squared residuals, for any \(\sigma\). Under independent Gaussian noise of constant variance, the maximum-likelihood estimator is the least-squares estimator. The assumptions are doing real work in that sentence. Heavier-tailed Laplace noise would make the MLE minimize absolute error. Input-dependent variance would make it weighted least squares. Correlated noise with covariance \(\Sigma\) would make it generalized least squares with weights \(\Sigma^{-1}\). Maximizing over \(\sigma^2\) as well, set \(\partial \ell / \partial \sigma^2 = -n/(2\sigma^2) + \mathrm{RSS}/(2\sigma^4) = 0\) to get \(\hat\sigma^2 = \mathrm{RSS}/n\), the average squared residual, biased low by the factor \((n-d)/n\) relative to the unbiased estimator because the same data chose \(\theta\).

Gauss-Markov, and what it does not say

The Gauss-Markov theorem states that if \(\E[\varepsilon] = 0\), \(\Cov[\varepsilon] = \sigma^2 I\), and the model is correctly specified as linear in \(\theta\), then among all linear unbiased estimators \(\tilde\theta = C y\) with \(\E[\tilde\theta] = \theta\) for every \(\theta\), the least-squares estimator has the smallest covariance in the positive-semidefinite order. The proof is short enough to sketch honestly. Write \(C = (X\T X)^{-1}X\T + D\). Unbiasedness for all \(\theta\) forces \(CX = I\), hence \(DX = 0\). Then \(\Cov[\tilde\theta] = \sigma^2 C C\T = \sigma^2\big[(X\T X)^{-1} + D D\T\big]\), because the cross terms contain \(DX = 0\). Since \(DD\T \succeq 0\), the extra covariance is positive semidefinite and vanishes only when \(D = 0\).

Every word in the hypothesis is load-bearing, and the theorem's reputation exceeds its reach. It says nothing about biased estimators, and the whole point of ridge and lasso is that a small bias buys a large variance reduction. The exact computation in Problem 5 shows a biased estimator beating the Gauss-Markov optimum by 11% in mean squared error. It says nothing about nonlinear estimators. It assumes homoskedastic uncorrelated noise and correct specification, neither of which typically holds. And it is a statement about the estimator's variance, not about prediction error on new data. Gauss-Markov is a useful sanity anchor and a bad design principle.

The squared condition number, or why nobody forms \(X\T X\)

The formula \(\theta = (X\T X)^{-1}X\T y\) is a statement about mathematics, not an algorithm. Two things are wrong with implementing it literally. The smaller problem is cost, since an explicit inverse is roughly three times the work of a factorization and solve, and is never needed. The larger problem is conditioning.

The condition number of a matrix in the 2-norm is \(\kappa(X) = \sigma_{\max}(X)/\sigma_{\min}(X)\), the ratio of largest to smallest singular value. It measures how much a relative perturbation of the input can be amplified in the output. Since the singular values of \(X\T X\) are the squares of those of \(X\),

The moment the Gram matrix is formed, the conditioning of the problem is squared. A backward-stable solve of a system with condition number \(\kappa\) returns an answer with relative error on the order of \(\kappa \cdot \epsilon_{\text{mach}}\), so the normal-equations route delivers roughly \(\kappa(X)^2 \epsilon\) while a method that never forms the Gram matrix delivers \(\kappa(X)\epsilon\). In double precision (\(\epsilon \approx 2.22 \times 10^{-16}\)) that is the difference between a usable answer and noise as soon as \(\kappa(X)\) reaches about \(10^{8}\).

The alternative is to solve the least-squares problem through a factorization of \(X\) itself. For QR, write \(X = QR\) with \(Q \in \R^{n \times d}\) having orthonormal columns and \(R \in \R^{d \times d}\) upper triangular. Because \(Q\T Q = I\),

and the second term does not involve \(\theta\), so the minimizer solves the triangular system \(R\theta = Q\T y\) by back substitution in \(O(d^2)\). The whole computation costs about \(2nd^2 - \tfrac23 d^3\) flops and never squares the conditioning. Cholesky is the middle option. Form \(G = X\T X + \lambda I\), factor \(G = LL\T\) with \(L\) lower triangular, and solve two triangular systems. It costs about \(nd^2 + \tfrac13 d^3\), roughly half of QR, and it is what production ridge solvers actually use, because the explicit \(\lambda I\) floors the smallest eigenvalue at \(\lambda\) and caps \(\kappa(G)\) at \((\sigma_{\max}^2 + \lambda)/\lambda\). Cholesky on an unregularized ill-conditioned Gram matrix is exactly as bad as the normal equations, and worse in one respect. It fails outright, since a matrix that is positive definite in exact arithmetic can lose that property to rounding.

This is measurable, so it was measured. Designs \(X \in \R^{400 \times 30}\) were built with prescribed condition numbers by choosing the singular values directly, \(y = X\theta^\ast\) exactly (no noise, so the only error is numerical), and the same problem solved four ways. The table reports \(\lVert\hat\theta - \theta^\ast\rVert / \lVert\theta^\ast\rVert\).

Read the second column against the third. The normal-equations error tracks \(\kappa(X)^2 \epsilon\) almost exactly, rising by four orders of magnitude for every two orders in \(\kappa(X)\), and it reaches 100% relative error, which is to say total garbage, at \(\kappa(X) = 10^{8}\). At the same point QR is still accurate to nine digits. At \(\kappa(X) = 10^{10}\) the Cholesky factorization aborts because the Gram matrix is no longer numerically positive definite. In single precision (\(\epsilon \approx 1.19 \times 10^{-7}\)) the normal equations are already worthless at \(\kappa(X) = 10^{4}\), which is an unremarkable condition number for a design matrix with two nearly collinear columns. The practical rules are short. Use QR or SVD when accuracy matters and you cannot regularize. Use Cholesky on \(X\T X + \lambda I\) when you can, because regularization fixes the conditioning and Cholesky is twice as fast. Never call inv. On this machine, solving a \(500{,}000 \times 1024\) ridge problem in float64 on an H100 takes 0.055 s by Gram-plus-Cholesky and 0.174 s by QR, a factor of 3.2, with the two answers differing by \(8.6 \times 10^{-9}\).

Locally weighted regression and the parametric/non-parametric line

A straight line through all the data is a strong global commitment. Locally weighted linear regression keeps the linear machinery but makes the commitment local. To predict at a query point \(x\), solve a weighted least-squares problem that cares most about training points near \(x\),

then predict \(\theta\T x\) with the locally fitted \(\theta\). The derivation is the same calculus with a diagonal weight matrix \(W = \diag(w_1, \dots, w_n)\). The objective is \(\tfrac12 (X\theta - y)\T W (X\theta - y)\), the gradient is \(X\T W X \theta - X\T W y\), and the solution is \(\theta = (X\T W X)^{-1} X\T W y\). The bandwidth \(\tau\) controls the size of the neighborhood. Small \(\tau\) tracks local wiggles (low bias, high variance), while large \(\tau\) recovers plain least squares (high bias, low variance).

The distinction this draws is worth being precise about, because "non-parametric" is often used loosely. A parametric method summarizes the training data into a fixed-size parameter vector and then discards the data. Prediction cost is independent of \(n\), and the capacity of the model is fixed before seeing the data. A non-parametric method keeps the data (or a summary whose size grows with \(n\)) and lets the effective complexity grow with the sample. Locally weighted regression, \(k\)-nearest neighbours, kernel ridge regression, and Gaussian processes are all in this family. The price of locality is that nothing is learned once and reused. The whole training set must be kept, and every prediction costs a fresh \(O(nd^2 + d^3)\) solve. The benefit is that the model can represent any smooth function as \(n\) grows, with no approximation error in the limit. The same weighted solve, with weights supplied by the model itself rather than by distance to a query, is the inner loop of Newton's method for every GLM below, which is why it earns its derivation here.

Shrinkage, from ridge and lasso to the bias-variance trade in closed form

Ridge regression and its Bayesian reading

Ridge regression adds a quadratic penalty to least squares.

The regularized system is always solvable. \(X\T X\) is positive semidefinite with eigenvalues \(d_i^2 \ge 0\), so \(X\T X + \lambda I\) has eigenvalues \(d_i^2 + \lambda \ge \lambda > 0\) and is invertible regardless of collinearity or of \(d > n\). Substituting the SVD \(X = U D V\T\) makes the effect transparent.

Each singular direction is shrunk by the factor \(d_i^2/(d_i^2 + \lambda)\), which is near 1 for strong directions and near 0 for weak ones. That is precisely the right thing to do, because the variance of the unregularized coefficient along direction \(i\) is \(\sigma^2/d_i^2\), so the directions ridge crushes are exactly the ones where the data carries almost no signal and the estimate is almost all noise. In a measured example with \(n = 500\), \(d = 40\), and two nearly duplicated columns, the singular values run from \(32.16\) down to \(1.486 \times 10^{-5}\). At \(\lambda = 1\) the shrinkage factor is \(0.99903\) on the strongest direction and \(2.2 \times 10^{-10}\) on the weakest. Ridge deleted the degenerate direction and left everything else alone. The trace \(\sum_i d_i^2/(d_i^2+\lambda)\) is the effective degrees of freedom, a continuous version of "number of parameters" that decreases smoothly from \(d\) to 0 as \(\lambda\) grows, and it is what AIC-style criteria should use for a penalized fit.

The Bayesian reading makes \(\lambda\) interpretable. Put a prior \(\theta \sim \N(0, \tau^2 I)\) and keep the Gaussian noise model \(y_i \mid x_i, \theta \sim \N(\theta\T x_i, \sigma^2)\). The maximum a posteriori estimate maximizes the log posterior, which by Bayes' rule is the log likelihood plus the log prior plus a constant.

Multiplying through by \(-\sigma^2\), which preserves the argmax up to flipping to argmin, the MAP problem is exactly ridge with

Noisy data or a tight prior push \(\lambda\) up, clean data or a diffuse prior push it toward zero, and \(\tau \to \infty\) recovers plain maximum likelihood. Because the Gaussian likelihood and Gaussian prior are conjugate, the full posterior is available in closed form, \(\theta \mid \D \sim \N\big((X\T X + \lambda I)^{-1}X\T y, \sigma^2 (X\T X + \lambda I)^{-1}\big)\), and the ridge estimate is simultaneously the posterior mean and the posterior mode. One caution is worth stating. MAP with the prior treated as a tuning knob is not "being Bayesian", it is choosing a regularizer. The honest Bayesian object is the whole posterior, and the honest frequentist tool for choosing \(\lambda\) is cross-validation.

The lasso through subgradients and soft thresholding

Replacing the \(\ell_2\) penalty with \(\ell_1\) gives the lasso of Tibshirani (1996).

The objective is convex but not differentiable at any point where a coordinate is zero, so the stationarity condition uses subgradients. The subdifferential of \(|t|\) is \(\{\mathrm{sign}(t)\}\) for \(t \ne 0\) and the interval \([-1, 1]\) at \(t = 0\). A convex function is minimized exactly where \(0\) belongs to its subdifferential, so \(\hat\theta\) is optimal if and only if, coordinatewise,

These are the KKT conditions of the lasso, and they already explain sparsity. A coordinate can sit exactly at zero over a whole interval of correlations, because zero is optimal as long as the corresponding gradient component stays inside \([-\lambda, \lambda]\). The \(\ell_2\) penalty has no such interval, since its gradient \(\lambda\theta_j\) vanishes at zero and cannot balance a nonzero data gradient. Ridge shrinks everything and zeroes nothing.

The orthogonal case, solved exactly. Suppose the design is orthogonal, \(X\T X = I\) (rescale if it is \(cI\)). Then \(\lVert X\theta - y \rVert^2 = \lVert\theta\rVert^2 - 2\theta\T X\T y + \lVert y\rVert^2\), and writing \(z = X\T y\) for the OLS solution the objective separates across coordinates.

Solve one coordinate. If \(\theta_j > 0\), the derivative is \(\theta_j - z_j + \lambda = 0\), giving \(\theta_j = z_j - \lambda\), which is consistent with \(\theta_j > 0\) only when \(z_j > \lambda\). If \(\theta_j < 0\), the derivative is \(\theta_j - z_j - \lambda = 0\), giving \(\theta_j = z_j + \lambda\), consistent only when \(z_j < -\lambda\). Otherwise the optimality interval at zero applies, \(-z_j \in [-\lambda, \lambda]\), that is \(|z_j| \le \lambda\), and \(\theta_j = 0\). Collecting the three cases gives the soft-thresholding operator

where \((t)_+ = \max(t, 0)\). Compare with ridge on the same orthogonal design, whose solution is \(\hat\theta_j = z_j/(1+\lambda)\), a proportional shrink, and with best-subset selection, whose solution is hard thresholding \(z_j \mathbb{1}[|z_j| > \sqrt{2\lambda}]\). The lasso is the convex relaxation sitting between them. It sets small coefficients exactly to zero like subset selection, but shrinks the survivors by a constant \(\lambda\) rather than leaving them alone, which is the source of its well-known bias on large coefficients and the motivation for the relaxed lasso and for the elastic net of Zou and Hastie (2005), which adds a small ridge term to handle correlated groups of predictors that plain lasso picks among arbitrarily.

This was verified numerically. On an orthogonal design with \(X\T X = nI\), \(n = 200\), \(d = 8\), the soft-thresholding formula and scikit-learn's coordinate-descent lasso agree to \(8.9 \times 10^{-16}\) at every penalty level tested, and the number of nonzero coefficients falls 5, 4, 3, 3 as \(\alpha\) goes 0.05, 0.2, 0.5, 1.0. The general (non-orthogonal) lasso has no closed form, and the standard solvers are coordinate descent, which applies exactly this soft-threshold to one coordinate at a time holding the others fixed, and proximal gradient methods, whose accelerated form (FISTA, Beck and Teboulle at the Technion and Tel Aviv, 2009) reaches \(O(1/k^2)\) convergence using the same soft-threshold as its proximal operator. Both work because \(\mathcal{S}_\lambda\) is the proximal operator of the \(\ell_1\) norm.

Logistic regression

Where the sigmoid comes from

For binary classification, \(y \in \{0, 1\}\), the natural model of the conditional is a Bernoulli with some mean \(\phi = p(y{=}1)\). Write the Bernoulli mass function in exponential form.

The coefficient multiplying the sufficient statistic \(y\) is the natural parameter \(\eta = \log\frac{\phi}{1-\phi}\), the log-odds. Inverting it, \(e^\eta = \phi/(1-\phi)\), so \(\phi\,(1 + e^\eta) = e^\eta\) and

The sigmoid is not an aesthetic choice among S-shaped curves. It is the unique function that maps the natural parameter of a Bernoulli back to its mean. Logistic regression is the decision to make the natural parameter linear in the features, \(\eta = \theta\T x\), giving \(h_\theta(x) = \sigma(\theta\T x) = p(y{=}1 \mid x; \theta)\). Everything else, the gradient, the Hessian, the convexity, follows mechanically, and the general version of this construction for any exponential-family output is the GLM section below.

Gradient and Hessian of the log-likelihood

Both labels can be written in one expression, \(p(y \mid x; \theta) = h^y (1-h)^{1-y}\) with \(h = \sigma(\theta\T x)\), so the negative log-likelihood over the dataset is the cross-entropy

The one identity needed is the sigmoid's derivative. From \(\sigma(z) = (1+e^{-z})^{-1}\),

Differentiate one term of \(\L\) with respect to \(\theta\) by the chain rule, with \(z_i = \theta\T x_i\), so that \(\nabla_\theta h_i = h_i(1-h_i)x_i\).

Expanding the bracket, \(\frac{y_i}{h_i}\, h_i(1-h_i) = y_i(1-h_i)\) and \(\frac{1-y_i}{1-h_i}\, h_i(1-h_i) = (1-y_i)h_i\), so the whole expression is \(-\big(y_i - y_i h_i - h_i + y_i h_i\big)x_i = (h_i - y_i)x_i\). The \(y_i h_i\) terms cancel, which is the small miracle that makes the gradient so clean. Summing and stacking into matrices,

the "prediction minus label, times input" form that will turn out to be universal across GLMs. Differentiating once more, only \(h_i\) depends on \(\theta\), so

Convexity. For any \(v \in \R^d\),

valid because every \(h_i \in (0,1)\) makes \(h_i(1-h_i) > 0\) so the square root is real. The Hessian is positive semidefinite everywhere, the negative log-likelihood is convex, and any stationary point is a global minimum. It is positive definite, hence the minimum unique, exactly when \(X\) has full column rank. One honest caveat applies. If the classes are linearly separable, no stationary point exists, because scaling any separating \(\theta\) up forever increases the likelihood and the MLE runs off to infinity. Adding \(\tfrac{\lambda}{2}\lVert\theta\rVert^2\) restores a finite unique optimum, since the penalty eventually dominates, which is why production solvers regularize by default and why scikit-learn's LogisticRegression ships with \(C = 1\) rather than \(C = \infty\).

Newton's method and IRLS

Newton's method minimizes a function by repeatedly jumping to the minimum of its local quadratic model. Taylor-expanding \(\L\) around the current iterate \(\theta_t\),

with \(g = \nabla \L(\theta_t)\) and \(H = \nabla^2 \L(\theta_t)\). Setting the gradient of the right-hand side to zero gives \(g + H(\theta - \theta_t) = 0\), that is the update

Near the optimum the quadratic model is nearly exact and the error squares at every step. The worked problem below shows the digit-doubling live. The update has a second reading that explains the name iteratively reweighted least squares. Substitute \(g = X\T(h-y)\) and factor \(X\T W X\) out of both terms.

Each Newton step is exactly a weighted least-squares solve, the locally weighted machinery from the previous section, against a working response \(z\), the current linear predictions corrected by the residual scaled up by the inverse variance. Points with \(h_i\) near 0 or 1 have tiny weight \(h_i(1-h_i)\) and correspondingly amplified working residuals. Points near the decision boundary carry the most weight, which is the statistical statement that a confidently classified point contains little information about where the boundary is. Every classical GLM fitter, from the original Nelder-Wedderburn formulation of 1972 through statsmodels today, is this loop.

The convergence is fast and it is worth having a number for it. On a synthetic problem with \(n = 4000\) and \(d = 12\) plus an \(\ell_2\) penalty, Newton converged in 7 iterations with gradient norms \(1292 \to 308.6 \to 74.5 \to 8.22 \to 0.131 \to 3.49\times10^{-5} \to 0\), the number of correct digits roughly doubling per step once the iterate is close, and the final gradient norm \(1.06 \times 10^{-13}\). The same problem solved by scikit-learn's L-BFGS at its tightest tolerance agreed to \(5.0 \times 10^{-6}\) in the coefficients, which is L-BFGS's stopping accuracy, not ours. The contrast with first-order methods on an ill-conditioned design is stark. On a problem with \(n = 200{,}000\), \(d = 128\), and \(\kappa(X\T X) = 6.7 \times 10^{8}\), Newton reached the optimum in 6 iterations and 0.036 s on an H100, while full-batch gradient descent at the optimal fixed step size \(1/L\) needed 180,350 iterations and 71.8 s to reach an objective within \(10^{-3}\) of the same optimum, a factor of 2000 in wall-clock time. That is the whole argument for second-order methods when \(d\) is small enough that a \(d \times d\) solve is affordable, and the whole reason they are abandoned when it is not.

Multinomial logistic regression and the identifiability subtlety

With \(K\) classes, model \(p(y = k \mid x) = \softmax_k(\Theta x)\) where \(\Theta \in \R^{K \times d}\) and

The negative log-likelihood over the data, with one-hot labels \(Y \in \{0,1\}^{n\times K}\), is

Differentiate with respect to \(z_{ik}\). The first term contributes \(-Y_{ik}\), and the log-sum-exp term contributes \(\sum_{k'} Y_{ik'} \cdot \frac{e^{z_{ik}}}{\sum_m e^{z_{im}}} = p_{ik}\) because the one-hot row sums to 1. Hence \(\partial \L / \partial z_{ik} = p_{ik} - Y_{ik}\), and by the chain rule through \(z_i = \Theta x_i\),

again "probabilities minus one-hot, times input". This is the gradient that every deep learning framework computes for a classification head, and it is why the softmax and the cross-entropy are almost always fused into one op. Computed separately, the softmax produces probabilities that are then logged, losing precision, whereas the fused version evaluates \(\log\sum_m e^{z_m}\) with the max subtracted and returns \(p - Y\) directly.

The identifiability subtlety. The softmax is invariant to adding any constant to all its inputs, \(\softmax(z + c\mathbf{1}) = \softmax(z)\). Therefore replacing \(\Theta\) by \(\Theta + \mathbf{1}c\T\) for any \(c \in \R^d\), that is adding the same row vector to every class's weight row, changes no predicted probability at all. The parameters are identified only up to this \(d\)-dimensional shift, the Hessian is singular with a \(d\)-dimensional null space, and the MLE is a flat manifold, not a point. Measured directly, shifting a fitted \(\Theta\) by a random row vector of norm 2.69 changed every predicted probability by at most \(8.9 \times 10^{-16}\) and left the negative log-likelihood identical to eight decimals at \(0.55818989\).

Three standard fixes, each with a different flavour. Pin one class's row to zero, which recovers binary logistic regression when \(K = 2\) and is what econometrics software does, because then the coefficients read as log-odds against a reference category. Add an \(\ell_2\) penalty, which is strictly convex and so selects the unique minimum-norm representative, the one whose weight rows sum to zero. This is what scikit-learn and every neural network with weight decay do implicitly. Or leave it alone and accept that the individual coefficients are meaningless while the predictions are not, which is what deep learning does. The practical warning is to never interpret the raw magnitude of a softmax weight in an unregularized multi-class fit, and to be suspicious of any feature-importance method that does. Incidentally, gradient descent from a zero initialization stays in the zero-row-sum subspace forever, because the gradient \((P-Y)\T X\) has columns summing to zero (rows of \(P\) and of \(Y\) both sum to 1). The measurement confirms the fitted column sums are \(0\) to machine precision even without a penalty. Optimizers hide the problem rather than solving it.

Generalized linear models, the unification

This section is the centerpiece of the page. Least squares and logistic regression have been derived separately, and they looked different, one solving a linear system, the other iterating. They are the same model. The object that makes them the same is the exponential family, and the theorem that does the work is one line of calculus applied to a normalization constant.

The exponential family

A family of distributions is an exponential family if its density or mass function can be written as

with natural parameter \(\eta\), sufficient statistic \(T(y)\) (usually just \(y\)), log partition function \(a(\eta)\), and base measure \(b(y)\). The role of \(a\) is to normalize,

with the integral replaced by a sum for discrete \(y\). Fixing \(T\), \(a\), and \(b\) and varying \(\eta\) traces out a family of distributions, and different choices of \(T, a, b\) give different families.

The Bernoulli was put in this form above, with \(T(y) = y\), \(\eta = \log\frac{\phi}{1-\phi}\), \(b(y) = 1\), and \(a(\eta) = -\log(1-\phi) = \log(1 + e^\eta)\), the softplus. The Gaussian with unit variance also fits. Expanding the square,

so \(\eta = \mu\), \(T(y) = y\), \(a(\eta) = \eta^2/2\). For the Poisson, \(p(y; \lambda) = \frac{\lambda^y e^{-\lambda}}{y!} = \frac{1}{y!}\exp(y \log\lambda - \lambda)\), so \(\eta = \log \lambda\), \(T(y) = y\), \(b(y) = 1/y!\), and \(a(\eta) = \lambda = e^\eta\). The multinomial with \(K\) categories is worked in Problem 4.

The log partition function generates the moments

The function \(a\) is not bookkeeping, it is a moment-generating machine. Start from the fact that the density integrates to one, for every \(\eta\).

Differentiate both sides with respect to \(\eta\) (scalar case, and differentiation under the integral is legitimate here because exponential families are smooth in the interior of their natural parameter space). The right side gives 0. On the left, the chain rule brings down \(T(y) - a'(\eta)\),

because the integrand is the density times \((T(y) - a'(\eta))\) and \(a'(\eta)\) is a constant with respect to \(y\). Therefore

Differentiate the same identity a second time. The derivative of the integral of \(p \cdot (T - a')\) is \(\int p \,(T - a')^2 \,dy - \int p\, a''\, dy = \E[(T - \E T)^2] - a''(\eta)\), so

The mean is the first derivative of the log partition function, the variance the second. Check it against each member. For the Bernoulli, \(a(\eta) = \log(1+e^\eta)\), so \(a'(\eta) = \frac{e^\eta}{1+e^\eta} = \sigma(\eta) = \phi\) and \(a''(\eta) = \sigma(1-\sigma) = \phi(1-\phi)\), both correct. For the Gaussian, \(a(\eta) = \eta^2/2\), so \(a' = \eta = \mu\) and \(a'' = 1 = \sigma^2\), correct. For the Poisson, \(a(\eta) = e^\eta\), so \(a' = a'' = e^\eta = \lambda\), which recovers the Poisson's defining property that its mean equals its variance.

Two structural consequences follow immediately and are used repeatedly below. First, since \(a''\) is a variance it is nonnegative, so \(a\) is convex. Second, the map \(\eta \mapsto \mu = a'(\eta)\) is strictly increasing wherever the variance is strictly positive, hence invertible. Its inverse \(\eta = (a')^{-1}(\mu)\) is the canonical link function. For the Bernoulli that inverse is the logit, for the Poisson it is the log, and for the Gaussian it is the identity. The link is not chosen, it is computed from the distribution.

The GLM construction and the one gradient they share

A generalized linear model makes three assumptions.

1. Given \(x\), the target follows an exponential-family distribution with natural parameter \(\eta\).
2. The natural parameter is linear in the features, \(\eta = \theta\T x\). This is the canonical link choice. A non-canonical link inserts a function \(g\) so that \(g(\mu) = \theta\T x\) with \(g \ne (a')^{-1}\).
3. The prediction is the conditional mean, \(h_\theta(x) = \E[y \mid x] = a'(\theta\T x)\).

The negative log-likelihood of one observation, dropping the \(\theta\)-free \(\log b(y)\), is

and the chain rule, using \(\nabla_\theta (\theta\T x_i) = x_i\), gives its gradient in one line.

For every canonical-link GLM, the gradient of the negative log-likelihood is "prediction minus target, times input". Least squares (\(a' = \mathrm{id}\)), logistic regression (\(a' = \sigma\)), Poisson regression (\(a' = \exp\)), and softmax regression (\(a' = \softmax\)) all share it, and nothing changes between them except the mean function. The stationarity condition \(X\T(\hat y - y) = 0\) is the same orthogonality statement that defined least squares by projection, now holding for every member of the family. The residual is orthogonal to the column space of the design.

Differentiate once more.

positive semidefinite because \(a''\) is a variance. So every canonical GLM is a convex optimization problem, and Newton's method on it is IRLS with weights equal to the conditional variances, exactly as derived for logistic regression but now with the variance function supplied by \(a''\). This also explains a fact about deep learning that is often stated as folklore. The reason softmax cross-entropy plus a final linear layer produces the tidy "probabilities minus one-hot" gradient in every framework is that the last layer of a classifier is a canonical GLM, whatever nonlinear feature extractor sits below it. The network learns \(x \mapsto \phi(x)\), and the head does GLM regression on \(\phi(x)\).

One implementation follows one derivation. In the code section a single 12-line glm_irls routine takes a family name, looks up the pair \((a', a'')\), and fits all three of Gaussian, Bernoulli, and Poisson. Verified against independent references, the Gaussian fit matches the OLS solution to \(2.2 \times 10^{-11}\) (it converges in a single step, because the quadratic model of a quadratic is exact), the Bernoulli fit matches scikit-learn's LogisticRegression to \(1.4 \times 10^{-6}\), and the Poisson fit matches scikit-learn's PoissonRegressor to \(1.0 \times 10^{-8}\), each in 7 IRLS iterations.

Poisson regression, worked

Counts (arrivals, insurance claims, page views, disease incidence) get \(y \mid x \sim \text{Poisson}(e^{\theta\T x})\). The exponential mean function guarantees positivity and makes effects multiplicative. A unit increase in feature \(j\) multiplies the predicted rate by \(e^{\theta_j}\), which is why insurance and epidemiology report rate ratios rather than differences. The variance function \(a'' = \mu\) means the model asserts that variance equals mean. When real counts are more dispersed than that, which is the common case, the standard remedies are a quasi-Poisson dispersion parameter or a negative binomial model, which adds one shape parameter and nests Poisson as a limit.

Generative against discriminative

Logistic regression models \(p(y \mid x)\) directly and never commits to how \(x\) is distributed. A generative classifier models the class-conditional densities \(p(x \mid y)\) and the prior \(p(y)\), then classifies through Bayes' rule, \(p(y \mid x) \propto p(x \mid y)\, p(y)\). The two approaches answer the same question from opposite directions, and the comparison is one of the cleanest theory results in the subject.

Gaussian discriminant analysis

The model takes \(y \sim \text{Bernoulli}(\phi)\) and \(x \mid y{=}k \sim \N(\mu_k, \Sigma)\) with a shared covariance. The joint log-likelihood over the data is

Maximizing over \(\phi\), the first sum is \(\sum_i [y_i \log\phi + (1-y_i)\log(1-\phi)]\), whose derivative \(\frac{n_1}{\phi} - \frac{n - n_1}{1-\phi}\) vanishes at \(\hat\phi = n_1 / n\), the class-1 fraction, with \(n_1 = \sum_i y_i\). Maximizing over \(\mu_1\), only class-1 terms matter, each contributing \(-\tfrac12 (x_i - \mu_1)\T \Sigma^{-1} (x_i - \mu_1)\), so the gradient with respect to \(\mu_1\) is \(\Sigma^{-1} \sum_{i: y_i = 1} (x_i - \mu_1)\), which vanishes at the class mean \(\hat\mu_1 = \frac{1}{n_1}\sum_{i:y_i=1} x_i\), and symmetrically for \(\hat\mu_0\). For \(\Sigma\), parameterize by the precision \(\Lambda = \Sigma^{-1}\). The Gaussian terms contribute \(\tfrac{n}{2}\log\lvert\Lambda\rvert - \tfrac12 \sum_i (x_i - \mu_{y_i})\T \Lambda\, (x_i - \mu_{y_i})\). Two matrix derivatives, both derivable entrywise, do the rest, \(\partial \log\lvert\Lambda\rvert / \partial \Lambda = \Lambda^{-1}\) from the cofactor expansion of the determinant, and \(\partial (v\T \Lambda v) / \partial \Lambda = v v\T\). Setting the gradient to zero,

the pooled within-class covariance. Every estimate is a counting or averaging operation, no iteration, no solver, one pass over the data. That is the practical appeal of generative training, and it is why linear discriminant analysis is still the fastest reasonable classifier to fit on wide data.

The posterior of shared-covariance GDA is exactly a sigmoid

Compute what GDA believes about \(y\) given \(x\).

where dividing numerator and denominator by the numerator identifies \(A(x) = \log\frac{\phi\,\N(x;\mu_1,\Sigma)}{(1-\phi)\,\N(x;\mu_0,\Sigma)}\). The normalizing constants \((2\pi)^{-d/2}|\Sigma|^{-1/2}\) are identical in both Gaussians and cancel, leaving

Expand both quadratic forms. The pure-\(x\) pieces are \(-\tfrac12 x\T\Sigma^{-1}x\) in each, with opposite signs, and they cancel because the covariance is shared. The cross terms give \(+\mu_1\T\Sigma^{-1}x - \mu_0\T\Sigma^{-1}x\), and the constant terms give \(-\tfrac12\mu_1\T\Sigma^{-1}\mu_1 + \tfrac12\mu_0\T\Sigma^{-1}\mu_0\). Therefore

which is affine in \(x\), so the GDA posterior is exactly \(\sigma(\theta\T x + \theta_0)\) with \(\theta = \Sigma^{-1}(\mu_1 - \mu_0)\). If the covariances are allowed to differ between classes the \(x\T\Sigma^{-1}x\) terms no longer cancel and the boundary becomes quadratic. That model is QDA, and it costs \(K d(d+1)/2\) covariance parameters instead of \(d(d+1)/2\).

For a one-dimensional numerical check of the algebra, take \(\mu_0 = -1\), \(\mu_1 = 2\), \(\sigma^2 = 1.5\), \(\phi = 0.4\). The formulas give \(\theta_1 = (2-(-1))/1.5 = 2\) and \(\theta_0 = \frac{1 - 4}{2(1.5)} + \log\frac{0.4}{0.6} = -1 - 0.4054651 = -1.4054651\). Evaluating both the direct Bayes ratio and \(\sigma(2x - 1.4054651)\) at \(x = 0.7\) yields \(0.4986337264\) from each, agreeing to all ten printed digits.

The converse fails, and that asymmetry is the point. A logistic-form posterior does not imply Gaussian class-conditionals. Poisson class-conditionals with a shared rate structure produce one too, as does any pair of exponential-family densities from the same family with a shared dispersion. GDA makes strictly stronger assumptions than logistic regression, so the set of problems on which GDA is correct is a strict subset of those on which logistic regression is correct.

The efficiency/robustness tradeoff, measured

The classical results quantify the exchange. Efron (1975) showed that when the Gaussian assumption holds, discriminant analysis is asymptotically more efficient than logistic regression, which pays a substantial efficiency penalty at moderate class separation. Ng and Jordan, then at Berkeley, sharpened the other side in 2001. The generative model converges to its asymptotic error at rate \(O(\log d)\) in sample size against \(O(d)\) for the discriminative one, so the generative model can win small-sample races even when its assumptions are wrong enough to lose in the limit. The picture is two learning curves that cross.

The crossing was measured directly. With \(d = 10\), 200 replicates per point, and a 20,000-point test set, LDA (shared-covariance GDA) and unregularized logistic regression were fit on samples of increasing size, once on genuinely Gaussian class-conditionals with a shared covariance, and once on heavy-tailed skewed class-conditionals built from \(t_2\) plus exponential noise, where the Gaussian assumption is badly wrong.

Both halves of the theory are visible. On the left, where the generative model is correct, GDA is ahead at every sample size and the advantage is largest when data is scarce (0.2315 against 0.2550 at \(n=40\), a 9% relative reduction) and shrinks to nothing by \(n = 2560\), where both estimators have essentially found the Bayes rule. On the right, where the generative model is wrong, GDA is still ahead in the middle of the range, at \(n = 40\) through \(160\), because its variance advantage outweighs its bias, and then logistic regression overtakes it and stays ahead, 0.1980 against 0.2023 at \(n = 2560\), a gap that will not close, because GDA's asymptotic error is simply higher. The practitioner's reading is that on small, wide problems a strong generative assumption is often worth its bias, and on large problems it is not. This is the same tradeoff that makes Naive Bayes a good baseline on a thousand documents and a bad model on ten million.

Naive Bayes and Laplace smoothing

For high-dimensional discrete \(x\), documents represented as word counts being the canonical case, modeling \(p(x \mid y)\) in full is hopeless, since a vocabulary of 50,000 binary features would need a table with \(2^{50000}\) entries. Naive Bayes assumes the features are conditionally independent given the class, \(p(x \mid y) = \prod_j p(x_j \mid y)\), collapsing the table to per-word probabilities. The assumption is almost never true, and the classifier works anyway, for a reason worth being precise about. The estimated posteriors are badly wrong (they are usually pushed toward 0 or 1 because correlated evidence is counted repeatedly), but the ordering of the posteriors is often right, and the Bayes rule only needs the argmax. Naive Bayes is a good classifier and a terrible probability estimator.

The MLE for each word probability is a count ratio, and therein lies a trap. A word never seen in a class gets probability zero, and one such word in a test document drives the entire class's likelihood to zero regardless of all other evidence. Laplace smoothing adds one phantom occurrence of every word to every class. This is not a hack. Placing a uniform Dirichlet prior \(\text{Dir}(1, \dots, 1)\) on the word distribution and taking the posterior mean gives

exactly, where \(N_{jk}\) counts word \(j\) in class \(k\), \(N_k\) is the class's total token count, and \(V\) the vocabulary size. A general \(\text{Dir}(\alpha, \dots, \alpha)\) prior gives \((N_{jk}+\alpha)/(N_k + \alpha V)\), and \(\alpha\) becomes a tuning parameter that trades trust in the counts against trust in the uniform prior.

The size of the effect is easy to underestimate. On a synthetic two-topic corpus with a 5000-word vocabulary, 25 tokens per document, and only 300 training documents, roughly 3900 of the 5000 words are unseen in each class and 99.95% of test documents contain at least one word unseen in one of the classes. Sweeping \(\alpha\) produced the following.

Effectively unsmoothed Naive Bayes gets 67.6% accuracy and assigns the true class a mean log-probability of \(-10.68\), which is to say a probability of about \(2 \times 10^{-5}\). It is confidently and badly wrong on most documents. Laplace smoothing at \(\alpha = 1\) lifts accuracy to 83.6%, and the best value here, \(\alpha = 5\), reaches 85.5% and beats a regularized logistic regression on the same data (82.6%). Too much smoothing hurts again at \(\alpha = 20\), as the prior begins to drown the counts. A single pseudocount is the difference between a usable model and an unusable one.

Kernels

The dual of ridge regression

Ridge regression has a second form that looks nothing like the first and is often far cheaper. Start from the stationarity condition \(X\T X\theta - X\T y + \lambda \theta = 0\) and solve for \(\theta\) without inverting anything.

The solution is a linear combination of the training inputs, with one coefficient per training point. That is already the representer theorem in the special case at hand, and it fell out of the stationarity condition without any extra assumption. Substituting \(\theta = X\T\alpha\) back into the definition of \(\alpha\),

and predictions at a new point are \(\hat y(x) = \theta\T x = \alpha\T X x = \sum_{i=1}^{n}\alpha_i \langle x_i, x\rangle\).

Two things changed. The system to solve is \(n \times n\) instead of \(d \times d\), which is a win exactly when \(d > n\). And, decisively, the data enters only through inner products, both in \(K_{ij} = \langle x_i, x_j\rangle\) during training and in \(\langle x_i, x\rangle\) at prediction time. Nothing in the derivation used the fact that \(\langle \cdot,\cdot\rangle\) was the Euclidean inner product on \(\R^d\). Replace it with the inner product in any feature space, \(k(x, x') = \langle \phi(x), \phi(x')\rangle\), and the same two equations fit a linear model in \(\phi\)-space without ever computing \(\phi\). That substitution is the kernel trick. As verification, on a 200-point problem the primal solution \((X\T X + \lambda I)^{-1}X\T y\) and the dual solution \(X\T(K+\lambda I)^{-1}y\) agree to \(3.3 \times 10^{-15}\), and an RBF kernel ridge regression implemented from the dual matches scikit-learn's KernelRidge to \(1.6 \times 10^{-15}\) in test predictions.

The representer theorem

The general statement covers any loss and any strictly increasing regularizer. Let \(\mathcal{H}\) be a reproducing kernel Hilbert space with kernel \(k\), and consider

with \(\Omega\) strictly increasing. The claim is that any minimizer has the form \(f(\cdot) = \sum_{i=1}^{n}\alpha_i\, k(x_i, \cdot)\).

Proof. Let \(V = \mathrm{span}\{k(x_1,\cdot), \dots, k(x_n,\cdot)\} \subseteq \mathcal{H}\) and decompose any candidate \(f = f_\parallel + f_\perp\) with \(f_\parallel \in V\) and \(f_\perp \perp V\). The reproducing property says \(f(x_i) = \langle f, k(x_i, \cdot)\rangle_{\mathcal{H}}\), and since \(k(x_i, \cdot) \in V\) while \(f_\perp \perp V\),

So the loss term depends on \(f_\perp\) not at all. Meanwhile, by orthogonality, \(\lVert f \rVert^2 = \lVert f_\parallel\rVert^2 + \lVert f_\perp \rVert^2 \ge \lVert f_\parallel\rVert^2\), with equality only if \(f_\perp = 0\), and \(\Omega\) is strictly increasing, so the regularizer is strictly larger unless \(f_\perp = 0\). Dropping \(f_\perp\) therefore leaves the loss unchanged and strictly decreases the objective unless it was already zero. Any minimizer has \(f_\perp = 0\), which is the claim. \(\square\)

The content is a dimension reduction that is exact rather than approximate. An optimization over a space that may be infinite-dimensional collapses to an optimization over \(n\) numbers. It also draws the boundary of the method. Strict monotonicity of \(\Omega\) is essential. Drop it and the theorem fails. And the resulting \(n \times n\) system is why kernel methods stopped scaling, since \(K\) needs \(O(n^2)\) memory and \(O(n^3)\) time to factor, so at \(n = 10^6\) the Gram matrix alone would be 8 TB in float64.

Mercer's condition and what counts as a kernel

Which functions \(k(x, x')\) are inner products in some feature space? The answer is exactly the symmetric positive semidefinite ones. If \(k(x,x') = \phi(x)\T\phi(x')\), then for any points \(x_1, \dots, x_m\) and any \(c \in \R^m\),

so the Gram matrix is PSD, and necessity is that easy. Mercer's theorem supplies the converse. If \(k\) is continuous, symmetric, and produces a PSD Gram matrix for every finite point set, then there is a feature map \(\phi\) (into \(\ell^2\), possibly infinite-dimensional) with \(k(x,x') = \langle\phi(x),\phi(x')\rangle\). The construction goes through the eigen-decomposition of the integral operator with kernel \(k\), and the full proof is in Schölkopf and Smola's Learning with Kernels. The practical value of the condition is that it lets kernels be built compositionally without ever writing down \(\phi\). Sums, products, and positive scalar multiples of kernels are kernels, and \(k(x,x')g(x)g(x')\) is a kernel for any function \(g\).

The polynomial kernel, expanded explicitly

Take \(x, z \in \R^2\) and \(k(x,z) = (x\T z + c)^2\). Expanding,

Every term is a product of something depending only on \(x\) and something depending only on \(z\), so the feature map can be read off directly.

and one checks \(\phi(x)\T\phi(z) = (x\T z + c)^2\) term by term. Numerically, with \(x = (1,2)\), \(z = (3,-1)\), and \(c = 1\), the direct evaluation is \((1\cdot3 + 2\cdot(-1) + 1)^2 = 2^2 = 4\), and the feature map gives \(\phi(x) = (1, 4, 2\sqrt2, \sqrt2, 2\sqrt2, 1)\), \(\phi(z) = (9, 1, -3\sqrt2, 3\sqrt2, -\sqrt2, 1)\), whose inner product is \(9 + 4 - 12 + 6 - 8 + 1 = 4\). They agree exactly.

The counting generalizes. \((x\T z + c)^p\) with \(x \in \R^d\) corresponds to a feature map into \(\binom{d+p}{p}\) dimensions, all monomials up to degree \(p\) with binomial weights. At \(d = 100\), \(p = 5\) that is 96,560,646 features, computed by the kernel at the cost of one dot product in 100 dimensions and one exponentiation. That ratio is the entire economic argument for kernels.

The RBF kernel is an infinite-dimensional map

The Gaussian or RBF kernel is \(k(x,z) = \exp(-\gamma\lVert x - z\rVert^2)\). Expanding the square, \(\lVert x-z\rVert^2 = \lVert x\rVert^2 + \lVert z\rVert^2 - 2x\T z\), so

Each term \((x\T z)^m\) is the polynomial kernel of degree \(m\), which has an explicit finite feature map as above. The RBF kernel is therefore a convergent infinite weighted sum of polynomial kernels of every degree, so its feature map is the concatenation of all those monomial maps with weights \(\sqrt{(2\gamma)^m/m!}\) times the scalar \(e^{-\gamma\lVert x\rVert^2}\), genuinely infinite-dimensional. This makes the representer theorem indispensable rather than merely convenient, since the primal parameter vector does not exist as a finite object. It also explains the role of \(\gamma\). Large \(\gamma\) puts weight on high-degree terms, making \(k(x,z)\) decay fast with distance and the model very local (low bias, high variance), while small \(\gamma\) makes every pair of points look similar and the model nearly linear.

Random features, undoing the kernel trick on purpose

Rahimi and Recht (2007, then at Intel Research Berkeley) observed that the \(O(n^2)\) Gram matrix can be avoided by going in the opposite direction. Instead of computing an infinite-dimensional inner product implicitly, approximate it with an explicit low-dimensional random feature map, and then run ordinary linear methods.

The derivation rests on Bochner's theorem, which says a continuous shift-invariant kernel \(k(x,z) = \kappa(x-z)\) with \(\kappa(0) = 1\) is positive definite if and only if \(\kappa\) is the Fourier transform of a probability measure \(p(\omega)\). Then

For real-valued kernels the imaginary parts cancel, and a standard trigonometric rewriting removes the complex arithmetic. Using \(\cos(a - b) = 2\,\E_{b \sim U[0,2\pi]}[\cos(a+b)\cos(b'+b)]\) in the appropriate form, define \(z_\omega(x) = \sqrt{2}\cos(\omega\T x + b)\) with \(b \sim \mathrm{Uniform}[0, 2\pi]\). Then \(\E_{\omega, b}[z_\omega(x) z_\omega(z)] = \kappa(x - z) = k(x,z)\). Drawing \(D\) independent pairs \((\omega_m, b_m)\) and stacking gives the feature map

The estimator is an average of \(D\) i.i.d. bounded terms, so Hoeffding gives \(|\psi(x)\T\psi(z) - k(x,z)| = O_p(D^{-1/2})\) pointwise, with a uniform bound over a compact domain costing only a \(\sqrt{\log}\) factor. For the RBF kernel with \(k(x,z) = \exp(-\gamma\lVert x-z\rVert^2)\), the corresponding spectral density is Gaussian, \(\omega \sim \N(0, 2\gamma I)\).

The \(D^{-1/2}\) rate is exactly what the measurement shows. Approximating an RBF kernel with \(\gamma = 0.5\) on 300 points in \(\R^8\) gives the following.

The RMS error tracks \(1/\sqrt{D}\) to two significant figures at every size, halving each time \(D\) quadruples. Downstream, kernel ridge regression run on the random features instead of the exact kernel converges to the exact solution at the same rate, with RMS differences in test predictions of 0.245, 0.093, 0.065, and 0.031 at \(D = 128, 512, 2048, 8192\). The engineering payoff is the cost model. Exact kernel ridge is \(O(n^2 d + n^3)\) and \(O(n^2)\) memory, random features is \(O(nDd + D^3)\) and \(O(nD)\) memory, so at \(n = 10^6\) and \(D = 10^4\) the second is feasible and the first is not. Random features are also the conceptual bridge to wide neural networks, since a one-hidden-layer network with frozen random first-layer weights is a random-feature model, which is the starting point for the neural tangent kernel analysis of Jacot, Gabriel, and Hongler at EPFL.

Support vector machines

The geometric margin

Use labels \(y \in \{-1, +1\}\) and a hyperplane \(\{x : w\T x + b = 0\}\). The functional margin of a point is \(\hat\gamma_i = y_i(w\T x_i + b)\), which is positive when the point is on the correct side. It is not a distance, because scaling \((w,b)\) by any positive constant scales it without moving the hyperplane. The geometric margin fixes that. The unit normal to the hyperplane is \(w/\lVert w\rVert\), so the point \(x_i - \gamma_i\, y_i\, w/\lVert w\rVert\) obtained by walking from \(x_i\) toward the hyperplane by distance \(\gamma_i\) lies on it,

which is invariant to rescaling and is a genuine Euclidean distance. The margin of the dataset is \(\gamma = \min_i \gamma_i\), and the maximum-margin classifier maximizes it.

The primal problem

Maximizing \(\min_i \hat\gamma_i / \lVert w\rVert\) is awkward because of the scale freedom, so fix it by requiring the closest points to have functional margin exactly 1. Then the geometric margin is \(1/\lVert w\rVert\), maximizing it is minimizing \(\lVert w\rVert\), and the problem becomes a convex quadratic program.

The objective is convex quadratic and the constraints are affine, so this is a well-conditioned problem with a unique solution when the data is separable, and Slater's condition holds whenever a strictly feasible point exists, giving strong duality.

The Lagrangian dual, with KKT

Introduce multipliers \(\alpha_i \ge 0\) for the constraints, written as \(1 - y_i(w\T x_i + b) \le 0\).

Minimize over the primal variables by setting derivatives to zero.

The first is the representer theorem again, arrived at from a completely different direction. Substitute both back into \(\mathcal{L}\). The quadratic term becomes \(\tfrac12\sum_{i,j}\alpha_i\alpha_j y_iy_j x_i\T x_j\), the term \(\sum_i\alpha_i y_i w\T x_i\) becomes \(\sum_{i,j}\alpha_i\alpha_j y_iy_j x_i\T x_j\), the \(b\) term vanishes by the second condition, and \(\sum_i \alpha_i\) survives. Combining the two quadratic pieces leaves

The data appears only through inner products, so every kernel from the previous section applies verbatim. The remaining KKT conditions are stationarity (used above), primal and dual feasibility, and complementary slackness,

This single equation gives the support-vector characterization. If a point is strictly inside the correct side, \(y_i(w\T x_i + b) > 1\), the bracket is nonzero so \(\alpha_i = 0\) and the point contributes nothing to \(w = \sum_i \alpha_i y_i x_i\). Only points sitting exactly on the margin, \(y_i(w\T x_i+b) = 1\), can have \(\alpha_i > 0\). The solution depends only on the support vectors. Deleting every other training point and refitting gives the identical hyperplane. That is a sparsity property of the solution, not of the algorithm, and it is what makes kernel SVMs usable at prediction time even when \(n\) is large, provided the number of support vectors is not.

Soft margin, C, and the hinge loss

Real data is not separable, and one mislabeled point makes the hard-margin problem infeasible. Introduce slack variables \(\xi_i \ge 0\) that allow constraint violations and pay for them linearly.

Running the same Lagrangian computation, with multipliers \(\alpha_i\) for the margin constraints and \(\mu_i\) for \(\xi_i \ge 0\), the stationarity condition in \(\xi_i\) is \(C - \alpha_i - \mu_i = 0\), and \(\mu_i \ge 0\) then forces \(\alpha_i \le C\). The dual is identical to the separable case except that the constraint \(\alpha_i \ge 0\) becomes the box constraint \(0 \le \alpha_i \le C\). That is the whole change, and it gives \(C\) its meaning as an upper bound on how much influence any single training point can have. Small \(C\) means a wide margin tolerating many violations, while large \(C\) approaches the hard-margin problem. Complementary slackness now classifies points into three groups, \(\alpha_i = 0\) (outside the margin, ignored), \(0 < \alpha_i < C\) (exactly on the margin, \(\xi_i = 0\)), and \(\alpha_i = C\) (inside the margin or misclassified, \(\xi_i > 0\)).

Measured on 400 points of overlapping Gaussians, sweeping \(C\) traces this out precisely. At \(C = 0.01\) the margin is 1.300 with 255 support vectors, at \(C = 1\) it is 0.818 with 184, and at \(C = 100\) it is 0.797 with 182, having essentially converged to the hard-margin limit. The margin shrinks and the support set thins as \(C\) grows, exactly as the theory says.

The unconstrained form. At the optimum, \(\xi_i\) is as small as the constraints allow, so \(\xi_i = \max(0,\, 1 - y_i(w\T x_i + b))\). Substituting eliminates both the slacks and the constraints.

This is regularized empirical risk minimization with the hinge loss \(\ell(m) = (1-m)_+\) as a function of the margin \(m = y f(x)\), and it puts the SVM in exactly the same frame as everything else on this page. Compare the three margin losses. Hinge \((1-m)_+\) is zero once the margin exceeds 1, which is where sparsity comes from. Logistic \(\log(1 + e^{-m})\) is never zero, so every point pulls forever, which is why logistic regression has no support vectors. Exponential \(e^{-m}\), used by AdaBoost, punishes negative margins hardest of the three and is correspondingly least robust to label noise. All three are convex upper bounds on the 0-1 loss \(\mathbb{1}[m \le 0]\), which is what makes them tractable surrogates, and all three are classification-calibrated in the sense of Bartlett, Jordan, and McAuliffe at Berkeley, meaning that driving the surrogate risk to its minimum drives the 0-1 risk to the Bayes rate.

Because the hinge form is unconstrained, it can be attacked with subgradient descent, and for large \(n\) that is often better than solving the dual. As verification, Pegasos-style stochastic subgradient descent on the hinge objective, run for 200,000 single-example steps on \(n = 800\) points, reached a primal objective of 178.2495 against libsvm's dual solution at 178.1805, a relative gap of \(3.9 \times 10^{-4}\), with the two weight vectors at cosine similarity 0.999997. First-order methods do not need the dual.

SMO in one paragraph

The dual is a quadratic program with \(n\) variables, box constraints, and one linear equality constraint \(\sum_i\alpha_iy_i = 0\). General QP solvers need the full \(n \times n\) kernel matrix in memory. Platt's sequential minimal optimization (1998) exploits the observation that the equality constraint makes it impossible to change one \(\alpha_i\) alone, but a pair \((\alpha_i, \alpha_j)\) can be optimized exactly in closed form. With all others fixed, \(\alpha_i y_i + \alpha_j y_j\) is a constant, so \(\alpha_j\) determines \(\alpha_i\), the objective becomes a one-dimensional quadratic, and its unconstrained optimum is clipped to the box. Each step is arithmetic with no inner solver, and the algorithm's intelligence lives entirely in the heuristics for choosing which pair to update, which pick the pair that most violates the KKT conditions. The result needs \(O(n)\) memory and is what libsvm, and therefore scikit-learn's SVC, still runs today.

Model selection

The bias-variance decomposition, exactly

Fix a query point \(x\). Data are generated as \(y = f(x) + \varepsilon\) with \(\E[\varepsilon] = 0\), \(\Var[\varepsilon] = \sigma^2\), and a learning procedure trained on a random dataset \(\D\) produces a predictor \(\hat h_\D\). The expected squared error of the prediction at \(x\), over both the noise in the test label and the randomness of the training set, decomposes with no approximation. Write \(\bar h(x) = \E_\D[\hat h_\D(x)]\) for the average prediction over training sets. First split off the test noise.

because \(\varepsilon\) is independent of \(\D\) and has mean zero, killing the cross term \(2\,\E[\varepsilon]\,\E_\D[f - \hat h_\D] = 0\). Now insert and subtract \(\bar h(x)\) inside the remaining square.

The last expectation is \(\bar h - \bar h = 0\) by the definition of \(\bar h\), so the cross term vanishes identically and

Nothing here is an inequality or an asymptotic statement. For squared loss the decomposition is an identity, made famous for neural networks by Geman, Bienenstock, and Doursat (1992). Bias is the systematic error of the procedure's average answer, and variance is how much the answer sloshes across training sets. Two warnings usually omitted. The decomposition is a property of squared loss and does not carry over cleanly to 0-1 loss, where the analogues are messy and the "variance" term can be negative in the usual constructions. And it is a statement about a fixed learning procedure averaged over datasets, not about a single fitted model, so "this model has high variance" is a category error unless the procedure is meant.

Measured, to see the identity hold and the U-curve appear. Target \(f(x) = \sin(2\pi x)\) on \([0,1]\), \(n = 25\) points, noise \(\sigma = 0.35\) so \(\sigma^2 = 0.1225\), 800 replicate datasets, polynomial least squares of varying degree, evaluated at \(x = 0.35\).

Bias falls by a factor of 74 from degree 0 to degree 3 and is numerically zero by degree 5 (a degree-5 polynomial approximates \(\sin\) on \([0,1]\) to within the noise). Variance rises monotonically by a factor of 16 across the table, and their sum plus the noise is minimized in the middle, at degree 5, with total 0.1440 against the irreducible floor of 0.1225. The last two columns agree to within Monte-Carlo error at every row, confirming the identity numerically. This is the classical U-curve, and the last part of this section explains why it is not the whole story.

Cross-validation, and the leave-one-out shortcut

\(k\)-fold cross-validation splits the data into \(k\) parts, trains on \(k-1\) and evaluates on the held-out part, and averages. What it estimates is the expected test error of the learning procedure trained on \(n(k-1)/k\) points, not of the specific model fit on all \(n\). That mismatch is the source of its bias, since smaller \(k\) means smaller training sets, hence a pessimistic estimate. It was measured on ridge regression with \(n = 80\), \(d = 20\), 200 replicate datasets, against the true expected test error of the same procedure trained on all 80 points (0.33056).

The bias story is unambiguous and matches the theory. Two-fold CV overestimates the error by 50% relative, and the bias shrinks monotonically to under 3% relative at LOO. The variance story deserves care, because the folk claim "LOOCV has high variance" is stated far more confidently than the evidence supports. Two different variances are being conflated. Across datasets (fourth column), the LOO estimate here has the lowest spread, not the highest. Conditional on a fixed dataset, over the randomness of the fold assignment (fifth column), the spread falls monotonically to exactly zero at LOO, because LOO is deterministic, there being only one way to leave one out. The real argument against LOO is different and narrower. The \(n\) training sets overlap in \(n-2\) points, so the \(n\) fold errors are highly correlated, and their average is close to an average of one thing, which limits how much averaging can reduce the estimator's variance as an estimate of the risk. That is a statement about correlation, not about a large observed spread, and it is why \(k = 5\) or \(10\) remains the default, being cheaper and empirically about as accurate.

For least squares LOO is free, which is worth knowing. Let \(\hat\theta^{(-i)}\) be the fit with point \(i\) deleted. A rank-one update of \(X\T X\) via the Sherman-Morrison identity gives

so the entire leave-one-out error is a byproduct of the single full fit.

High-leverage points (\(h_{ii}\) near 1) have their residuals inflated the most, which is exactly right, since those are the points the model fitted itself to. Problem 1 verified the identity numerically to below \(10^{-12}\). Generalized cross-validation replaces each \(h_{ii}\) by their average \(d/n\), giving \(\mathrm{GCV} = \frac{\mathrm{RSS}/n}{(1-d/n)^2}\), which is what smoothing-spline and ridge software uses when the leverages are expensive.

AIC and BIC

Information criteria trade fit against complexity without resampling. Both have the form \(-2\log\hat{L} + \text{penalty}\), where \(\hat L\) is the maximized likelihood. AIC (Akaike, 1974) uses penalty \(2k\) for \(k\) parameters. It comes from an asymptotic bias correction, since the maximized log-likelihood overestimates the expected log-likelihood on new data by approximately \(k\), so subtracting \(k\) (doubled by convention) is an unbiased-ish estimate of out-of-sample deviance. AIC is therefore an estimate of predictive performance and is not consistent for model identification. It keeps a nonzero probability of overfitting even as \(n \to \infty\). BIC (Schwarz, 1978) uses penalty \(k\log n\). It comes from a Laplace approximation to the marginal likelihood \(p(\D \mid \text{model})\), so it estimates posterior model probability under a unit-information prior, and it is consistent when the true model is in the candidate set. Since \(\log n > 2\) for \(n > 7\), BIC always penalizes more and selects smaller models.

Here are both criteria on a measured example, polynomial regression where the truth is cubic, \(n = 120\), \(\sigma = 0.6\), so \(\log n = 4.79\).

Both pick degree 3, correctly. Notice how little the RSS improves past degree 3 (38.47 to 37.95, a 1.4% reduction for six extra parameters) and how sharply BIC punishes it, rising 27 points where AIC rises only 10. Two cautions apply in practice. \(k\) must be the effective degrees of freedom, not the nominal parameter count, for any penalized fit, and both criteria assume a correctly specified likelihood, so neither is meaningful for a model fitted by an objective that is not a log-likelihood.

Regularization paths

Choosing \(\lambda\) means fitting the model at many \(\lambda\), so it pays to compute the whole path at once. Ridge gets it almost free. After one SVD of \(X\), every \(\theta_\lambda = V\diag(d_i/(d_i^2+\lambda))U\T y\) costs \(O(d)\), so a hundred values of \(\lambda\) cost barely more than one. The lasso path is piecewise linear in \(\lambda\), which the LARS algorithm of Efron, Hastie, Johnstone, and Tibshirani (2004) exploits to compute all of the knots exactly. In practice glmnet-style coordinate descent over a decreasing grid with warm starts is faster and is what everything uses, since the solution at \(\lambda_{t}\) is an excellent initialization for \(\lambda_{t+1}\). The grid convention worth copying is to start at \(\lambda_{\max} = \max_j |x_j\T y|\), which by the KKT conditions is the smallest penalty that zeroes every coefficient, and descend geometrically to \(\lambda_{\max}/1000\).

Double descent, the modern revision of the U-curve

The classical picture says test error is U-shaped in model complexity and that interpolating the training data is a disaster. Modern practice contradicts this daily. Networks with far more parameters than samples, trained to zero training error, generalize well. Belkin, Hsu, Ma, and Mandal (2019) reconciled the two with the double-descent curve. As complexity increases past the interpolation threshold, where the model has just enough capacity to fit the training data exactly, test error spikes, and then falls again, often below the classical minimum. Nakkiran and coauthors (2020) demonstrated the same shape for deep networks in model size, in training epochs, and in dataset size, the last being the counterintuitive one, where more data can hurt at fixed model size.

The mechanism is visible in linear algebra, so it was measured on random ReLU features with \(n = 60\) training points and minimum-norm interpolation (the pseudoinverse solution) as \(p\), the number of features, sweeps past \(n\).

The peak is real and large. Test error rises from 0.609 at \(p = 40\) to 16.2 at \(p = 65\), a factor of 27, and then descends to 0.101 at \(p = 1000\), six times better than the best underparameterized model. The last row explains why. Near \(p = n\) there is exactly one interpolating solution and it has no freedom to be small, so the coefficient norm explodes to 10.8. Past \(p = n\) there are infinitely many interpolating solutions and the minimum-norm one gets steadily smaller, down to 0.19 at \(p = 3000\). The variance spike at the threshold is the smallest singular value of the feature matrix passing through zero, and the second descent is the implicit regularization of the minimum-norm choice taking over.

The row that matters most in practice is the third. A tiny explicit ridge penalty, \(\lambda = 10^{-2}\), cuts the peak from 6.919 to 2.167 at \(p = n\) and from 16.196 to 2.851 at \(p = 65\), while changing nothing at all away from the threshold (0.307, 0.158, 0.101 in both rows). Double descent is largely an artifact of not regularizing. The peak is where the estimator is most ill-conditioned, which is exactly where a penalty helps most, and with optimally tuned ridge the curve is monotone. This is the position argued by Nakkiran and coauthors and consistent with the benign-overfitting analyses of Bartlett, Long, Lugosi, and Tsigler, which characterize precisely when the minimum-norm interpolant generalizes. It requires the covariance spectrum to have many small directions that can absorb the noise harmlessly. The correct summary is not "overfitting is fine" but "the classical U-curve was drawn for a parameter count that stops being the right complexity measure once models interpolate, and the effective complexity of a minimum-norm interpolant decreases as you add parameters."

Unsupervised learning

k-means as alternating minimization

The k-means objective over assignments \(c_i \in \{1,\dots,K\}\) and centroids \(\mu_k\) is

Minimizing jointly is NP-hard, but minimizing over each block with the other fixed is trivial. With \(\mu\) fixed, \(J\) separates over \(i\) and the best \(c_i\) is the nearest centroid, \(c_i = \argmin_k \lVert x_i - \mu_k\rVert^2\). With \(c\) fixed, \(J\) separates over \(k\) and \(\partial J/\partial \mu_k = -2\sum_{i: c_i = k}(x_i - \mu_k) = 0\) gives the cluster mean \(\mu_k = \frac{1}{|C_k|}\sum_{i \in C_k} x_i\). Lloyd's algorithm alternates the two.

Convergence. Each step minimizes \(J\) over one block holding the other fixed, so \(J\) is non-increasing at every half-step and the sequence \(J_0 \ge J_1 \ge \cdots \ge 0\) is monotone and bounded below, hence convergent. More is true. The assignment \(c\) takes finitely many values (\(K^n\)), and given \(c\) the optimal \(\mu\) is determined, so the pair \((c,\mu)\) visits finitely many configurations. Since \(J\) strictly decreases whenever \(c\) changes, no configuration can repeat, and the algorithm must halt in finitely many steps. What it halts at is a local minimum, and nothing in the argument says which one.

How much that matters was measured on six Gaussian clusters with standard deviation 0.45, 120 points each. The objective is monotone (\(849.3 \to 436.9 \to 352.0 \to 318.4 \to 301.9 \to \cdots \to 273.4\), decreasing at both the assignment and the update half-step), the best objective found is 273.351, and the worst random restart lands at 388.897, 42% worse. Only 29.5% of 200 random restarts reach the global optimum. Careful seeding helps less than its reputation suggests here. k-means++ style seeding reached the optimum on 31.5% of runs, a marginal improvement, though its real value is the \(O(\log K)\) approximation guarantee of Arthur and Vassilvitskii, which bounds the damage of a bad run rather than making good runs more common. The practical rule is unchanged. Always use multiple restarts and keep the best, which is why n_init defaults to 10 in every implementation.

Gaussian mixtures and EM, from the lower bound

A Gaussian mixture models the density as \(p(x) = \sum_{k=1}^{K}\pi_k\, \N(x;\mu_k,\Sigma_k)\). The log-likelihood

has a logarithm of a sum, which does not separate and has no closed-form maximizer. EM gets around this by optimizing a lower bound instead, and the derivation via that bound is worth doing in full because it generalizes far beyond mixtures.

Introduce for each point an arbitrary distribution \(q_i\) over the latent assignment \(z_i \in \{1,\dots,K\}\), and multiply and divide inside the sum,

by Jensen's inequality, since \(\log\) is concave and \(q_i\) is a probability distribution, so the expression is \(\log \E_{q_i}[\cdot] \ge \E_{q_i}[\log \cdot]\). The bound \(\mathcal{F}\) is the ELBO. Jensen is tight exactly when the argument is constant over the support, that is when \(p(x_i, z_i{=}k)/q_i(k)\) does not depend on \(k\), which forces \(q_i(k) \propto p(x_i, z_i{=}k)\), and since \(q_i\) must normalize,

the posterior over the latent variable, called the responsibility. The two steps follow.

• E step. Set \(q_i = p(z_i \mid x_i; \vartheta_t)\), which makes the bound touch the log-likelihood, \(\mathcal{F}(q, \vartheta_t) = \ell(\vartheta_t)\). Concretely \(r_{ik} = \pi_k \N(x_i;\mu_k,\Sigma_k) / \sum_m \pi_m \N(x_i;\mu_m,\Sigma_m)\).
• M step. Maximize \(\mathcal{F}(q, \vartheta)\) over \(\vartheta\) with \(q\) held fixed. Since \(\mathcal{F} = \sum_i \sum_k r_{ik}\log p(x_i, z_i{=}k;\vartheta) + \text{const}(q)\), this is a weighted complete-data maximum likelihood problem, which for Gaussians is solved in closed form.

The monotonic improvement proof is now three inequalities,

where (1) is the tightness of Jensen at the E-step choice of \(q\), (2) is the definition of the M step as a maximization over \(\vartheta\), and (3) is Jensen's inequality again, now at the new parameters where the bound need not be tight. Therefore \(\ell(\vartheta_{t+1}) \ge \ell(\vartheta_t)\), and EM never decreases the log-likelihood. Since the likelihood is bounded above for a well-posed model, the sequence converges. It converges to a stationary point of \(\ell\), not necessarily a global maximum, and the same restart advice as k-means applies.

Carry out the M step for Gaussians. With weights \(r_{ik}\) and \(N_k = \sum_i r_{ik}\), maximizing \(\sum_{i,k} r_{ik}[\log\pi_k + \log\N(x_i;\mu_k,\Sigma_k)]\) subject to \(\sum_k \pi_k = 1\) gives, by exactly the same three derivative computations used for GDA but with soft counts in place of hard ones,

The parallel with GDA is exact and worth stating. EM for a Gaussian mixture is GDA with the hard labels replaced by soft posterior responsibilities, iterated because the responsibilities depend on the parameters they are used to estimate. Setting each \(r_{ik}\) to 0 or 1 by its argmax gives "hard EM", and further fixing \(\Sigma_k = \sigma^2 I\) with \(\sigma \to 0\) recovers k-means exactly, which is why k-means is the right mental picture and the wrong model whenever clusters have different shapes or sizes.

Numerical stability is not optional here. Responsibilities are ratios of exponentially small numbers, and computing them by exponentiating first underflows. In the measured run a test point far from every component has log joint densities \((-4286.8, -5838.3, -12129.9)\). Exponentiating gives \((0, 0, 0)\) in float64 and the normalization is \(0/0\). Computing \(r_{ik} = \exp(\log p_{ik} - \mathrm{logsumexp}_k \log p_{ik})\) instead returns \((1, 0, 0)\) correctly, because the subtraction happens in log space where the numbers are representable. Every serious implementation does this, and the covariance update needs a floor \(\Sigma_k \mathrel{+}= \epsilon I\) as well, since a component that captures a single point drives its covariance to zero and the likelihood to \(+\infty\), a genuine singularity of the objective rather than a numerical accident.

The measured run used 1500 points, 3 components, and deterministic initialization. EM converged in 63 iterations, the mean log-likelihood rose monotonically from \(-6.7193\) to \(-3.81889\) with the smallest single-iteration increment \(+8.8\times10^{-13}\) (never negative), and the converged value matches scikit-learn's GaussianMixture from the same initialization to \(4.4\times10^{-16}\), with identical mixing weights \((0.2447, 0.3539, 0.4014)\) to six decimals.

EM and variational inference are the same algorithm

The lower-bound derivation above never used the fact that \(q_i\) was chosen to be the exact posterior. It only used that \(q_i\) is a distribution. Writing the gap explicitly,

which is an identity, obtained by expanding \(\mathcal{F} = \E_q[\log p(x,z)] + \mathbb{H}[q]\) and \(p(x,z) = p(z\mid x)p(x)\). This rewrites the whole algorithm as coordinate ascent on \(\mathcal{F}(q,\vartheta)\). The E step maximizes over \(q\) (achieving \(\KL = 0\) by setting \(q = p(z\mid x)\)) and the M step maximizes over \(\vartheta\). When the exact posterior is intractable, restrict \(q\) to a family \(\mathcal{Q}\) (factorized, or parameterized by a neural network) and maximize \(\mathcal{F}\) over that family instead. The bound no longer touches, so the E step becomes approximate, but the M step and the monotonicity argument for \(\mathcal{F}\) survive untouched. That is variational EM, and it is exactly the ELBO that a variational autoencoder maximizes, with \(q\) an amortized encoder network and the M step folded into the same gradient update. Mean-field variational inference, stochastic variational inference, and the VAE are all this one diagram with different choices of \(\mathcal{Q}\) and different optimizers.

PCA, derived twice

Center the data so \(\sum_i x_i = 0\), and let \(S = \frac{1}{n}\sum_i x_i x_i\T\) be the sample covariance.

Derivation 1, maximum variance. Find the unit direction \(u\) along which the projected data varies most. The projection of \(x_i\) onto \(u\) is \(u\T x_i\), and its sample variance is

so the problem is \(\max_{\lVert u\rVert = 1} u\T S u\), the maximization of a Rayleigh quotient. Form the Lagrangian \(u\T S u - \lambda(u\T u - 1)\) and set the gradient to zero.

Stationary points are eigenvectors, and at an eigenvector the objective value is \(u\T S u = \lambda\), so the maximum is attained at the eigenvector of the largest eigenvalue. Subsequent directions come from the same problem restricted to the orthogonal complement, yielding the top \(k\) eigenvectors.

Derivation 2, minimum reconstruction error. Choose an orthonormal \(U_k \in \R^{d\times k}\) and represent each point by its projection \(\hat x_i = U_kU_k\T x_i\). Minimize the reconstruction error.

where the middle and last terms combined to \(-\lVert U_k\T x_i\rVert^2\) and the sum over \(i\) turned it into a trace against \(S\). The first term does not depend on \(U_k\). Therefore minimizing reconstruction error is exactly maximizing retained variance \(\tr(U_k\T S U_k)\), and the two derivations produce the same \(U_k\). The identity behind the equivalence is Pythagorean, total variance = retained variance + reconstruction error, which was checked numerically and holds to \(3\times10^{-14}\) at every \(k\) from 1 to 10 in the measured run.

Via the SVD, which is how it is computed. Never form \(S\), since it squares the condition number for the same reason the normal equations do. Instead take the thin SVD of the centered data matrix, \(X_c = U D V\T\). Then \(S = \frac{1}{n}X_c\T X_c = \frac{1}{n} V D^2 V\T\), so the principal directions are the right singular vectors \(V\), the eigenvalues are \(d_i^2/n\), and the scores are \(X_c V = UD\). Verified against scikit-learn's PCA, the components are identical to \(0.0\) after sign alignment, with explained-variance ratios agreeing to \(1.1\times10^{-16}\).

Probabilistic PCA. Tipping and Bishop (1999) gave PCA a generative model, \(z \sim \N(0, I_q)\) and \(x \mid z \sim \N(Wz + \mu, \sigma^2 I_d)\), which marginalizes to \(x \sim \N(\mu, WW\T + \sigma^2 I)\). Maximizing the marginal likelihood has a closed form. With \(\lambda_1 \ge \cdots \ge \lambda_d\) the eigenvalues of \(S\),

for any orthogonal \(R\) (the model is identified only up to rotation of the latent space). The noise variance is the average discarded eigenvalue, which is a satisfying statement. The variance PCA throws away is exactly what the probabilistic model calls noise. Measured with true injected noise variance \(0.4^2 = 0.16\), the estimate is \(0.165103\), and the recovered subspace is within a principal angle of \(0.99^\circ\) of the true one. The probabilistic version buys a likelihood (so PCA can be compared to other models, and \(q\) chosen by BIC), principled handling of missing data through EM, and mixtures of PPCA for nonlinear structure. As \(\sigma^2 \to 0\) it recovers classical PCA exactly.

ICA, when uncorrelated is not enough

PCA finds uncorrelated directions. Independent component analysis finds independent ones, which is strictly stronger and solves a different problem. The model has sources \(s \in \R^d\) with independent components, observed after an invertible linear mixing \(x = As\). Let \(W = A^{-1}\), so \(s = Wx\). The change-of-variables formula for densities gives

using independence to factor \(p_s\). The log-likelihood over the data is therefore

Differentiating needs \(\nabla_W \log|\det W| = W^{-\mathsf{T}}\), which follows from Jacobi's formula. Choosing a source density with heavier-than-Gaussian tails, the standard convenient choice being the logistic density whose log-derivative is \(1 - 2\sigma(u)\), gives the update

Why Gaussian sources are impossible. If \(s \sim \N(0, I)\) then for any orthogonal \(R\), \(Rs\) has the same distribution, so \(x = As\) and \(x = (AR\T)(Rs)\) are observationally identical and \(A\) cannot be recovered beyond a rotation. Non-Gaussianity is not a technical convenience, it is the identifying assumption, and this is why ICA is sometimes described as a search for maximally non-Gaussian projections.

Measured on two mixed non-Gaussian signals (a square wave and an amplitude-modulated sine) through the mixing matrix \(\begin{pmatrix}1 & 0.7\\ 0.4 & 1\end{pmatrix}\), ICA recovers the sources with mean absolute correlation 0.9937 (individual correlations 0.998 and 0.9894), while PCA on the same data reaches only 0.7431, because the principal axes are the wrong axes. With the same mixing applied to genuinely Gaussian sources, ICA's recovery falls to 0.9057, and what remains is the accidental alignment expected from a two-dimensional rotation rather than genuine identification, exactly as the theory predicts.

Trees and ensembles

Recursive partitioning and impurity

A decision tree partitions the input space into axis-aligned boxes and fits a constant in each. Finding the optimal partition is NP-hard, so trees are grown greedily. At each node, choose the feature and threshold that most reduces an impurity measure, then recurse. For classification with class proportions \(p_k\) at a node, the two standard impurities are

and a split into children \(L, R\) is scored by the weighted impurity drop \(\Delta = I(\text{parent}) - \frac{n_L}{n}I(L) - \frac{n_R}{n}I(R)\). Gini has an interpretation worth knowing. It is the probability of misclassifying a point drawn from the node if it were labelled by drawing from the node's own class distribution, and it is also the expected 0-1 loss of a randomized classifier. Entropy is the expected number of bits to encode the label. Both are concave in \(p\), which is what guarantees \(\Delta \ge 0\) for every split by Jensen's inequality, so a greedy step can never increase impurity. Their practical difference is negligible, though entropy is slightly more willing to isolate a small pure group because \(-p\log p\) has infinite slope at 0. For regression the impurity is the within-node variance and the same machinery applies.

A single tree is a high-variance, low-bias predictor. Grown to purity it interpolates the training data, and changing one point can change a high split and reshape everything below it. Pruning by cost-complexity, minimizing \(\text{RSS} + \alpha|T|\) over subtrees with \(\alpha\) chosen by cross-validation, controls this, but the far more effective answer is to average many trees.

The variance-reduction argument for bagging, and its limit

Suppose \(B\) predictors each have variance \(\sigma^2\) and pairwise correlation \(\rho\). The variance of their average is

The second term vanishes as \(B \to \infty\). The first does not. Averaging removes only the uncorrelated part of the variance, and \(\rho\sigma^2\) is a floor no amount of averaging can break. Breiman's bagging (1996) creates the ensemble by bootstrap resampling the training data, which decorrelates the trees a little because each sees a different sample (about 63.2% of the distinct points, since \(1 - (1-1/n)^n \to 1 - e^{-1}\)). The out-of-bag remainder gives a free validation estimate. Bias is essentially unchanged, since each tree is nearly unbiased and the average of unbiased predictors is unbiased.

Random forests (Breiman, 2001) attack \(\rho\) directly. At every split, consider only a random subset of \(m\) of the \(d\) features. This makes individual trees worse (their \(\sigma^2\) goes up, because they are sometimes forced to split on a mediocre feature) in exchange for making them less alike (their \(\rho\) goes down). The formula above says the trade is worth it whenever the reduction in \(\rho\sigma^2\) exceeds the increase in \(\sigma^2\), which is an empirical question with an optimum in the interior. It was measured on a 600-point regression with 5 informative and 5 pure-noise features.

The two competing effects are separately visible. \(\rho\) rises monotonically from 0.237 to 0.687 as \(m\) grows, and individual tree quality improves monotonically from 27.2 to 11.7. The ensemble error is minimized in the interior at \(m = 7\), where the product is best, and pure bagging (\(m = 10\)) is slightly worse. Note that a single unrestricted tree scores 11.27 on this problem and the ensemble scores 4.14, so averaging cut the error by 63% without touching the bias. The commonly cited defaults, \(m = \sqrt{d}\) for classification and \(d/3\) for regression, are reasonable starting points and, as this table shows, worth tuning, since the curve is flat near the optimum but steep at small \(m\).

Boosting as gradient descent in function space

Bagging reduces variance by averaging independently-built models. Boosting reduces bias by building models sequentially, each fixing what the previous ones got wrong. Friedman's (2001) formulation makes the connection to optimization exact.

Seek an additive model \(F(x) = \sum_{m=1}^{M}\nu\, f_m(x)\) minimizing \(\sum_i L(y_i, F(x_i))\). Treat the vector of training predictions \(\mathbf{F} = (F(x_1), \dots, F(x_n)) \in \R^n\) as the optimization variable. Gradient descent on the loss in this \(n\)-dimensional space would step along

but a step defined only at the training points is not a function and cannot be evaluated anywhere else. Boosting fixes this by projecting the gradient onto the class of base learners. Fit a regression tree \(f_m\) to the targets \(-g_i\) by least squares, which finds the base learner most aligned with the negative gradient, and take a step along it, \(F \leftarrow F + \nu f_m\). Boosting is gradient descent where the step direction is constrained to lie in the span of the base learners, and the learning rate \(\nu\) is the step size. For squared loss \(-g_i = y_i - F(x_i)\), the residual, which recovers the folk description of boosting as "fit the residuals" as the special case it is.

AdaBoost and the exponential loss. Freund and Schapire's AdaBoost (1997) predates this view and looks different. Reweight misclassified examples, fit a weak classifier, weight it by \(\alpha_m = \tfrac12\log\frac{1-\epsilon_m}{\epsilon_m}\) where \(\epsilon_m\) is its weighted error, repeat. Friedman, Hastie, and Tibshirani (2000) showed it is forward stagewise additive modelling with the exponential loss \(L(y, F) = e^{-yF}\) for \(y \in \{-1,+1\}\). For the derivation, at step \(m\), choose \((\alpha, f)\) to minimize \(\sum_i e^{-y_i(F_{m-1}(x_i) + \alpha f(x_i))} = \sum_i w_i^{(m)} e^{-y_i\alpha f(x_i)}\) with \(w_i^{(m)} = e^{-y_iF_{m-1}(x_i)}\). Splitting the sum by whether \(f(x_i) = y_i\) gives \(e^{-\alpha}\sum_{\text{correct}} w_i + e^{\alpha}\sum_{\text{wrong}} w_i\). For any fixed \(\alpha > 0\) this is minimized by the \(f\) with least weighted error, and differentiating in \(\alpha\) yields exactly Freund and Schapire's \(\alpha_m = \tfrac12\log\frac{1-\epsilon_m}{\epsilon_m}\). The exponentially growing weights on misclassified points are just \(w_i \propto e^{-y_iF(x_i)}\), the exponential loss evaluated at the current margin. Note that the population minimizer of exponential loss is \(F^\ast(x) = \tfrac12\log\frac{p(y=1\mid x)}{p(y=-1\mid x)}\), half the log-odds, so AdaBoost is estimating a logistic model through an unusual route.

Why boosting keeps improving after zero training error. This was the puzzle of the 1990s, and the resolution is the margin explanation of Schapire, Freund, Bartlett, and Lee (1998). The training error stops improving but the margin distribution \(y_iF(x_i)/\sum_m|\alpha_m|\) keeps shifting to the right, and generalization bounds that depend on the margin rather than on the number of rounds continue to tighten. The exponential loss keeps pushing on points that are already correct, unlike hinge loss, which is exactly why margins keep growing. The same loss is also why AdaBoost is fragile under label noise. A mislabeled point acquires exponentially large weight and the ensemble contorts itself around it. LogitBoost and gradient boosting with the logistic loss are the standard fix, and every modern implementation defaults to a logistic or Huberized loss for precisely this reason.

XGBoost, the second-order objective and closed-form leaf weights

Chen and Guestrin, then at the University of Washington, published XGBoost in 2016 with a refinement that turns out to matter a great deal in practice. Use a second-order Taylor expansion of the loss rather than a first-order one, and put an explicit regularizer on the tree. Write the objective at round \(t\), where \(f_t\) is the tree being added.

with \(T\) the number of leaves and \(w_j\) the value in leaf \(j\). Taylor-expand each loss term to second order around \(F^{(t-1)}(x_i)\), writing \(g_i = \partial_F L\) and \(h_i = \partial^2_F L\) evaluated there.

A tree is constant on each leaf, so \(f_t(x_i) = w_{q(i)}\) where \(q(i)\) is the leaf that \(x_i\) falls into. Group the sum by leaf, defining \(G_j = \sum_{i \in I_j} g_i\) and \(H_j = \sum_{i \in I_j} h_i\).

Each leaf now appears in its own independent one-dimensional quadratic. Differentiating, \(G_j + (H_j+\lambda)w_j = 0\), so the optimal leaf value and the resulting objective are

The second formula is a score for a tree structure, and it gives the split criterion by differencing. Splitting a leaf with statistics \((G, H)\) into \((G_L,H_L)\) and \((G_R,H_R)\) changes the objective by

and the \(-\gamma\) makes it possible for a split to have negative gain, which is pre-pruning built into the criterion rather than bolted on afterwards. Three things are better here than in first-order gradient boosting. The leaf value is the exact minimizer of the local quadratic model rather than a line-search approximation, so each round makes more progress. The Hessian \(H_j\) acts as a per-example confidence weight. For logistic loss \(h_i = p_i(1-p_i)\), so confidently-predicted examples contribute little to both the leaf value and the split score, which is the IRLS weighting from the GLM section reappearing inside a tree. And the regularizer is on the tree, so it prunes and shrinks simultaneously. For squared loss, \(g_i = F_i - y_i\) and \(h_i = 1\), so \(w_j^\ast = -G_j/(n_j + \lambda)\) reduces to the mean residual in the leaf, shrunk toward zero by \(\lambda\). Verified numerically, the closed form and the mean negative residual agree at \(-0.00140195\) versus \(-0.00140195\) when \(\lambda = 0\).

An implementation from scratch of exactly this, 200 rounds of depth-3 trees with \(\eta = 0.1\) and \(\lambda = 1\), reaches test MSE 2.1807 on a 2000-point held-out set against scikit-learn's first-order GradientBoostingRegressor at 2.2829 with identical hyperparameters, a 4.5% improvement from the second-order leaf weights alone, with training MSE falling 24.07, 3.99, 1.59, 0.81, 0.49 over rounds 1, 25, 50, 100, 200 and test MSE 21.27, 5.12, 3.01, 2.37, 2.18. The PyTorch and JAX versions in the implementation section reproduce 2.1814 and 2.1777, the small differences coming from split tie-breaking.

What the production systems add. XGBoost's approximate split finding buckets each feature into quantile bins so a split scan costs \(O(\#\text{bins})\) instead of \(O(n)\), and its sparsity-aware algorithm learns a default direction for missing values at each node rather than imputing. LightGBM, from Ke and coauthors at Microsoft Research (2017), adds histogram subtraction (a child's histogram is the parent's minus its sibling's, halving the work), leaf-wise rather than level-wise growth, gradient-based one-side sampling that keeps large-gradient examples and subsamples the rest, and exclusive feature bundling for sparse data. The result is typically 5-10x faster training at comparable accuracy. CatBoost, from Prokhorenkova and coauthors at Yandex (2018), targets a subtler bug. Naive target-based encoding of categorical features leaks the label of the example being encoded, and their ordered boosting scheme computes statistics using only examples that precede the current one in a random permutation, eliminating the resulting prediction shift. Categorical handling is why CatBoost often wins on data with high-cardinality categories.

Making a model usable, from calibration to imbalance

Calibration

A model is calibrated if among the cases where it says 0.7, about 70% are positive. This is a different property from accuracy and from ranking quality, and it is the property that matters whenever the output feeds a decision with an explicit cost, a threshold set by policy, or an expected-value computation. The standard summary is the expected calibration error. Bin the predictions, and

the sample-weighted average gap between the empirical frequency and the mean predicted probability within each bin. ECE is a diagnostic, not a loss. It can be made small by a constant predictor, so it should always be read alongside a proper scoring rule such as the Brier score \(\frac{1}{n}\sum_i (p_i - y_i)^2\) or the log loss, both of which are minimized only by the true conditional probability and both of which decompose into calibration plus refinement terms.

Two standard recalibration methods, both fit on a held-out calibration set, never on the training data. Platt scaling (1999) fits a one-dimensional logistic regression to the model's score, \(\hat p = \sigma(A s + B)\), two parameters, which cannot change the ranking at all (it is monotone) and therefore leaves AUC exactly unchanged while fixing the probability scale. Isotonic regression (Zadrozny and Elkan, 2002) fits the best non-decreasing step function by the pool-adjacent-violators algorithm, which is far more flexible, cannot change the ranking either, and needs more calibration data (a few thousand points) or it overfits.

Measured on a 20-feature synthetic classification problem with a 2000-point calibration set and a 2000-point test set, 15 bins.

Several things worth reading off this table. Isotonic regression beats Platt scaling on all three metrics here, as expected with 2000 calibration points. With 200 the ordering usually reverses. Passing an SVM margin through a sigmoid, which people do constantly, gives an ECE of 0.164, by far the worst number in the table, because the margin is on an arbitrary scale that has nothing to do with log-odds. Platt scaling cuts it by a factor of 4.5 for the cost of fitting two parameters. And logistic regression, which optimizes a proper scoring rule directly, is well calibrated out of the box at ECE 0.0334 while being the worst model in the table by Brier score, which is the cleanest possible demonstration that calibration and accuracy are orthogonal. Guo, Pleiss, Sun, and Weinberger at Cornell (2017) showed the modern version of this problem. Deep networks became substantially less calibrated as they got deeper and more accurate, and a single temperature parameter fitted on validation data, which is Platt scaling with \(B\) fixed at 0, fixes most of it.

Imbalance, and why ROC-AUC can mislead

When positives are rare, accuracy is useless (predict "negative" always and score 99%) and ROC-AUC is subtly misleading. The reason is structural rather than empirical. The ROC curve plots the true positive rate against the false positive rate, and the false positive rate has the negatives in its denominator. When negatives outnumber positives 99 to 1, moving from 1,000 false positives to 2,000 changes the FPR from 0.0101 to 0.0202, a visually invisible step on the ROC curve, while it changes precision from 0.5 to 0.33, a catastrophe for anyone who has to act on the flagged cases. The precision-recall curve puts predicted positives in the denominator instead and therefore stays sensitive exactly where it matters. The baselines also differ. A random classifier has ROC-AUC 0.5 regardless of prevalence, while its average precision equals the prevalence, so 0.01 on a 1% problem.

A worked example on 100,000 cases with 1,000 positives (1% prevalence) and two scorers. Classifier A concentrates its skill at the very top of the ranking, while classifier B lifts all positives uniformly.

The two metrics rank the classifiers oppositely, and neither is wrong. They answer different questions. If a fraud team can investigate 100 cases a day, A finds 100 real frauds and B finds 62, and A is obviously the better system, which is what average precision reports. If the goal is to eventually retrieve as many positives as possible with a large review budget, B is slightly better at the 10,000 mark, which is where B's ROC-AUC advantage comes from. The rule is short. Report ROC-AUC when the operating point is unknown and both errors matter symmetrically. Report precision-recall, average precision, and precision at the actual review budget when positives are rare and only the top of the ranking will ever be acted on. For a hand-checkable illustration of the same sensitivity, with 10 items and 2 positives ranked 1st and 3rd, ROC-AUC is 0.9375 and average precision \(\tfrac12(1 + \tfrac23) = 0.8333\). Move the positives to ranks 2 and 4 and ROC-AUC falls only to 0.8125 while average precision falls to \(\tfrac12(\tfrac12 + \tfrac12) = 0.5\), a 40% drop against ROC's 13%.

On handling imbalance itself, the honest summary is that the methods are less impressive than their literature. Class weights (equivalently, reweighting the loss) are principled and cost nothing. They change the implied decision threshold and, for a proper loss, are equivalent to changing the prior. Random undersampling of the majority throws away data but is fast and often works. SMOTE and its descendants synthesize minority examples by interpolating between neighbours. They help metrics that reward recall and frequently damage calibration, since the training prevalence no longer matches deployment. The approach that most often wins is duller. Train on the natural distribution with a proper loss, keep the model calibrated, and then choose the threshold to optimize the actual business objective on a validation set. Resampling changes the model to move the threshold. Setting the threshold moves the threshold.

Gradient-boosted trees still win on tabular data

The claim is repeatedly rediscovered and repeatedly true. Grinsztajn, Oyallon, and Varoquaux at INRIA (2022) benchmarked tree ensembles against tuned neural architectures across 45 tabular datasets and found trees ahead, and, more usefully, identified the three properties of tabular data responsible. Targets are often irregular and non-smooth functions of the inputs, where a network's smoothness bias hurts. Uninformative features are common, and trees ignore them almost for free while networks must learn to. And the data is not rotation-invariant, meaning individual columns carry meaning, whereas MLPs are rotationally equivariant in their first layer and so cannot exploit that. Shwartz-Ziv and Armon reached the same conclusion independently in 2022, and the Yandex group's own attempts to design better tabular networks (FT-Transformer and relatives) narrowed the gap without closing it.

Measured here, RMSE on held-out data across three problems of different character.

The pattern is exactly the one Grinsztajn and coauthors describe, and it is a warning against blanket claims in either direction. On a small near-linear problem, plain ridge beats everything, and reaching for a tree ensemble is an error. On a smooth analytic target with no junk features, the MLP matches gradient boosting (1.269 against 1.286), because the neural network's inductive bias is correct for that target. On the irregular problem, with step functions and 15 uninformative columns out of 20, gradient boosting reaches 0.541 against an irreducible noise floor of 0.500, essentially solving it, while the best of nine MLP configurations manages only 1.311, more than twice the error, and ridge is worse still. The gap is not small and it is not a tuning artifact. It is the smoothness bias and the rotation-invariance issue showing up exactly where the theory says they will.

Worked problems

Six problems are distributed through the derivations above, at the points where they do the most good. Two more are collected here. One closes a gap, proving the leave-one-out identity that the model-selection section used without justification, and one is a numerical computation of the boosting split criterion that a reader can check with a calculator.

Implementation

Six algorithms, each implemented from scratch in PyTorch and JAX and checked against an independent reference. Everything runs in float64, because several of these blocks are specifically about numerical accuracy and float32 would hide exactly the effects being demonstrated. The listings share the preamble shown in the first block, and each later block continues the same script. All of them were executed on an H100 80GB HBM3 and the outputs quoted are what they printed.

Ridge regression by QR, and the conditioning experiment

The first block is the numerical-linear-algebra argument in executable form. Designs with a prescribed condition number are built by choosing the singular values directly, the exact solution is known by construction, and the same problem is solved through the normal equations, through QR on the stacked ridge system, and through the SVD. The ridge solver itself is the stacked-matrix trick. Appending \(\sqrt{\lambda} I\) to \(X\) and zeros to \(y\) turns ridge into an ordinary least-squares problem, so one QR handles both. On this machine the PyTorch version prints relative errors of 3.9e-13, 2.8e-9, 9.6e-6, and 1.0e-1 for the normal equations at \(\kappa(X) = 10^2, 10^4, 10^6, 10^8\) against 9.5e-16, 7.1e-14, 9.6e-12, and 1.3e-9 for QR. The JAX version, with its own seeds, prints 4.7e-13, 5.4e-9, 2.3e-5, 3.1e-1 against 3.6e-15, 2.5e-13, 4.1e-12, 9.7e-10. Same conclusion, eight orders of magnitude apart at the hard end.

"""PyTorch blocks for the page. Every block must run and print sane numbers."""
import math
import torch

torch.set_default_dtype(torch.float64)
dev = "cuda" if torch.cuda.is_available() else "cpu"

def ridge_qr(X, y, lam):
    """Ridge by QR on the stacked system. X: (n, d), y: (n,) -> (d,)."""
    n, d = X.shape
    Xa = torch.cat([X, math.sqrt(lam) * torch.eye(d, dtype=X.dtype, device=X.device)], 0)
    ya = torch.cat([y, torch.zeros(d, dtype=y.dtype, device=y.device)])
    Q, R = torch.linalg.qr(Xa)                       # Q: (n+d, d), R: (d, d)
    return torch.linalg.solve_triangular(R, (Q.T @ ya).unsqueeze(1), upper=True).squeeze(1)

def ridge_cholesky(X, y, lam):
    d = X.shape[1]
    G = X.T @ X + lam * torch.eye(d, dtype=X.dtype, device=X.device)
    L = torch.linalg.cholesky(G)
    return torch.cholesky_solve((X.T @ y).unsqueeze(1), L).squeeze(1)

def design_with_condition(n, d, kappa, gen):
    U, _ = torch.linalg.qr(torch.randn(n, d, generator=gen, device=dev))
    V, _ = torch.linalg.qr(torch.randn(d, d, generator=gen, device=dev))
    s = torch.logspace(0, -math.log10(kappa), d, device=dev, dtype=torch.float64)
    return (U * s) @ V.T

g = torch.Generator(device=dev).manual_seed(7)
n, d = 400, 30
theta_star = torch.randn(d, generator=g, device=dev)
print(f"{'cond(X)':>9} {'normal eq':>12} {'QR':>12} {'SVD':>12}")
for e in [2, 4, 6, 8]:
    X = design_with_condition(n, d, 10.0 ** e, g)
    y = X @ theta_star
    rel = lambda t: ((t - theta_star).norm() / theta_star.norm()).item()
    t_ne = torch.linalg.solve(X.T @ X, X.T @ y)
    t_qr = ridge_qr(X, y, 0.0)
    U, S, Vh = torch.linalg.svd(X, full_matrices=False)
    t_svd = Vh.T @ ((U.T @ y) / S)
    print(f"     1e{e} {rel(t_ne):12.3e} {rel(t_qr):12.3e} {rel(t_svd):12.3e}")

# near-collinear design: QR and Cholesky must still agree once lambda > 0
g = torch.Generator(device=dev).manual_seed(1)
X = torch.randn(500, 40, generator=g, device=dev)
X[:, 5] = X[:, 4] + 1e-6 * torch.randn(500, generator=g, device=dev)
y = X @ torch.randn(40, generator=g, device=dev) + 0.5 * torch.randn(500, generator=g, device=dev)
for lam in [1e-3, 1.0, 100.0]:
    d = (ridge_qr(X, y, lam) - ridge_cholesky(X, y, lam)).abs().max().item()
    print(f"lam={lam:<8g} QR vs Cholesky max|diff| = {d:.3e}")

"""JAX blocks for the page."""
import functools, math
import jax, jax.numpy as jnp
from jax import random, jit, lax
from jax.scipy.special import logsumexp

jax.config.update("jax_enable_x64", True)

@functools.partial(jit, static_argnums=())
def ridge_qr(X, y, lam):
    """Ridge by QR on the stacked system. X: (n, d), y: (n,) -> (d,)."""
    d = X.shape[1]
    Xa = jnp.concatenate([X, jnp.sqrt(lam) * jnp.eye(d)], axis=0)   # (n+d, d)
    ya = jnp.concatenate([y, jnp.zeros(d)])
    Q, R = jnp.linalg.qr(Xa)                                        # Q:(n+d,d) R:(d,d)
    return jax.scipy.linalg.solve_triangular(R, Q.T @ ya, lower=False)

@jit
def ridge_cholesky(X, y, lam):
    d = X.shape[1]
    L = jnp.linalg.cholesky(X.T @ X + lam * jnp.eye(d))
    return jax.scipy.linalg.cho_solve((L, True), X.T @ y)

def design_with_condition(key, n, d, kappa):
    k1, k2 = random.split(key)
    U, _ = jnp.linalg.qr(random.normal(k1, (n, d)))
    V, _ = jnp.linalg.qr(random.normal(k2, (d, d)))
    s = jnp.logspace(0, -math.log10(kappa), d)
    return (U * s) @ V.T

key = random.PRNGKey(7)
key, sk = random.split(key)
theta_star = random.normal(sk, (30,))
print(f"{'cond(X)':>9} {'normal eq':>12} {'QR':>12} {'SVD':>12}")
for e in [2, 4, 6, 8]:
    key, sk = random.split(key)
    X = design_with_condition(sk, 400, 30, 10.0 ** e)
    y = X @ theta_star
    rel = lambda t: float(jnp.linalg.norm(t - theta_star) / jnp.linalg.norm(theta_star))
    t_ne = jnp.linalg.solve(X.T @ X, X.T @ y)
    t_qr = ridge_qr(X, y, 0.0)
    U, S, Vh = jnp.linalg.svd(X, full_matrices=False)
    t_svd = Vh.T @ ((U.T @ y) / S)
    print(f"     1e{e} {rel(t_ne):12.3e} {rel(t_qr):12.3e} {rel(t_svd):12.3e}")

Logistic regression by Newton's method

The IRLS loop, written as the Hessian solve it actually is. Note two implementation points that matter more than they look. The Hessian is formed as (X * W[:, None]).T @ X rather than by building a diagonal matrix, which turns an \(O(n^2)\) allocation into an \(O(nd)\) one. And the linear system is solved by Cholesky rather than a general solve, which is valid because the regularized Hessian is positive definite and is twice as fast. The JAX version wraps the whole iteration in lax.scan so the entire fit is one compiled kernel launch rather than one per iteration. PyTorch converges in 7 iterations with gradient norms 1.2e3, 2.9e2, 6.6e1, 6.3, 7.3e-2, 1.0e-5, 2.0e-13, while JAX, at eight fixed iterations, ends at 9.2e-14. The doubling of correct digits per step is visible in both.

def logistic_newton(X, y, l2=1.0, tol=1e-12, max_iter=50):
    """X: (n, d), y: (n,) in {0,1}. Returns theta (d,) and the gradient-norm trace."""
    n, d = X.shape
    theta = torch.zeros(d, dtype=X.dtype, device=X.device)
    I = torch.eye(d, dtype=X.dtype, device=X.device)
    trace = []
    for _ in range(max_iter):
        h = torch.sigmoid(X @ theta)                 # (n,)
        grad = X.T @ (h - y) + l2 * theta            # (d,)
        W = h * (1 - h)                              # (n,) Bernoulli variance
        H = (X * W.unsqueeze(1)).T @ X + l2 * I      # (d, d)
        step = torch.cholesky_solve(grad.unsqueeze(1), torch.linalg.cholesky(H)).squeeze(1)
        theta = theta - step
        trace.append(grad.norm().item())
        if step.norm() < tol:
            break
    return theta, trace

g = torch.Generator(device=dev).manual_seed(11)
X = torch.randn(4000, 12, generator=g, device=dev)
b = torch.randn(12, generator=g, device=dev) * 0.8
y = (torch.rand(4000, generator=g, device=dev) < torch.sigmoid(X @ b)).double()
theta, tr = logistic_newton(X, y)
print("iterations:", len(tr))
print("gradient norms:", [f"{v:.3e}" for v in tr])

def logistic_newton(X, y, l2=1.0, iters=8):
    """X: (n, d), y: (n,) in {0,1}. lax.scan keeps the whole loop inside one jit."""
    d = X.shape[1]
    I = jnp.eye(d)

    def step(theta, _):
        h = jax.nn.sigmoid(X @ theta)                       # (n,)
        grad = X.T @ (h - y) + l2 * theta                   # (d,)
        W = h * (1.0 - h)                                   # (n,)
        H = (X * W[:, None]).T @ X + l2 * I                 # (d, d)
        L = jnp.linalg.cholesky(H)
        return theta - jax.scipy.linalg.cho_solve((L, True), grad), jnp.linalg.norm(grad)

    return lax.scan(step, jnp.zeros(d), None, length=iters)

key = random.PRNGKey(11)
k1, k2, k3 = random.split(key, 3)
X = random.normal(k1, (4000, 12))
b = random.normal(k2, (12,)) * 0.8
y = (random.uniform(k3, (4000,)) < jax.nn.sigmoid(X @ b)).astype(jnp.float64)
theta, trace = jit(logistic_newton)(X, y)
print("gradient norms:", [f"{v:.3e}" for v in trace.tolist()])

One GLM, three links

This is the unification made concrete, a single fitting routine, with the family entering only as the pair \((a', a'')\), the mean function and the variance function. Swapping identity/ones for sigmoid/mu(1-mu) for exp/mu turns least squares into logistic regression into Poisson regression, and nothing else in the code changes. The Gaussian case converges in one step because the quadratic model of a quadratic is exact, and it agrees with the direct least-squares solution to 2.8e-11 in PyTorch and 1.7e-11 in JAX. Verified against independent implementations in the NumPy run, the Bernoulli fit matches scikit-learn's LogisticRegression to 1.4e-6 and the Poisson fit matches PoissonRegressor to 1.0e-8, both in 7 iterations. The step damping is the globalization discussed in Problem 4. Without it, the exponential link overshoots on the first step and can overflow.

LINKS = {
    "gaussian":  (lambda eta: eta,                 lambda mu: torch.ones_like(mu)),
    "bernoulli": (torch.sigmoid,                   lambda mu: mu * (1 - mu)),
    "poisson":   (torch.exp,                       lambda mu: mu),
}

def glm_irls(X, y, family, l2=1e-8, max_iter=50, tol=1e-12):
    """Canonical-link GLM by IRLS. The only family-specific objects are the
    mean function a'(eta) and the variance function a''(eta) = V(mu)."""
    mean_fn, var_fn = LINKS[family]
    n, d = X.shape
    theta = torch.zeros(d, dtype=X.dtype, device=X.device)
    I = torch.eye(d, dtype=X.dtype, device=X.device)
    for _ in range(max_iter):
        eta = X @ theta
        mu = mean_fn(eta)                            # (n,) = a'(eta)
        grad = X.T @ (mu - y) + l2 * theta           # the universal GLM gradient
        W = var_fn(mu)                               # (n,) = a''(eta)
        H = (X * W.unsqueeze(1)).T @ X + l2 * I
        step = torch.linalg.solve(H, grad.unsqueeze(1)).squeeze(1)
        theta = theta - step * min(1.0, 5.0 / max(1.0, step.norm().item()))
        if step.norm() < tol:
            break
    return theta

g = torch.Generator(device=dev).manual_seed(21)
X = torch.randn(800, 6, generator=g, device=dev)
b = torch.randn(6, generator=g, device=dev)
yg = X @ b + 0.5 * torch.randn(800, generator=g, device=dev)
yb = (torch.rand(800, generator=g, device=dev) < torch.sigmoid(X @ b)).double()
bp = torch.randn(6, generator=g, device=dev) * 0.3
yp = torch.poisson(torch.exp(X @ bp), generator=g)
Uo, So, Vho = torch.linalg.svd(X, full_matrices=False)
ols = Vho.T @ ((Uo.T @ yg) / So)
print("gaussian vs lstsq:", (glm_irls(X, yg, "gaussian") - ols).abs().max().item())
print("bernoulli coefs :", glm_irls(X, yb, "bernoulli")[:3].tolist())
print("poisson coefs   :", glm_irls(X, yp, "poisson")[:3].tolist(), "true", bp[:3].tolist())

# Each family is exactly two functions: the mean a'(eta) and the variance a''(eta).
LINKS = {
    "gaussian":  (lambda eta: eta,       lambda mu: jnp.ones_like(mu)),
    "bernoulli": (jax.nn.sigmoid,        lambda mu: mu * (1.0 - mu)),
    "poisson":   (jnp.exp,               lambda mu: mu),
}

def glm_irls(X, y, family, l2=1e-8, iters=12):
    mean_fn, var_fn = LINKS[family]
    d = X.shape[1]
    I = jnp.eye(d)

    def step(theta, _):
        mu = mean_fn(X @ theta)                             # a'(X theta)
        grad = X.T @ (mu - y) + l2 * theta                  # the one GLM gradient
        H = (X * var_fn(mu)[:, None]).T @ X + l2 * I        # a''-weighted Gram
        delta = jnp.linalg.solve(H, grad)
        # damp the first steps so exp-link fits cannot overshoot into overflow
        scale = jnp.minimum(1.0, 5.0 / jnp.maximum(1.0, jnp.linalg.norm(delta)))
        return theta - scale * delta, jnp.linalg.norm(grad)

    return lax.scan(step, jnp.zeros(d), None, length=iters)

key = random.PRNGKey(21)
k1, k2, k3, k4, k5 = random.split(key, 5)
X = random.normal(k1, (800, 6))
b = random.normal(k2, (6,))
yg = X @ b + 0.5 * random.normal(k3, (800,))
yb = (random.uniform(k4, (800,)) < jax.nn.sigmoid(X @ b)).astype(jnp.float64)
bp = random.normal(k5, (6,)) * 0.3
yp = random.poisson(random.PRNGKey(99), jnp.exp(X @ bp)).astype(jnp.float64)
U, S, Vh = jnp.linalg.svd(X, full_matrices=False)
ols = Vh.T @ ((U.T @ yg) / S)
tg, _ = glm_irls(X, yg, "gaussian")
tb, _ = glm_irls(X, yb, "bernoulli")
tp, _ = glm_irls(X, yp, "poisson")
print("gaussian vs OLS  :", float(jnp.abs(tg - ols).max()))
print("bernoulli coefs  :", [round(float(v), 5) for v in tb[:3]])
print("poisson coefs    :", [round(float(v), 5) for v in tp[:3]],
      "true", [round(float(v), 5) for v in bp[:3]])

Gaussian mixtures by EM, with log-sum-exp

The E step is where implementations go wrong, so it is worth reading closely. Densities are never exponentiated before normalizing. Log joint densities are formed via the Cholesky factor of each covariance (which gives both the Mahalanobis distance by triangular solve and the log determinant as twice the sum of the log diagonal, at no extra cost), then logsumexp normalizes in log space and the responsibilities come out of a single exp of a quantity that is guaranteed to be at most zero. The covariance update adds \(\epsilon I\) to prevent a component from collapsing onto a single point, which is a genuine singularity of the likelihood and not a numerical artifact. Both versions verify that the log-likelihood is monotone at every iteration, and the NumPy run additionally matches scikit-learn's GaussianMixture from the same initialization to 4.4e-16.

def gmm_em(X, K, iters=300, reg=1e-6, tol=1e-12):
    """X: (n, D). Returns pi (K,), mu (K, D), Sigma (K, D, D), log-likelihood trace."""
    n, D = X.shape
    pi = torch.full((K,), 1.0 / K, dtype=X.dtype, device=X.device)
    mu = X[:K].clone()
    Sigma = torch.cov(X.T).expand(K, D, D).clone()
    lls = []
    for _ in range(iters):
        L = torch.linalg.cholesky(Sigma)                       # (K, D, D)
        diff = X.unsqueeze(0) - mu.unsqueeze(1)                # (K, n, D)
        sol = torch.linalg.solve_triangular(L, diff.transpose(1, 2), upper=False)
        maha = (sol ** 2).sum(1)                               # (K, n)
        logdet = 2 * torch.log(torch.diagonal(L, dim1=1, dim2=2)).sum(1)   # (K,)
        logp = (torch.log(pi).unsqueeze(1) - 0.5 * (maha + logdet.unsqueeze(1)
                + D * math.log(2 * math.pi))).T                # (n, K)
        ll = torch.logsumexp(logp, dim=1)                      # (n,) stable
        lls.append(ll.mean().item())
        R = torch.exp(logp - ll.unsqueeze(1))                  # (n, K) responsibilities
        Nk = R.sum(0)                                          # (K,)
        pi = Nk / n
        mu = (R.T @ X) / Nk.unsqueeze(1)
        cen = X.unsqueeze(0) - mu.unsqueeze(1)                 # (K, n, D)
        Sigma = torch.einsum("kn,knd,kne->kde", R.T, cen, cen) / Nk.view(K, 1, 1)
        Sigma = Sigma + reg * torch.eye(D, dtype=X.dtype, device=X.device)
        if len(lls) > 1 and abs(lls[-1] - lls[-2]) < tol:
            break
    return pi, mu, Sigma, lls

g = torch.Generator(device=dev).manual_seed(5)
mus = torch.tensor([[0.0, 0.0], [4.0, 1.0], [1.0, 5.0]], device=dev)
z = torch.randint(0, 3, (1500,), generator=g, device=dev)
X = mus[z] + torch.randn(1500, 2, generator=g, device=dev)
pi, mu, Sig, lls = gmm_em(X, 3)
print("EM iterations:", len(lls), " monotone:",
      all(lls[i + 1] >= lls[i] - 1e-12 for i in range(len(lls) - 1)))
print("log-likelihood:", f"{lls[0]:.6f} -> {lls[-1]:.6f}")
print("mixing weights:", sorted(round(v, 4) for v in pi.tolist()))

def gmm_em(X, K, iters=120, reg=1e-6):
    n, D = X.shape

    def e_step(params):
        pi, mu, Sigma = params
        L = jnp.linalg.cholesky(Sigma)                              # (K, D, D)
        diff = X[None, :, :] - mu[:, None, :]                       # (K, n, D)
        sol = jax.vmap(lambda l, dd: jax.scipy.linalg.solve_triangular(
            l, dd.T, lower=True))(L, diff)                          # (K, D, n)
        maha = jnp.sum(sol ** 2, axis=1)                            # (K, n)
        logdet = 2.0 * jnp.sum(jnp.log(jnp.diagonal(L, axis1=1, axis2=2)), axis=1)
        logp = (jnp.log(pi)[:, None]
                - 0.5 * (maha + logdet[:, None] + D * jnp.log(2 * jnp.pi))).T   # (n, K)
        ll = logsumexp(logp, axis=1)                                # (n,)
        return jnp.exp(logp - ll[:, None]), ll.mean()

    def step(params, _):
        R, ll = e_step(params)                                      # (n, K)
        Nk = R.sum(0)                                               # (K,)
        pi = Nk / n
        mu = (R.T @ X) / Nk[:, None]
        cen = X[None, :, :] - mu[:, None, :]                        # (K, n, D)
        Sigma = jnp.einsum("nk,knd,kne->kde", R, cen, cen) / Nk[:, None, None]
        return (pi, mu, Sigma + reg * jnp.eye(D)), ll

    init = (jnp.full((K,), 1.0 / K), X[:K], jnp.tile(jnp.cov(X.T), (K, 1, 1)))
    (pi, mu, Sigma), lls = lax.scan(step, init, None, length=iters)
    return pi, mu, Sigma, lls

key = random.PRNGKey(5)
k1, k2 = random.split(key)
mus = jnp.array([[0.0, 0.0], [4.0, 1.0], [1.0, 5.0]])
z = random.randint(k1, (1500,), 0, 3)
X = mus[z] + random.normal(k2, (1500, 2))
pi, mu, Sigma, lls = jit(gmm_em, static_argnums=(1,))(X, 3)
d = jnp.diff(lls)
print("monotone:", bool(jnp.all(d >= -1e-12)),
      f" log-likelihood {float(lls[0]):.6f} -> {float(lls[-1]):.6f}")
print("mixing weights:", sorted(round(float(v), 4) for v in pi))

PCA by SVD, and probabilistic PCA

PCA is computed from the SVD of the centered data, never from an explicit covariance matrix, for the same squared-condition-number reason that governs least squares. The block also checks the equivalence of the two derivations numerically. Retained variance plus reconstruction error minus total variance is 2.8e-14 in PyTorch and 1.4e-14 in JAX, which is the Pythagorean identity holding to machine precision. Probabilistic PCA is three extra lines on top of the same factorization, since its MLE is closed-form, and it recovers the injected noise variance of 0.16 as 0.1606 (PyTorch) and 0.1581 (JAX).

def pca_svd(X, k):
    """X: (n, d). Returns components (k, d), explained variance (k,), scores (n, k)."""
    mean = X.mean(0)
    Xc = X - mean
    U, S, Vh = torch.linalg.svd(Xc, full_matrices=False)       # Vh: (d, d)
    var = S ** 2 / (X.shape[0] - 1)
    return Vh[:k], var[:k], Xc @ Vh[:k].T

def ppca(X, q):
    """Tipping-Bishop closed-form MLE: sigma^2 is the mean discarded eigenvalue."""
    n, d = X.shape
    Xc = X - X.mean(0)
    U, S, Vh = torch.linalg.svd(Xc, full_matrices=False)
    lam = S ** 2 / n
    sigma2 = lam[q:].mean()
    W = Vh[:q].T @ torch.diag(torch.sqrt(torch.clamp(lam[:q] - sigma2, min=0.0)))
    return W, sigma2

g = torch.Generator(device=dev).manual_seed(9)
Wt = torch.randn(10, 3, generator=g, device=dev)
Z = torch.randn(600, 3, generator=g, device=dev)
X = Z @ Wt.T + 0.4 * torch.randn(600, 10, generator=g, device=dev) + 5.0
comp, var, scores = pca_svd(X, 3)
total = ((X - X.mean(0)) ** 2).sum() / (X.shape[0] - 1)
print("explained variance ratio:", [round(v, 5) for v in (var / total).tolist()])
Xc = X - X.mean(0)
recon = ((Xc - (Xc @ comp.T) @ comp) ** 2).sum(1).mean()
kept = ((Xc @ comp.T) ** 2).sum(1).mean()
tot = (Xc ** 2).sum(1).mean()
print(f"reconstruction + kept - total = {(recon + kept - tot).item():.3e}  (must be 0)")
W, s2 = ppca(X, 3)
print("PPCA sigma^2:", round(s2.item(), 6), "(true 0.16)")

@functools.partial(jit, static_argnums=(1,))
def pca_svd(X, k):
    Xc = X - X.mean(0)
    U, S, Vh = jnp.linalg.svd(Xc, full_matrices=False)
    return Vh[:k], S[:k] ** 2 / (X.shape[0] - 1), Xc @ Vh[:k].T

@functools.partial(jit, static_argnums=(1,))
def ppca(X, q):
    n = X.shape[0]
    Xc = X - X.mean(0)
    _, S, Vh = jnp.linalg.svd(Xc, full_matrices=False)
    lam = S ** 2 / n
    sigma2 = lam[q:].mean()
    return Vh[:q].T @ jnp.diag(jnp.sqrt(jnp.clip(lam[:q] - sigma2, 0.0))), sigma2

key = random.PRNGKey(9)
k1, k2, k3 = random.split(key, 3)
Wt = random.normal(k1, (10, 3))
Z = random.normal(k2, (600, 3))
X = Z @ Wt.T + 0.4 * random.normal(k3, (600, 10)) + 5.0
comp, var, scores = pca_svd(X, 3)
Xc = X - X.mean(0)
total = (Xc ** 2).sum(1).mean()
kept = ((Xc @ comp.T) ** 2).sum(1).mean()
recon = ((Xc - (Xc @ comp.T) @ comp) ** 2).sum(1).mean()
print("explained variance ratio:",
      [round(float(v), 5) for v in var / ((Xc ** 2).sum() / (X.shape[0] - 1))])
print(f"reconstruction + kept - total = {float(recon + kept - total):.3e}  (must be 0)")
W, s2 = ppca(X, 3)
print("PPCA sigma^2:", round(float(s2), 6), "(true 0.16)")

Gradient boosting with second-order leaf weights

The XGBoost objective implemented directly, with leaf values \(-G_j/(H_j+\lambda)\), split gain from the difference of \(G^2/(H+\lambda)\) scores, and \(\gamma\) as a minimum gain threshold. The split scan is vectorized by sorting each feature once and taking cumulative sums of the gradients and Hessians, so all \(m-1\) candidate thresholds for a feature are scored in one pass rather than in a loop. That single change is the difference between a usable implementation and a toy. The recursion itself stays in Python, because tree structure is data-dependent control flow that neither framework can trace. Squared loss makes \(h_i = 1\), so this reduces to shrunken residual fitting, but changing two lines to \(g_i = p_i - y_i\), \(h_i = p_i(1-p_i)\) gives logistic boosting with no other edits. Test MSE after 200 rounds is 2.1814 (PyTorch) and 2.1777 (JAX) against scikit-learn's first-order implementation at 2.2829 with identical hyperparameters.

def fit_tree(X, gvec, hvec, depth, lam, gamma, min_child=1.0):
    """Exact greedy tree on second-order statistics. X: (n, d)."""
    n, d = X.shape

    def score(G, H):
        return G * G / (H + lam)

    def build(idx, dep):
        G, H = gvec[idx].sum(), hvec[idx].sum()
        leaf = {"leaf": True, "w": (-G / (H + lam)).item()}
        if dep >= depth or idx.numel() < 2:
            return leaf
        base, best = score(G, H), (0.0, None, None)
        for j in range(d):
            order = idx[torch.argsort(X[idx, j])]
            gs, hs = gvec[order].cumsum(0), hvec[order].cumsum(0)
            vals = X[order, j]
            GL, HL = gs[:-1], hs[:-1]
            GR, HR = G - GL, H - HL
            gain = 0.5 * (score(GL, HL) + score(GR, HR) - base) - gamma
            ok = (vals[:-1] != vals[1:]) & (HL >= min_child) & (HR >= min_child)
            gain = torch.where(ok, gain, torch.full_like(gain, -1.0))
            k = int(gain.argmax())
            if gain[k].item() > best[0]:
                best = (gain[k].item(), j, 0.5 * (vals[k] + vals[k + 1]).item())
        if best[1] is None:
            return leaf
        _, j, thr = best
        m = X[idx, j] <= thr
        return {"leaf": False, "j": j, "thr": thr,
                "L": build(idx[m], dep + 1), "R": build(idx[~m], dep + 1)}

    return build(torch.arange(n, device=X.device), 0)

def tree_predict(tree, X):
    out = torch.zeros(X.shape[0], dtype=X.dtype, device=X.device)

    def walk(node, idx):
        if node["leaf"]:
            out[idx] = node["w"]
            return
        m = X[idx, node["j"]] <= node["thr"]
        walk(node["L"], idx[m])
        walk(node["R"], idx[~m])

    walk(tree, torch.arange(X.shape[0], device=X.device))
    return out

def boost(Xtr, ytr, Xte, rounds=200, eta=0.1, depth=3, lam=1.0, gamma=0.0):
    F = torch.full_like(ytr, ytr.mean().item())
    Fe = torch.full((Xte.shape[0],), ytr.mean().item(), dtype=ytr.dtype, device=ytr.device)
    for _ in range(rounds):
        gvec = F - ytr                       # squared loss: g = dL/dF
        hvec = torch.ones_like(gvec)         #               h = d2L/dF2 = 1
        tree = fit_tree(Xtr, gvec, hvec, depth, lam, gamma)
        F = F + eta * tree_predict(tree, Xtr)
        Fe = Fe + eta * tree_predict(tree, Xte)
    return F, Fe

from sklearn.datasets import make_friedman1
import numpy as np
Xtr_n, ytr_n = make_friedman1(n_samples=800, noise=1.0, random_state=0)
Xte_n, yte_n = make_friedman1(n_samples=2000, noise=1.0, random_state=1)
Xtr = torch.tensor(Xtr_n, device=dev); ytr = torch.tensor(ytr_n, device=dev)
Xte = torch.tensor(Xte_n, device=dev); yte = torch.tensor(yte_n, device=dev)
F, Fe = boost(Xtr, ytr, Xte)
print("train MSE:", round(((F - ytr) ** 2).mean().item(), 4),
      " test MSE:", round(((Fe - yte) ** 2).mean().item(), 4))

# Trees are data-dependent control flow, so the tree structure is built in Python
# while every split scan is a vectorised jnp reduction.

@jit
def best_split(Xs, g, h, lam, gamma):
    """Xs: (m, d) node rows, g/h: (m,). Returns (gain, feature, threshold)."""
    order = jnp.argsort(Xs, axis=0)                      # (m, d)
    gs = jnp.cumsum(jnp.take_along_axis(g[:, None], order, 0), 0)   # (m, d)
    hs = jnp.cumsum(jnp.take_along_axis(h[:, None], order, 0), 0)
    vals = jnp.take_along_axis(Xs, order, 0)
    G, H = gs[-1], hs[-1]
    GL, HL = gs[:-1], hs[:-1]
    GR, HR = G - GL, H - HL
    score = lambda a, b: a * a / (b + lam)
    gain = 0.5 * (score(GL, HL) + score(GR, HR) - score(G, H)) - gamma
    gain = jnp.where(vals[:-1] != vals[1:], gain, -jnp.inf)
    flat = jnp.argmax(gain)
    i, j = jnp.unravel_index(flat, gain.shape)
    return gain[i, j], j, 0.5 * (vals[i, j] + vals[i + 1, j])

def fit_tree(X, g, h, depth, lam, gamma):
    def build(idx, dep):
        G, H = float(g[idx].sum()), float(h[idx].sum())
        leaf = {"leaf": True, "w": -G / (H + lam)}
        if dep >= depth or idx.size < 2:
            return leaf
        gain, j, thr = best_split(X[idx], g[idx], h[idx], lam, gamma)
        if not (float(gain) > 0.0):
            return leaf
        j, thr = int(j), float(thr)
        m = X[idx, j] <= thr
        return {"leaf": False, "j": j, "thr": thr,
                "L": build(idx[m], dep + 1), "R": build(idx[~m], dep + 1)}
    import numpy as np
    return build(np.arange(X.shape[0]), 0)

def tree_predict(tree, X):
    import numpy as np
    out = np.zeros(X.shape[0])

    def walk(node, idx):
        if node["leaf"]:
            out[idx] = node["w"]
            return
        m = np.asarray(X[idx, node["j"]] <= node["thr"])
        walk(node["L"], idx[m])
        walk(node["R"], idx[~m])

    walk(tree, np.arange(X.shape[0]))
    return jnp.asarray(out)

from sklearn.datasets import make_friedman1
import numpy as np
Xtr, ytr = make_friedman1(n_samples=800, noise=1.0, random_state=0)
Xte, yte = make_friedman1(n_samples=2000, noise=1.0, random_state=1)
Xtr, ytr = jnp.asarray(Xtr), jnp.asarray(ytr)
Xte, yte = jnp.asarray(Xte), jnp.asarray(yte)
F = jnp.full(ytr.shape, ytr.mean())
Fe = jnp.full(yte.shape, ytr.mean())
for _ in range(200):
    gvec = F - ytr                       # squared loss: first derivative
    hvec = jnp.ones_like(gvec)           #               second derivative
    tree = fit_tree(np.asarray(Xtr), np.asarray(gvec), np.asarray(hvec), 3, 1.0, 0.0)
    F = F + 0.1 * tree_predict(tree, np.asarray(Xtr))
    Fe = Fe + 0.1 * tree_predict(tree, np.asarray(Xte))
print("train MSE:", round(float(((F - ytr) ** 2).mean()), 4),
      " test MSE:", round(float(((Fe - yte) ** 2).mean()), 4))

How it is done in practice

The gap between the derivations above and a deployed system is mostly numerical linear algebra, memory movement, and the discipline of not leaking information between training and evaluation. A few things that separate working systems from notebook code.

Nobody solves the problem you derived. They solve a preconditioned version of it. Features get centered and scaled before any penalized fit, because a ridge or lasso penalty is not scale-invariant. Doubling a column halves the coefficient that the same penalty tolerates, so unscaled features mean an arbitrary, silent, per-feature regularization strength. The intercept is excluded from the penalty for the same reason. Scaling parameters must be estimated on the training fold only and applied to the validation fold, which is why every serious pipeline wraps the scaler and the model in a single object that is fitted as a unit. Scaling before splitting is the single most common source of optimistic cross-validation results.

The solver is chosen by shape, not by taste. For \(n \gg d\), form the \(d \times d\) Gram matrix and Cholesky-solve. The Gram formation is a single large matmul that runs at near-peak throughput. Measured on an H100 80GB HBM3 in float64, the full ridge fit (Gram, regularize, Cholesky, triangular solves) takes 8.63 ms at \(n = 200{,}000, d = 1024\), 132 ms at \(n = 500{,}000, d = 2048\), and 1.20 s at \(n = 10^6, d = 4096\), sustaining 48.6, 31.8, and 28.0 effective TFLOP/s respectively. The float32 versions of the same three are 8.82 ms, 81.1 ms, and 695 ms. That fp64 and fp32 are close at \(d = 1024\) is not a bug. The H100's fp64 tensor cores and its non-TF32 fp32 path have similar peak rates, so at these shapes fp64 costs almost nothing and there is no reason to give up the precision. For \(d \gg n\), work in the dual and solve the \(n \times n\) system instead. For both large, use an iterative method (LSQR or conjugate gradients on the normal equations) that only needs matrix-vector products, or subsample.

Second-order methods are used exactly where they fit. A logistic regression with \(n = 10^6\) and \(d = 256\) fits by Newton's method in 7 iterations and 0.101 s on an H100 in float32, 0.110 s in float64, recovering the generating coefficients at correlation 0.9999, and each iteration is dominated by one \((X\T W)X\) matmul. That is faster than any first-order method will manage, and it is why statsmodels and every GLM package use IRLS rather than SGD. The moment \(d\) reaches \(10^5\) the \(d \times d\) Hessian no longer fits and the calculus reverses. This is the whole reason deep learning uses first-order methods, and the various quasi-Newton and Kronecker-factored approximations (L-BFGS, K-FAC, Shampoo) exist to recover some of the curvature information without materializing the matrix.

Boosted-tree training is memory-bound, not compute-bound. The inner loop is a histogram accumulation over gradients and Hessians, so the engineering that matters is data layout (column-major, pre-binned to uint8), histogram subtraction to halve the work per level, and cache-friendly access patterns. This is why LightGBM's speedups come from algorithmic bookkeeping rather than floating-point throughput, and why GPU implementations of GBDT are a smaller win than one might expect, since the operation is scattered accumulation, not dense matmul.

Evaluation discipline is where most projects actually fail. Any tuning done on a validation set makes that set optimistic, so a final untouched test set is not optional. Time-series data needs forward-chaining splits, and grouped data (multiple rows per customer, patient, or document) needs group-aware splits, or the model memorizes the group and the score is fiction. Target encoding of categorical variables must be computed out-of-fold or it leaks the label directly, which is exactly what CatBoost's ordered statistics eliminate. And outputs should be monitored for calibration drift as well as accuracy drift, because an input shift can break a threshold long before it breaks a ranking.

The current research frontier

Why interpolation works. The double-descent observation of Belkin, Hsu, Ma, and Mandal has matured into a fairly precise theory. Bartlett, Long, Lugosi, and Tsigler (2020) characterized benign overfitting for minimum-norm linear regression. The interpolant generalizes when the covariance spectrum has a small number of large eigenvalues carrying signal and a long tail of small ones that absorb the noise harmlessly, and they give matching conditions on the effective ranks. Hastie, Montanari, Rosset, and Tibshirani analysed the same regime under proportional asymptotics (\(d/n \to \gamma\)) and showed that optimally tuned ridge is monotone, which is the technical version of the practical claim above that double descent is largely a consequence of not regularizing. Work at ETH Zurich by Bühlmann's and Yang's groups on high-dimensional inference and the interpolation regime continues in this direction, as does van de Geer's line on lasso theory that supplies the sharp oracle inequalities.

Kernels as the tractable model of neural networks. The neural tangent kernel of Jacot, Gabriel, and Hongler at EPFL made infinitely wide networks exactly equivalent to kernel regression with a specific kernel, which turned a class of questions about deep learning into questions about the classical theory on this page. The follow-up literature is largely about the gap. Finite networks learn features and the NTK does not, and the mean-field and \(\mu\)P parameterizations of Yang and coauthors describe regimes where feature learning survives the width limit. Random features remain the practical residue. The Performer attention approximation and a family of scalable Gaussian-process methods are Rahimi-Recht applied to modern architectures.

Distribution-free uncertainty. Conformal prediction has moved from a niche technique to a standard tool, largely through work by Vovk at Royal Holloway, Ramdas at CMU, and Barber, Candès, and coauthors. The appeal is a theorem with almost no assumptions. From a held-out calibration set, take the \(\lceil(n+1)(1-\alpha)\rceil\)-th smallest conformity score and use it as a threshold, and the resulting prediction set covers the truth with probability at least \(1-\alpha\) for any model and any distribution, requiring only exchangeability. It is the correct answer to "my model is uncalibrated and I need a guarantee", and the active questions are conditional coverage (marginal coverage can hide systematic failure on subgroups) and validity under distribution shift, where exchangeability breaks.

Tabular deep learning, still trying. After Grinsztajn and coauthors and Shwartz-Ziv and Armon, the response has been architectural (FT-Transformer, TabNet, NODE, SAINT) and, more recently, in-context. TabPFN, from Hutter's group at Freiburg, pretrains a transformer on millions of synthetic tabular tasks so that inference on a new small dataset is a single forward pass with no gradient steps at all, and it is genuinely competitive with tuned gradient boosting on small problems. Whether this scales past a few thousand rows is the open question, and it is the most interesting thing to have happened to tabular learning in a decade.

Causality and distribution shift. A line of work running through Bühlmann at ETH Zurich (invariant causal prediction), Peters and Schölkopf at Max Planck, and Arjovsky and coauthors (invariant risk minimization) attacks the observation that empirical risk minimization has no defence against a shift in the environment. It will happily use a spurious correlation that holds in training and breaks in deployment. The proposal is to search for predictors whose conditional distribution is stable across environments, on the grounds that stability is a signature of causal structure. Whether the objectives proposed so far actually deliver this outside carefully constructed settings is contested, and the negative results are as informative as the positive ones.

Optimization theory catching up with practice. Work at the Technion on first-order, proximal, and online methods, and at Weizmann by Shamir on the limits of stochastic optimization, has produced sharp lower bounds explaining where adaptivity genuinely helps and where speedups are impossible. The practical residue is that the accelerated proximal machinery derived for the lasso now underpins most structured-sparsity and low-rank estimation.

Open source to read

Repositories worth reading rather than merely importing, with the file to open first.

scikit-learn/scikit-learn is the reference implementation of essentially everything on this page, and unusually readable for a library of its size. Open sklearn/linear_model/_ridge.py and follow the solver dispatch. _solve_cholesky, _solve_svd, and _solve_lsqr are the three routes from the numerical section, each chosen by problem shape, with the \(n < d\) case explicitly switching to the dual form. Then sklearn/mixture/_gaussian_mixture.py for EM done carefully with Cholesky precisions, and sklearn/calibration.py for Platt scaling and isotonic regression side by side.

dmlc/xgboost is the second-order derivation in production C++. Open src/tree/param.h, which contains CalcGain and CalcWeight. Those two functions are \(G^2/(H+\lambda)\) and \(-G/(H+\lambda)\), the closed forms derived above, with the weight-clipping and \(\alpha\)-regularization branches that the paper omits. Then src/tree/updater_quantile_hist.cc for the histogram-based split finding that makes it fast.

microsoft/LightGBM is the same objective with better data structures. Open src/treelearner/serial_tree_learner.cpp and look for ConstructHistograms. The histogram-subtraction trick, where a child's histogram is computed as the parent's minus its sibling's, is the single algorithmic idea behind most of the speed advantage. src/treelearner/feature_histogram.hpp has the leaf-output and gain formulas.

catboost/catboost is worth reading for the categorical handling and the target-leakage fix rather than the boosting. Open catboost/libs/algo/greedy_tensor_search.cpp for the split search over feature combinations, and follow how ordered target statistics are computed from a random permutation prefix so that an example never contributes to its own encoding.

statsmodels/statsmodels is the classical statistics side, with the inference that machine-learning libraries omit. Open statsmodels/genmod/generalized_linear_model.py, where the fit method is IRLS exactly as derived here, and the family classes in statsmodels/genmod/families/ are the \((a', a'')\) pairs from the GLM table written out one file at a time, including the non-canonical links and the deviance definitions.

google/jax is worth reading for how the linear algebra is actually plumbed. Open jax/_src/lax/linalg.py and follow cholesky, qr, and svd through their abstract evaluation rules, batching rules, and JVP rules. The derivative of a Cholesky factorization is not obvious and the code is the clearest available derivation of it.

pyro-ppl/pyro is the EM-to-variational-inference bridge in executable form. Open pyro/infer/svi.py for the optimization loop stripped to essentials, then pyro/infer/trace_elbo.py to see the exact ELBO from the EM section being estimated from execution traces, with the score-function and reparameterized gradient estimators as separate paths.

scikit-learn-contrib/imbalanced-learn collects every resampling method in one place, which is the fastest way to see how modest most of them are. Open imblearn/over_sampling/_smote/base.py. SMOTE's core is roughly fifteen lines of nearest-neighbour interpolation, and reading it makes clear both why it sometimes helps and why it damages calibration.

Common misconceptions

"The normal equations are how least squares is computed." They are how it is derived. Forming \(X\T X\) squares the condition number, and the measurement above shows the normal-equations solution reaching 100% relative error at \(\kappa(X) = 10^8\) where QR is still accurate to nine digits, and failing already at \(\kappa(X) = 10^4\) in single precision. Every production solver factors \(X\) directly or regularizes first.

"Regularization works by preventing the model from fitting the training data." Regularization works by trading bias for variance, and the trade can be computed exactly. In Problem 2, the optimal ridge penalty is \(\lambda^\ast = \sigma^2/(\theta^\ast)^2\), a ratio of noise power to signal power, and the biased estimator beats the unbiased one by 11%. The mechanism is shrinkage of the low-signal directions of the design, visible in the SVD as a factor \(d_i^2/(d_i^2+\lambda)\) that is 0.999 on the strongest direction and \(2\times10^{-10}\) on the weakest. Nothing about it is a restraint on effort.

"L1 gives sparsity because it penalizes small coefficients more." It does not penalize small coefficients more, since \(|t|\) is smaller than \(t^2\) for \(|t| < 1\). Sparsity comes from non-differentiability at zero. The subdifferential of \(|t|\) at 0 is the whole interval \([-1,1]\), so zero stays optimal for a whole range of data gradients, which is exactly the soft-thresholding condition \(|z_j| \le \lambda\). A smooth penalty has a single-valued gradient at zero and can never produce an exact zero.

"Logistic regression and linear regression are unrelated models." They are the same model with different log-partition functions. Both are canonical-link GLMs, both have gradient \(X\T(\hat y - y)\), both have Hessian \(X\T W X\) with \(W\) the conditional variances, and both are fit by the identical IRLS loop. The Gaussian case just happens to converge in one step because its \(W\) is constant. One 12-line routine in the implementation section fits both, plus Poisson, by swapping two functions.

"Naive Bayes works because features really are conditionally independent." They almost never are, and the estimated probabilities are correspondingly terrible, usually pinned near 0 or 1 because correlated evidence gets counted several times. It works when it works because the Bayes rule needs only the argmax of the posteriors, and the ordering survives distortions that destroy the magnitudes. Naive Bayes is a decent classifier and a bad probability estimator, and using its outputs as probabilities is a real error.

"Support vectors are the points on the margin." One direction only. Every point with \(\alpha_i > 0\) lies on the margin by complementary slackness, but a point can lie exactly on the margin and still have \(\alpha_i = 0\), contributing nothing. Problem 6 exhibits exactly this. \(x_2 = (2,2)\) has functional margin exactly 1 and dual variable exactly 0, and deleting it changes the solution not at all.

"More model capacity always eventually hurts." The classical U-curve is correct in the underparameterized regime and is not the whole picture. Past the interpolation threshold, test error can fall again and end below the classical minimum. In the measured sweep, error peaks at 16.2 near \(p = n\) and descends to 0.101 at \(p = 1000\), six times better than the best underparameterized model. The complementary correction is equally important. A tiny ridge penalty cuts that peak by a factor of three to six and makes the curve nearly monotone, so double descent is mostly a symptom of leaving an ill-conditioned problem unregularized.

"A model with 0.95 AUC is a good model." On imbalanced data, quite possibly not. The measured example has classifier B beating classifier A on ROC-AUC (0.907 against 0.882) while having less than half its average precision (0.215 against 0.417) and finding 62 real positives in its top 100 against A's 100. ROC-AUC's false positive rate is diluted by a large negative class, so it is nearly blind to differences at the top of the ranking, which is the only part anyone acts on when positives are rare.

"Cross-validation gives an unbiased estimate of test error." It gives a slightly pessimistic estimate of the error of a model trained on \(n(k-1)/k\) points, not of the model you ship, which was trained on all \(n\). As measured, 2-fold CV overstated the true error by 50% relative, 10-fold by 6%, and LOO by 2.5%. And any hyperparameter chosen using the CV score makes that score optimistic in turn, which is why nested cross-validation or a held-out test set exists.

Self-check

References

1. Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning, 2nd edition. Springer, 2009. doi:10.1007/978-0-387-84858-7. Chapters 3, 4, 7, 10, and 14 map onto the sections here.
2. Bishop, C. M. Pattern Recognition and Machine Learning. Springer, 2006. The cleanest treatment of the exponential family, EM from the lower bound, and probabilistic PCA, written at Microsoft Research Cambridge.
3. Bishop, C. M. and Bishop, H. Deep Learning, Foundations and Concepts. Springer, 2024. bishopbook.com. The GLM material recast as the output layers of networks.
4. Murphy, K. P. Probabilistic Machine Learning, An Introduction (2022) and Advanced Topics (2023). MIT Press. probml.github.io/pml-book. The best single source for the connections between the classical and modern literature.
5. Wasserman, L. All of Statistics, A Concise Course in Statistical Inference. Springer, 2004. doi:10.1007/978-0-387-21736-9. From Carnegie Mellon, the fastest route to the inference background assumed here.
6. Shalev-Shwartz, S. and Ben-David, S. Understanding Machine Learning, From Theory to Algorithms. Cambridge University Press, 2014. doi:10.1017/CBO9781107298019. The estimation-error half of the decomposition, and the Pegasos SVM solver used in the measurements.
7. Boyd, S. and Vandenberghe, L. Convex Optimization. Cambridge University Press, 2004. cvxbook. Chapter 5 supplies the Lagrangian duality and KKT conditions used in the SVM derivation.
8. McCullagh, P. and Nelder, J. A. Generalized Linear Models, 2nd edition. Chapman and Hall, 1989. doi:10.1201/9780203753736. Deviance, quasi-likelihood, and the non-canonical links this page only mentions.
9. Nelder, J. A. and Wedderburn, R. W. M. Generalized linear models. Journal of the Royal Statistical Society, Series A, 135(3):370-384, 1972. doi:10.2307/2344614. The paper that unified regression, logit, probit, and log-linear models, and introduced IRLS as their common fitting algorithm.
10. Cortes, C. and Vapnik, V. Support-vector networks. Machine Learning, 20(3):273-297, 1995. doi:10.1007/BF00994018. The soft-margin SVM as it is still used, from AT&T Bell Labs.
11. Ng, A. Y. and Jordan, M. I. On discriminative vs. generative classifiers, a comparison of logistic regression and naive Bayes. NeurIPS, 2001. NeurIPS proceedings. Work done at Berkeley, with the sample-complexity argument behind the crossing learning curves measured above.
12. Dempster, A. P., Laird, N. M., and Rubin, D. B. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39(1):1-38, 1977. doi:10.1111/j.2517-6161.1977.tb01600.x. From Harvard, the paper that named EM and proved its monotonicity in general.
13. Tibshirani, R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 58(1):267-288, 1996. doi:10.1111/j.2517-6161.1996.tb02080.x.
14. Zou, H. and Hastie, T. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society, Series B, 67(2):301-320, 2005. doi:10.1111/j.1467-9868.2005.00503.x.
15. Beck, A. and Teboulle, M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183-202, 2009. doi:10.1137/080716542. From the Technion and Tel Aviv, the accelerated proximal method whose proximal operator is exactly the soft-thresholding formula derived here.
16. Breiman, L. Bagging predictors. Machine Learning, 24(2):123-140, 1996. doi:10.1007/BF00058655, and Random forests. Machine Learning, 45(1):5-32, 2001. doi:10.1023/A:1010933404324. Both from Berkeley, and the second contains the correlation argument measured above.
17. Freund, Y. and Schapire, R. E. A decision-theoretic generalization of on-line learning and an application to boosting. Journal of Computer and System Sciences, 55(1):119-139, 1997. doi:10.1006/jcss.1997.1504, and Schapire, Freund, Bartlett, and Lee, Boosting the margin, Annals of Statistics 26(5):1651-1686, 1998, doi:10.1214/aos/1024691352, for the margin explanation of why boosting keeps improving after zero training error.
18. Friedman, J. H. Greedy function approximation, a gradient boosting machine. Annals of Statistics, 29(5):1189-1232, 2001. doi:10.1214/aos/1013203451. Boosting as gradient descent in function space.
19. Chen, T. and Guestrin, C. XGBoost, a scalable tree boosting system. KDD, 2016. arXiv:1603.02754. From the University of Washington, the second-order objective and closed-form leaf weights derived above.
20. Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., and Liu, T.-Y. LightGBM, a highly efficient gradient boosting decision tree. NeurIPS, 2017. NeurIPS proceedings. Microsoft Research, with histogram subtraction, GOSS, and exclusive feature bundling. See also Prokhorenkova et al., CatBoost, NeurIPS 2018, arXiv:1706.09516, for ordered boosting and target-leak-free categorical encoding.
21. Rahimi, A. and Recht, B. Random features for large-scale kernel machines. NeurIPS, 2007. NeurIPS proceedings. Done at Intel Research Berkeley, with the Bochner argument and the \(D^{-1/2}\) rate confirmed in the measurements.
22. Belkin, M., Hsu, D., Ma, S., and Mandal, S. Reconciling modern machine-learning practice and the classical bias-variance trade-off. PNAS, 116(32):15849-15854, 2019. doi:10.1073/pnas.1903070116, then Nakkiran, P., Kaplun, G., Bansal, Y., Yang, T., Barak, B., and Sutskever, I. Deep double descent. ICLR, 2020. arXiv:1912.02292, and Bartlett, P. L., Long, P. M., Lugosi, G., and Tsigler, A. Benign overfitting in linear regression. PNAS, 117(48):30063-30070, 2020. doi:10.1073/pnas.1907378117. The last, from Berkeley, gives the spectral conditions under which interpolation is harmless.
23. Grinsztajn, L., Oyallon, E., and Varoquaux, G. Why do tree-based models still outperform deep learning on typical tabular data? NeurIPS Datasets and Benchmarks, 2022. arXiv:2207.08815. From INRIA, the three-property diagnosis reproduced in the bake-off above.
24. Guo, C., Pleiss, G., Sun, Y., and Weinberger, K. Q. On calibration of modern neural networks. ICML, 2017. arXiv:1706.04599. From Cornell, temperature scaling and the ECE diagnostic. The classical antecedents are Platt (1999) for sigmoid scaling and Zadrozny and Elkan (KDD 2002) for isotonic regression.
25. Tipping, M. E. and Bishop, C. M. Probabilistic principal component analysis. Journal of the Royal Statistical Society, Series B, 61(3):611-622, 1999. doi:10.1111/1467-9868.00196, and Jacot, A., Gabriel, F., and Hongler, C. Neural tangent kernel, convergence and generalization in neural networks. NeurIPS, 2018. arXiv:1806.07572. The latter, from EPFL, is why the kernel material on this page became current again.