The classical U-curve and why it predicts the wrong thing
Fix a target f* and a noisy label model y = f*(x) + ε with E[ε] = 0, Var(ε) = σ^2. Train an estimator f̂ on a random dataset D of N points. Classical learning theory says model capacity is a dial: turn it up and training error falls monotonically, while test error first falls (less bias) then rises (more variance), giving the familiar U with a single optimum somewhere before you can fit the noise.
The prediction that breaks is what happens at the right end. The classical curve says the interpolation point — where f̂ hits every training label exactly — is a disaster, and beyond it things only get worse. Deep networks, boosting, and even plain linear regression with enough features say otherwise: cross the interpolation threshold and keep going, and test error comes back down, often below the classical optimum. To see why, we need the decomposition the U-curve is built from, made honest about randomness.
The bias-variance decomposition, stated precisely
For a fixed test input x, the expected squared error over both the label noise and the random draw of the training set D splits into three non-negative pieces:
E[ (y - f̂(x))^2 ]
= σ^2 (irreducible noise)
+ ( E_D[f̂(x)] - f*(x) )^2 (bias^2)
+ E_D[ ( f̂(x) - E_D[f̂(x)] )^2 ] (variance)The noise floor σ^2 is fixed; nothing you do to the model touches it. Bias is how far the average fitted function sits from the truth — it is large for models too rigid to represent f* and shrinks as capacity grows. Variance is how much f̂ jitters when you re-draw the training data. The classical story is entirely about these two trading off. Double descent is a statement about the variance term specifically: it is not monotone in capacity. It has a sharp peak, and the peak sits at one very particular place.
A solvable model: linear regression on random features
To compute those terms we need a model we can actually solve. Take random features: map each input through a fixed nonlinear featurization φ(x) ∈ R^p (think a random projection followed by a nonlinearity, the linearized stand-in for a wide network), then fit a linear map f̂(x) = φ(x)·β. Stack the training features into a design matrix Φ: [N, p] and labels y: [N]. The knob is p, the number of features (parameters); N is fixed.
The ratio γ = p / N is the whole story. γ < 1 is the under-parameterized regime: more equations than unknowns, an overdetermined system with no exact fit. γ > 1 is over-parameterized: more unknowns than equations, so infinitely many β hit the labels exactly. The boundary γ = 1, where p = N and Φ is square, is the interpolation threshold — and, we will see, the site of the catastrophe.
Two solutions on two sides of p = N
The estimator is ordinary ridgeless least squares, but its form flips at the threshold. For γ < 1 the system is overdetermined and we minimize residual error:
γ < 1 : β̂ = (Φ^T Φ)^-1 Φ^T y (ordinary least squares)
γ > 1 : β̂ = Φ^T (Φ Φ^T)^-1 y (minimum-norm interpolator)Both are the Moore-Penrose pseudoinverse β̂ = Φ^+ y, specialized to the tall and wide cases. Below the threshold we cannot fit exactly, so we take the best least-squares residual. Above it, infinitely many vectors interpolate, and the pseudoinverse picks the one with the smallest Euclidean norm ||β||. That min-norm choice is not a design decision you make consciously — it is what gradient descent from zero converges to on a least-squares loss. Hold on to it: the implicit preference for small norm is exactly what tames the second regime.
Where the variance comes from: the inverse Gram matrix
Plug the noisy labels y = Φβ* + ε into the estimator and the variance contribution is governed by the inverse of the Gram matrix. In the underparameterized case the parameter-space covariance of β̂ is
Cov(β̂) = σ^2 (Φ^T Φ)^-1Write Φ = U Σ V^T in singular values s_1 ≥ … ≥ s_r. Then (Φ^T Φ)^-1 has eigenvalues 1 / s_i^2, so the total variance scales like σ^2 Σ_i 1 / s_i^2. Every direction contributes an amount inversely proportional to the square of its singular value. A single small singular value — a direction in feature space the training data barely constrains — dumps an enormous term into the sum. The variance of the fit is therefore a magnifying glass on the smallest singular value of the design matrix. The question becomes: when is s_min smallest?
The condition number diverges exactly at the threshold
The answer comes from random matrix theory. For a tall or wide random Φ with aspect ratio γ = p/N, the eigenvalues of the normalized Gram matrix (1/N)Φ^TΦ follow the Marchenko-Pastur law, which pins the edges of the spectrum:
s_min^2 / N ≈ (1 - √γ)^2
s_max^2 / N ≈ (1 + √γ)^2
condition number κ = s_max / s_min ≈ (1 + √γ) / |1 - √γ|Look at what happens as γ → 1: (1 - √γ)^2 → 0, so s_min → 0 and κ → ∞. The design matrix becomes singular precisely when the number of parameters equals the number of samples. That is the whole mechanism: at p = N the system is square and generically invertible but barely so, its smallest singular value crushed to zero, and the σ^2 Σ 1/s_i^2 variance blows up. The interpolation-threshold peak in the double descent curve is a diverging condition number, nothing more mysterious.
Why the variance descends again past the peak
If γ = 1 is a singular-matrix catastrophe, why does γ = 4 behave well? Because past the threshold the estimator changes character. With p > N the system is underdetermined, so s_min is no longer forced toward zero — from Marchenko-Pastur, (1 - √γ)^2 grows again as γ moves away from 1, and the spectrum lifts off the origin. Just as important, the pseudoinverse now selects the minimum-norm interpolant.
Minimum norm is implicit regularization. Among the infinitely many β that fit the data, the smallest one spreads its weight thinly across many features and refuses to load huge coefficients onto the ill-conditioned directions that caused the spike. As p grows, there are ever more interpolating solutions to choose from, the min-norm one gets smoother, and both its norm and its variance fall. Bias stays low because the richer feature space can represent f*. Low bias plus falling variance is the second descent.
A worked numeric example
Take N = 100 training points, label noise σ^2 = 1, and sweep the feature count p. Using the edge formulas, watch the smallest singular value and the variance proxy σ^2 / s_min^2:
| p | γ = p/N | s_min^2 / N ≈ (1-√γ)^2 | variance ∝ 1/s_min^2 |
|---|---|---|---|
| 50 | 0.50 | 0.086 | moderate |
| 90 | 0.90 | 0.0025 | large |
| 99 | 0.99 | 0.000025 | exploding |
| 101 | 1.01 | 0.000025 | exploding |
| 200 | 2.00 | 0.172 | small |
| 400 | 4.00 | 1.000 | tiny |
The variance proxy climbs by four orders of magnitude as γ approaches 1 from either side, then collapses once γ is comfortably above 1. The test-error curve inherits that shape exactly: a gentle classical U on the left, a violent spike straddling p = N, and a long clean descent into the overparameterized regime. The peak is razor-thin, which is why you rarely see it unless you deliberately tune capacity through N.
Ridge regularization dissolves the peak
The decomposition also predicts the cure. Add an ℓ2 penalty λ||β||^2 and the estimator becomes β̂ = (Φ^TΦ + λI)^-1 Φ^T y. The inverse Gram eigenvalues turn from 1/s_i^2 into s_i^2 / (s_i^2 + λ)^2, which stays bounded even as s_i → 0.
In other words, ridge caps the contribution of the near-zero singular directions that manufacture the spike. With enough regularization the interpolation-threshold peak flattens out entirely, and the test-error curve reverts to a single monotone descent — optimally-tuned ridge regression has no double descent at all. This is a useful sanity check on the whole picture: the phenomenon is a property of ridgeless (or weakly regularized) interpolation. The same explicit regularization that a practitioner would reach for to fix an ill-conditioned linear system is exactly what erases the peak, confirming that the peak was an ill-conditioning artifact from the start.