Double Descent In Ridgeless Regression

Why Risk Can Rise and Fall as Model Size Grows

Increasing the number of fitted coefficients can first improve prediction, then make it dramatically worse near the point where the training data can just be interpolated, and then reduce risk again as still more coefficients are added. This pattern is called double descent. The interpolation threshold matters because the fitted linear system becomes poorly conditioned there: small noise in the responses can produce large changes in the estimated coefficients.

To see this mechanism without relying on a complicated model, we will calculate risk in a Gaussian linear model. The calculation also makes an important distinction: below the threshold, the usual estimator is least squares and need not interpolate; above it, the minimum-norm estimator interpolates exactly.

Definition: The prediction risk of an estimator $\hat\beta$, for a fresh covariate vector $x_0$, is $E[(x_0^\mathsf{T}\hat\beta-x_0^\mathsf{T}\beta)^2]$, where the expectation is over the training data and the independent test covariate. The interpolation threshold in a linear model with $n$ observations is the transition around $p=n$, where $p$ is the number of fitted coefficients.

A Gaussian Model for the Risk Curve

Let $x_i\in\mathbb R^D$ have independent $N(0,1)$ coordinates, independently across observations, and let $y_i=x_i^\mathsf{T}\beta+\epsilon_i$, where the coefficient vector $\beta\in\mathbb R^D$ is fixed and $\epsilon_i\sim N(0,\sigma^2)$ independently of the covariates. Fit only the first $p$ coordinates, with $p\leq D$. The test vector $x_0$ is an independent draw from the same covariate distribution.

Write $X_p$ for the $n\times p$ matrix of the first $p$ training covariates. When $p<n$, define $\hat\beta_{1:p}$ to be the least-squares estimator; if $X_p$ has full column rank, it is $(X_p^\mathsf{T}X_p)^{-1}X_p^\mathsf{T}y$. When $p>n$, $X_p$ has full row rank with probability one, and define the fitted coefficient vector to be the minimum-Euclidean-norm solution that satisfies $X_p\hat\beta_{1:p}=y$: $X_p^\mathsf{T}(X_pX_p^\mathsf{T})^{-1}y$. When $p=n$, define it as the unique solution $\hat\beta_{1:p}=X_p^{-1}y$, which exists with probability one. These two definitions should not be conflated. For $p<n$, noisy responses generally do not lie in the column space of $X_p$, so exact interpolation is generally impossible.

The estimator is embedded in $\mathbb R^D$ by setting its unfit coordinates to zero. Define the omitted signal energy $\tau_p^2=\sum_{j>p}\beta_j^2$, the included signal energy $b_p^2=\sum_{j\leq p}\beta_j^2$, and $v_p=\sigma^2+\tau_p^2$. The omitted signal contributes random variation to the responses: conditional on the fitted design, its contribution is independent of $X_p$, Gaussian, and has variance $\tau_p^2$ in each observation.

Theorem: In this model, the expected prediction risk is $$ R_p = \begin{cases} \displaystyle \tau_p^2+v_p\frac{p}{n-p-1}, & p<n-1,\\[6pt] \displaystyle \tau_p^2+b_p^2\left(1-\frac{n}{p}\right)+v_p\frac{n}{p-n-1}, & p>n+1. \end{cases} $$ At the remaining integer dimensions $p=n-1,n,n+1$, the expected risk is infinite if $v_p>0$. The displayed finite-risk formulas apply when their denominators are positive.

Deriving the Two Risk Formulas

Because $x_0\sim N(0,I_D)$, any fixed coefficient error $a\in\mathbb R^D$ satisfies $E[(x_0^\mathsf{T}a)^2\mid a]=\|a\|^2$. Thus prediction risk is the expected squared Euclidean coefficient error. Below the threshold, the response can be written $y=X_p\beta_{1:p}+\eta$, where $\eta$ is independent of $X_p$ and has covariance $v_p I_n$. The least-squares coefficient error on the included coordinates is $(X_p^\mathsf{T}X_p)^{-1}X_p^\mathsf{T}\eta$; the omitted coordinates contribute squared error $\tau_p^2$.

The following calculation gives the expected trace needed for both sides of the threshold. Let $Z$ be an $m\times q$ matrix with independent $N(0,1)$ entries, where $m>q+1$. For column $j$, condition on the other $q-1$ columns and let $P_{-j}$ be the orthogonal projection onto their span. The Schur-complement identity gives $((Z^\mathsf{T}Z)^{-1})_{jj}=[z_j^\mathsf{T}(I-P_{-j})z_j]^{-1}$. The projection $I-P_{-j}$ has rank $m-q+1$, and $z_j$ is an independent standard Gaussian vector. Therefore the quadratic form in brackets has a chi-squared distribution with $m-q+1$ degrees of freedom. For $U\sim\chi^2_k$, its density yields $E[U^{-1}]=1/(k-2)$ for $k>2$, and this expectation is infinite for $k\leq2$. Summing the diagonal expectations gives $$ E[\operatorname{tr}((Z^\mathsf{T}Z)^{-1})]=\frac{q}{m-q-1},\qquad m>q+1. $$

For $p<n$, conditional on $X_p$, the expected squared norm of the least-squares noise contribution is $v_p\operatorname{tr}((X_p^\mathsf{T}X_p)^{-1})$. Applying the trace calculation with $m=n,q=p$, then adding the omitted signal error, gives the first risk formula.

For $p>n$, put $P=X_p^\mathsf{T}(X_pX_p^\mathsf{T})^{-1}X_p$, the orthogonal projection onto the row space of $X_p$. The included-coordinate coefficient error is $-(I-P)\beta_{1:p}+X_p^\mathsf{T}(X_pX_p^\mathsf{T})^{-1}\eta$. The cross term has conditional expectation zero because $\eta$ is independent with mean zero. Its noise contribution has conditional covariance $v_p X_p^\mathsf{T}(X_pX_p^\mathsf{T})^{-2}X_p$; by cyclicity of trace, the trace of this matrix is $v_p\operatorname{tr}((X_pX_p^\mathsf{T})^{-1})$. Applying the trace calculation to $X_p^\mathsf{T}$, an $p\times n$ standard Gaussian matrix, gives expected noise contribution $v_p n/(p-n-1)$ when $p>n+1$.

Finally, Gaussian rotational invariance makes the row space a uniformly oriented $n$-dimensional subspace of $\mathbb R^p$. Consequently $E[P]=(n/p)I_p$: the expectation is invariant under every orthogonal change of coordinates and has trace $n$. Hence $E[\|(I-P)\beta_{1:p}\|^2]=b_p^2(1-n/p)$. Adding the omitted-coordinate error proves the overparameterized formula. The same chi-squared calculation shows that the relevant inverse-trace expectation diverges at $p=n-1,n,n+1$. At those dimensions the design matrices have the required ranks with probability one, but the expected risk is infinite when $v_p>0$. $\square$

Worked Example: Risk Before the Interpolation Threshold

Take $n=100$, $\sigma^2=1$, and let the signal have $\beta_j=1$ for $j=1,2,3,4$ and $\beta_j=0$ thereafter. Thus the total signal energy is 4. Compare $p=1$, $p=2$, $p=4$, and $p=50$, all below $n$.

For $p=1$, the omitted signal energy is $\tau_1^2=3$, so $v_1=1+3=4$. The risk is $3+4(1)/(100-1-1)=3+4/98\approx3.0408$.

For $p=2$, $\tau_2^2=2$, $v_2=3$, and the risk is $2+3(2)/(100-2-1)=2+6/97\approx2.0619$.

For $p=4$, all signal coordinates are included, so $\tau_4^2=0$ and $v_4=1$. The risk is $0+4/(100-4-1)=4/95\approx0.0421$.

For $p=50$, the risk is $50/(100-50-1)=50/49\approx1.0204$. Thus risk first falls as the fitted model captures signal, then rises as the number of noisy degrees of freedom grows. These are expected risks, not claims that every realized dataset follows this exact ordering.

Worked Example: The Peak and the Overparameterized Side

Continue with $n=100$, noise variance 1, and signal energy 4 concentrated in the first four coordinates. At $p=90$, the underparameterized formula gives $90/(100-90-1)=10$. At $p=99=n-1$, the expected risk is infinite because the denominator in the inverse-Wishart trace calculation is zero. At $p=100$, the square Gaussian design is invertible with probability one, but the expected squared coefficient error is still infinite. At $p=101$, the overparameterized inverse-trace expectation is infinite as well.

At $p=102$, the finite overparameterized formula applies. Its bias term is $4(1-100/102)=8/102\approx0.0784$, and its noise term is $100/(102-100-1)=100$. The total risk is approximately $100.0784$. At $p=110$, the terms are $4(1-100/110)=4/11\approx0.3636$ and $100/9\approx11.1111$, giving risk $11.4747$. At $p=150$, they are $4(1-100/150)=4/3$ and $100/49\approx2.0408$, giving risk $3.3741$. The noise-driven risk descends sharply after the threshold, although the bias from the minimum-norm solution's unobserved directions increases.

The Proportional-Limit Curve

The finite-sample formulas show explicitly why expected risk can be infinite at dimensions very close to $n$. A second view takes $n,p$ large together, with the ratio $\gamma=p/n$ approaching a fixed value. Suppose $p$ eventually includes all nonzero coordinates of the fixed signal, so $\tau_p^2=0$ and $b_p^2=B^2=\|\beta\|^2$. For fixed $0<\gamma<1$, the underparameterized formula tends to $\sigma^2\gamma/(1-\gamma)$. For fixed $\gamma>1$, the overparameterized formula tends to $B^2(1-1/\gamma)+\sigma^2/(\gamma-1)$. These follow by dividing the finite-sample denominators by $n$.

Key takeaway: In this Gaussian model, the limiting risk diverges as $\gamma$ approaches 1 from either side when $\sigma^2>0$. The curve rises toward the interpolation threshold and has a singular peak there; beyond the threshold, the noise contribution decreases, while the minimum-norm signal bias increases toward $B^2$. The balance of those terms determines the shape of the second descent.

Worked Example: Comparing Two Overparameterized Sizes

Let the signal energy be $B^2=4$, the noise variance be $\sigma^2=1$, and consider the proportional-limit risk for $\gamma>1$: $R(\gamma)=4(1-1/\gamma)+1/(\gamma-1)$. At $\gamma=2$, the bias term is $4(1-1/2)=2$, and the noise term is $1/(2-1)=1$, so $R(2)=3$.

At $\gamma=4$, the bias term is $4(1-1/4)=3$, and the noise term is $1/(4-1)=1/3$, so $R(4)=10/3\approx3.3333$. More parameters have reduced the noise term from 1 to $1/3$, but increased the bias term from 2 to 3; the total risk is higher at $\gamma=4$. Thus the post-threshold branch need not decrease monotonically or tend to zero.

Interpretation and a Common Mistake

The peak is tied to the instability of solving a nearly square random system, not to a universal rule that every model must have a particular risk curve. The precise curve depends on the signal, noise, design, and estimator. In this example, adding useful coordinates can reduce omitted-signal bias below the threshold; near the threshold, variance becomes large; and above it, the minimum-norm choice leaves some signal outside the row space while its noise amplification falls as $p$ grows.

Mistake: calling every below-threshold fit an interpolator. With $p<n$ and positive noise, $y$ is generally not in the column space of $X_p$. The estimator is least squares, which minimizes residual sum of squares; it does not generally satisfy $X_p\hat\beta=y$.
Mistake: treating the interpolation threshold as an ordinary finite-risk point. A full-rank square Gaussian matrix is invertible with probability one, but that does not imply a finite expected inverse squared norm. Near-square inverse moments can diverge, as the chi-squared calculation shows.
Mistake: assuming the overparameterized risk must keep falling. In the derived formula the noise term falls with $p$, but the bias term can rise. Their sum, not either term alone, determines prediction risk.

Check Your Understanding

Use the Gaussian model and distinguish finite-sample statements from proportional-limit statements.

Why is least squares, rather than exact interpolation, the usual estimator when $p<n$ and the responses contain continuous noise?
For $p<n-1$, identify the omitted-signal and estimation-variance contributions to the risk formula.
Why can the expected risk be infinite at a design dimension where the Gaussian design matrix is nevertheless full rank with probability one?
In the overparameterized formula, which term decreases with $p$ when the signal is fully included, and which term increases?
What assumptions are needed for the proportional-limit formulas stated in the tutorial?