Generalization: double descent and implicit regularization

Before this lesson: lesson 02 (training and test error, the training loop) and lesson 06 (early stopping acts like ridge). Math: row space and null space in linear-algebra lesson 02, the SVD, the pseudoinverse and minimum-norm least squares in linear-algebra lesson 08, and random matrices in linear-algebra lesson 18 (only for one folded box). From other tracks: the bias–variance decomposition and the VC dimension (ML fundamentals track), and the step-size limit of gradient descent (optimization track). Time: about 100 minutes, plus the notebook. This is a heavy lesson. The folded boxes can be skipped on a first reading.

---

Why this matters #

ResNets and transformers usually have more parameters than training examples. They are trained to (almost) zero training loss. Classical statistics predicts that such models overfit. In practice they often generalize well, and a larger model is often better. Two experimental facts make the puzzle sharp:

• Double descent. The test error, as a function of model size, can go down, up, and down again (Belkin et al. 2019; Nakkiran et al. 2019). It also happens over training time and over the number of training samples: near the peak, more data can hurt.
• Random labels. A standard CNN fits CIFAR-10 with completely random labels to 100% training accuracy (Zhang et al. 2017). The same network on the real labels reaches about 86% test accuracy.

This lesson proves what can be proved for linear models (they are also the linearization of very wide networks), and gives the experimental facts for deep networks. As of 2026 there is no complete theory for deep networks. Interviews often ask: "Why do over-parameterized networks not overfit?", "What is double descent?", "Does gradient descent find the minimum-norm solution?". In ML for science, a surrogate fitted on 100 simulation runs with about 100 features sits exactly at the dangerous point. There the fitting rule (ridge, early stopping, the start of training) decides the result more than the number of parameters.

1. Predict first #

2. The picture #

One equation, two unknowns #

Take one training example with two weights: the equation $w_1 + 2w_2 = 5$. In matrix form, $X = (1\ \ 2)$ and $y = 5$. Every point on the line $w_1 + 2w_2 = 5$ has zero loss $\mathcal{L}(w) = \tfrac12(w_1 + 2w_2 - 5)^2$. Which point does gradient descent find?

1. The gradient is $\nabla\mathcal{L}(w) = (w_1 + 2w_2 - 5)\cdot(1, 2)$. The first factor is a number, the residual $r$. So the gradient always points along $(1, 2)$, whatever $w$ is.
2. The vector $(1, 2)$ is the normal of the line (it is perpendicular to the line).
3. So gradient descent moves only along the normal. From the origin it moves on the ray through $(1, 2)$.
4. It stops where this ray meets the line: at $(1, 2)$. This is the point of the line closest to the origin, the solution with the smallest norm.

The loss does not prefer $(1, 2)$: every point of the line has loss 0. The algorithm prefers it. This is implicit regularization: the training procedure, not a penalty term, picks one of many equally good solutions. If we start at another point $w_0$, gradient descent again moves along the normal, so it stops at the point of the line closest to $w_0$.

Interactive figure: loading…

One sample with noise #

Now take one unknown and one sample: fit $y = wx$ with $x = 0.1$ and the label $y = 1 + \varepsilon$, where $\varepsilon$ is noise. The only solution is $w = y/x = 10(1 + \varepsilon)$. The noise is multiplied by $1/x = 10$. A small input number makes the fit very sensitive to noise.

The same happens with many features when the number of features $p$ equals the number of samples $n$. Then the square matrix $X$ must be inverted. Its smallest singular value is typically small, and the noise is divided by it. This point $p = n$ is the interpolation threshold: the model is just large enough to fit every training point, and it has exactly one way to do it.

Three regimes #

Keep $n$ fixed and let the number of features $p$ grow:

• $p \ll n$ (under-parameterized). Least squares averages the noise over many more equations than unknowns. More features lower the bias and raise the variance: the classical U-curve.
• $p \approx n$ (critical). There is one interpolant, and it is very sensitive to the noise. The test error peaks.
• $p \gg n$ (over-parameterized). There are many interpolants. The minimum-norm one spreads the fit over many directions. It has a small norm and amplifies the noise only a little. The test error falls again: the second descent.

Plotted against $p$, the test error goes down, up, and down: model-wise double descent (the name is from Belkin et al. 2019). Training longer moves a large model through the same three regimes in time: epoch-wise double descent (Nakkiran et al. 2019, section 3.4 below).

Interactive figure: loading…

3. Definitions and main results #

3.1 Setting #

We have $n$ training pairs $(x_i, y_i)$ with $x_i \in \mathbb{R}^p$. The $x_i$ can be raw inputs or features $\phi(x_i)$, for example random features. Stack them as rows of $X \in \mathbb{R}^{n\times p}$ (row $i$ is $x_i^\top$, as everywhere in this track), and let $y \in \mathbb{R}^n$. A linear model predicts $x^\top w$ with weights $w \in \mathbb{R}^p$. The training loss is $\mathcal{L}(w) = \tfrac12\|Xw - y\|^2$ (a sum, not a mean; lesson 06 divides by $n$).

• An interpolant is a $w$ with $Xw = y$, that is, with zero training loss. If $p > n$ and $\operatorname{rank} X = n$, interpolants exist for every $y$. They form the set $w^\star + N(X)$, a shifted copy of the null space $N(X) = \{z : Xz = 0\}$, of dimension $p - n$.
• The row space $C(X^\top)$ is the span of the rows $x_i$. It is perpendicular to $N(X)$ (linear-algebra lesson 02).
• The minimum-norm interpolant is $X^+y = X^\top(XX^\top)^{-1}y$, where $X^+$ is the pseudoinverse (linear-algebra lesson 08). It lies in the row space. It is also called ridgeless regression, because it is the limit of ridge regression as the penalty $\lambda \to 0$.
• With the SVD $X = \sum_{i=1}^{r} s_i u_i v_i^\top$ (singular values $s_1 \ge \dots \ge s_r > 0$, $r = \operatorname{rank} X$, left and right singular vectors $u_i \in \mathbb{R}^n$ and $v_i \in \mathbb{R}^p$):

Read the last formula slowly: along each direction $v_i$, the fit is the data coefficient $u_i^\top y$ divided by $s_i$. A small singular value makes a large coefficient. This one formula explains the peak of section 3.3.

Implicit regularization (or implicit bias) is the preference of a training algorithm for some solutions over others that have the same training loss, without any penalty term in the loss.

3.2 Gradient descent finds the minimum-norm interpolant #

Theorem. Let $X \in \mathbb{R}^{n\times p}$ with $p > n$ and $\operatorname{rank} X = n$. Run gradient descent $w_{k+1} = w_k - \eta X^\top(Xw_k - y)$ from $w_0 = 0$, with a step $0 < \eta < 2/s_1^2$. Then $w_k \to X^+y$, the unique interpolant with the smallest Euclidean norm.

Why this is true, in one sentence. Every gradient $X^\top(\cdot)$ is a combination of the rows of $X$, so gradient descent can never leave the row space; and the only interpolant inside the row space is the minimum-norm one.

Proof.

1. The iterates stay in the row space. $w_0 = 0$ lies in $C(X^\top)$. Every update adds $-\eta X^\top(Xw_k - y)$, which is $X^\top$ times a vector, so it lies in $C(X^\top)$. Hence $w_k = X^\top\alpha_k$ for some $\alpha_k \in \mathbb{R}^n$.
2. The residual goes to zero. Let $r_k = Xw_k - y$. Then $r_{k+1} = Xw_k - \eta XX^\top r_k - y = (I - \eta XX^\top)\,r_k$. The $n\times n$ matrix $XX^\top$ is symmetric positive definite (rank $n$), with eigenvalues $s_1^2 \ge \dots \ge s_n^2 > 0$. The eigenvalues of $I - \eta XX^\top$ are $1 - \eta s_i^2$. Because $0 < \eta s_i^2 \le \eta s_1^2 < 2$, they all lie in $(-1, 1)$. So $r_k \to 0$ geometrically.
3. The limit. $Xw_k = XX^\top\alpha_k \to y$, and $XX^\top$ is invertible, so $\alpha_k \to (XX^\top)^{-1}y$ and $w_k \to w_\infty = X^\top(XX^\top)^{-1}y$. This $w_\infty$ interpolates and lies in $C(X^\top)$.
4. It has the smallest norm. Any interpolant can be written $w = w_\infty + z$ with $Xz = 0$, so $z \in N(X)$. Since $N(X) \perp C(X^\top)$, Pythagoras gives $\|w\|^2 = \|w_\infty\|^2 + \|z\|^2 \ge \|w_\infty\|^2$, with equality only for $z = 0$. So $w_\infty = X^+y$. $\blacksquare$

The condition $\eta < 2/s_1^2$ is the usual step-size limit "$\eta < 2/$(largest eigenvalue of the Hessian)" of the optimization track: here the Hessian is $X^\top X$, and its largest eigenvalue is $s_1^2$.

Corollary (nonzero start). From any $w_0$, gradient descent converges to

the interpolant closest to $w_0$. Here $X^+X$ is the orthogonal projector onto the row space, and $(I - X^+X)w_0$ is the null-space part of $w_0$.

Proof. Split $w_0$ into its row-space part and its null-space part $z_0 = (I - X^+X)w_0$. Every update lies in the row space, so the null-space part of $w_k$ stays $z_0$ forever. The row-space part follows the proof above (step 2 does not use $w_0 = 0$), so the limit is $X^+y + z_0$. It is the interpolant closest to $w_0$, because $w_0 - w_\infty$ lies in the row space and is therefore perpendicular to $N(X)$, the directions along the solution set. $\blacksquare$

So the algorithm and its start choose the solution. Fine-tuning a pretrained model starts at the pretrained weights, not at 0; in the linearized picture it returns the solution closest to them.

What else changes the chosen solution.

• SGD. A mini-batch gradient is also a combination of rows $x_i$, so step 1 still holds. If SGD from 0 converges to an interpolant, it is $X^+y$ (Zhang et al. 2017, section 5).
• A fixed preconditioner. Gradient descent with a fixed symmetric positive definite matrix $P$, $w \leftarrow w - \eta PX^\top(Xw - y)$, is plain gradient descent in the variable $v = P^{-1/2}w$. So from 0 it converges to the interpolant that minimizes $w^\top P^{-1}w$, not $\|w\|^2$.
• Adam. Adam divides each coordinate of the step by its own running scale $\sqrt{\hat v}$. Its preconditioner changes at every step, so there is no general formula for its limit. When it converges to an interpolant, this is in general not the minimum-norm one. On the one-equation example it is exact: the gradient is always $r\,(1, 2)$, so the second coordinate of $m$ is always twice the first, and the second coordinate of $v$ four times the first. The ratio $m/\sqrt v$ is therefore the same in both coordinates: Adam moves along $(1, 1)$. From 0 it stops at $(5/3, 5/3)$, with squared norm $50/9 \approx 5.56 > 5$. (This ignores Adam's tiny $\epsilon$.)

Other parametrizations, wide networks and other losses

• Parametrization. Write the same linear model as $w = u\odot u - v\odot v$ ($\odot$ is the elementwise product) and run gradient descent on $(u, v)$ from $u = v = \alpha\mathbf{1}$, so that $w = 0$ at the start. As $\alpha \to 0$, the limit becomes the minimum-$\ell_1$ interpolant, which prefers sparse $w$; as $\alpha \to \infty$, the minimum-$\ell_2$ one (Woodworth et al. 2020). Same function class, same loss, different solution. In the notebook, with $n = 10$ samples, $p = 40$ weights and a truth with 3 nonzero weights: $\alpha = 10^{-3}$ finds $w$ within 0.03 of the truth; $\alpha = 3$ returns $X^+y$, which is 2.12 away from it.
• Wide networks. Near its initialization $\theta_0$, the outputs of a network on the training inputs are $f(X;\theta) \approx f(X;\theta_0) + J(\theta - \theta_0)$, where $J$ is the Jacobian with respect to the parameters. This is linear least squares in $\theta - \theta_0$. So gradient descent returns the smallest change $J^+\big(y - f(X;\theta_0)\big)$, by the corollary above. Very wide networks stay close to this linear model during training: the neural tangent kernel (NTK) regime.
• Logistic loss. On linearly separable data the logistic loss has no minimizer: $\|w_k\| \to \infty$. But the direction $w_k/\|w_k\|$ converges to the maximum-margin direction, the hard-margin SVM solution (Soudry et al. 2018).
• Finite steps. With a finite step $\eta$, gradient descent follows, to first order in $\eta$, the gradient flow of the modified loss $\mathcal{L} + \tfrac{\eta}{4}\|\nabla\mathcal{L}\|^2$ (Barrett and Dherin 2021; Prince, Understanding Deep Learning, eq. 9.8). A larger step adds a stronger preference for flat regions.

3.3 The peak at the interpolation threshold #

Model. Assume a true linear model $y = X\beta + \varepsilon$, where the noise $\varepsilon$ has mean 0 and covariance $\sigma^2 I$. Test inputs are $x \sim \mathcal N(0, I_p)$. The excess test risk of an estimate $\hat w$ is

because $\mathbb{E}[xx^\top] = I$. (The full test error adds the irreducible noise $\sigma^2$.)

Bias and variance of the minimum-norm fit, step by step.

1. Put $y = X\beta + \varepsilon$ into $\hat w = X^+y$: $\hat w = X^+X\beta + X^+\varepsilon$.
2. $X^+X$ is the orthogonal projector onto the row space. So $X^+X\beta$ is $\beta$ projected onto the row space, and $\hat w - \beta = -(I - X^+X)\beta + X^+\varepsilon$.
3. With the SVD, $X^+\varepsilon = \sum_i \frac{u_i^\top\varepsilon}{s_i}\,v_i$: the noise is divided by the singular values.
4. Average over the noise. The cross term has mean 0, and $\mathbb{E}\|X^+\varepsilon\|^2 = \sigma^2\operatorname{tr}(X^{+}X^{+\top}) = \sigma^2\sum_i 1/s_i^2$. So

• For $p < n$ (full column rank), $X^+X = I$: there is no bias, and the variance is $\sigma^2\operatorname{tr}\big((X^\top X)^{-1}\big)$.
• For $p > n$, the part of $\beta$ in $N(X)$ is never seen by the data: this is the bias. The variance is $\sigma^2\operatorname{tr}\big((XX^\top)^{-1}\big)$.
• The peak. The variance is a sum of $1/s_i^2$. For a random $n\times p$ matrix with i.i.d. $\mathcal N(0, 1)$ entries, the smallest nonzero singular value is about $|\sqrt p - \sqrt n|$ when $n$ and $p$ are large (linear-algebra lesson 18). At $p = n$ it is only of order $1/\sqrt n$. One tiny $s_i$ is enough to make the variance huge. Far from $p = n$, every $s_i$ stays away from 0, and the variance is small.

Interactive figure: loading…

Exact risk for Gaussian features

Let the entries of $X$ be i.i.d. $\mathcal N(0, 1)$, independent of $\varepsilon$. Two facts about inverse Wishart matrices give the means of the variance term: $\mathbb{E}[(X^\top X)^{-1}] = I_p/(n - p - 1)$ for $p \le n - 2$, and $\mathbb{E}[(XX^\top)^{-1}] = I_n/(p - n - 1)$ for $p \ge n + 2$. For the bias, use rotation invariance: the row space of $X$ is a uniformly random $n$-dimensional subspace of $\mathbb{R}^p$, so the projector $X^+X$ has mean $(n/p)I$, and $\mathbb{E}\|(I - X^+X)\beta\|^2 = (1 - n/p)\|\beta\|^2$. Together:

$$\mathbb{E} R = \begin{cases} \sigma^2\,\dfrac{p}{n - p - 1}, & p \le n - 2,\\[3mm] \|\beta\|^2\Big(1 - \dfrac{n}{p}\Big) + \sigma^2\,\dfrac{n}{p - n - 1}, & p \ge n + 2,\end{cases}$$

and $\mathbb{E} R = \infty$ for $p \in \{n - 1, n, n + 1\}$. For a fixed ratio $\gamma = p/n$ and $n \to \infty$, the risk tends to $\sigma^2\gamma/(1 - \gamma)$ for $\gamma < 1$, and to $\|\beta\|^2(1 - 1/\gamma) + \sigma^2/(\gamma - 1)$ for $\gamma > 1$ (Hastie, Montanari, Rosset and Tibshirani 2022). The same steps give the mean squared norm of the fit for $p \ge n + 2$: $\mathbb{E}\|\hat w\|^2 = (n/p)\|\beta\|^2 + \sigma^2 n/(p - n - 1)$. It falls as $p$ grows: more features give more ways to interpolate, and the shortest one gets shorter.

Three facts that keep the picture honest (Hastie et al. 2022):

1. In this well-specified model (the $p$ features contain everything the target uses), the best risk is in the under-parameterized regime, and for $p \to \infty$ the risk tends to the null risk $\|\beta\|^2$ (the risk of always predicting 0). The second descent exists, but it does not win.
2. When the model is misspecified (the target uses information that the $p$ features miss, as with random features), more features also lower the bias. Then the best risk can lie past the threshold, as in the widget above.
3. Ridge with a well-tuned $\lambda > 0$ is better than the minimum-norm fit at every $p$, and it has no peak. Ridge minimizes $\|Xw - y\|^2 + \lambda\|w\|^2$. In the SVD it replaces the factor $1/s_i$ by $s_i/(s_i^2 + \lambda)$, which is never larger than $1/(2\sqrt\lambda)$ (because $s_i^2 + \lambda \ge 2s_i\sqrt\lambda$). So a tiny $s_i$ can no longer blow up the noise.

3.4 Early stopping and epoch-wise double descent #

Lesson 06 showed that early stopping acts like ridge. Here is the same result in SVD form; it holds for any $n$ and $p$. (In lesson 06, $s$ was an eigenvalue of $X^\top X/n$; here $s_i$ is a singular value of $X$, so $s_i^2$ plays that role.)

1. Let $c_t = v_i^\top w_t$ be the coefficient of $w_t$ along the right singular vector $v_i$.
2. Since $v_i^\top X^\top = s_i u_i^\top$ and $u_i^\top Xw = s_i\,v_i^\top w$, the update $w_{t+1} = w_t - \eta X^\top(Xw_t - y)$ gives the scalar recursion $c_{t+1} = c_t - \eta s_i(s_i c_t - u_i^\top y)$.
3. Its factor is $1 - \eta s_i^2$ and its fixed point is $u_i^\top y/s_i$. From $c_0 = 0$:

The filter factor $1 - (1 - \eta s_i^2)^t$ grows from 0 to 1. Direction $i$ needs about $1/(\eta s_i^2)$ steps. Large singular directions are fitted first. The small ones, where $1/s_i$ amplifies the noise, come last. Stopping after $t$ steps acts like ridge with $\lambda \approx 1/(\eta t)$: both filters are near 1 for $s_i^2 \gg 1/(\eta t)$, and near $\eta t s_i^2$ for $s_i^2 \ll 1/(\eta t)$.

Epoch-wise double descent: one mechanism. For one group of directions, the test error over time is a U: it falls while the signal is fitted, and it rises while the noise is fitted. Now take two groups of features that learn at very different speeds. Two U-curves overlap in time:

1. The fast features fit their own signal: the test error goes down.
2. The fast features keep fitting the residual. It still holds the slow signal and the noise, which the fast weights cannot represent well. The fast weights grow in wrong directions: the test error goes up.
3. Much later, the slow directions are learned. They take the slow signal out of the residual, and the fast weights are corrected: the test error goes down again.

Interactive figure: loading…

Heckel and Yilmaz (2021) explain epoch-wise double descent in deep networks in the same way, with layers that learn at different speeds. They remove the bump by scaling the step size per layer. A practical lesson: keep the best checkpoint of the whole run, not the first local minimum of the validation curve.

3.5 Double descent in deep networks #

For a deep network, the horizontal axis is not exactly the number of parameters. Nakkiran et al. (2019) define the effective model complexity (EMC) of a training procedure (architecture, optimizer, number of epochs, data augmentation) as the largest number of samples $n$ on which it reaches an average training error of at most $\epsilon$ (they use $\epsilon = 0.1$). Their hypothesis: if the EMC is well below or well above $n$, anything that raises the EMC lowers the test error; if EMC $\approx n$ (the critical regime), it may raise it. Their experiments:

• Model-wise. ResNet18s of growing width on CIFAR-10 and CIFAR-100, 5-layer CNNs, and transformers on IWSLT'14 translation show a peak where the models just reach zero training error. Label noise (for example 15% of labels changed at random) makes the peak much stronger. Without label noise there is often only a plateau, but sometimes a peak. Data augmentation and more samples move the peak to larger models.
• Epoch-wise. A large model's test error falls, rises and falls again during training. Medium models follow a U in time, and models too small to interpolate only improve.
• Sample-wise. More data moves the threshold to larger models. So near the threshold, more data can hurt: for some transformer sizes on IWSLT'14, training on 18k samples is worse than on 4k samples.

Belkin et al. (2019) found the peak at $p = n$ random Fourier features on MNIST, and at about $nK$ parameters for a network with one hidden layer and $K$ outputs (the network must fit $nK$ numbers).

3.6 Why capacity bounds cannot explain generalization #

Uniform convergence bounds. Let $R(h)$ be the test (0–1) error of a classifier $h$ with values in $\{-1, +1\}$, and $\hat R(h)$ its training error on $n$ samples. A standard bound (Mohri, Rostamizadeh and Talwalkar, Foundations of Machine Learning, chapter 3): with probability at least $1 - \delta$ over the training set, for every $h$ in the class $\mathcal H$,

where

$S$ is the training set and the $\sigma_i = \pm 1$ are independent random signs. The Rademacher complexity $\mathfrak R_n$ measures how well the class can fit random $\pm 1$ labels on $n$ training inputs. VC bounds have the same form, with a term of order $\sqrt{d_{VC}\log(n)/n}$, where $d_{VC}$ is the VC dimension of the class. These bounds depend on the class, not on the labels and not on the algorithm.

The random-label experiment (Zhang et al. 2017). An Inception network (1.65 M parameters) on CIFAR-10 ($n = 50{,}000$ images, 10 classes), without data augmentation or weight decay:

labels	training accuracy	test accuracy
true	100%	85.75%
random	100%	9.78% (chance: 10%)

Same architecture, same optimizer, same hyperparameters. The networks also fit images of random pixels. Inception V3 reaches 95.2% top-1 training accuracy on ImageNet with random labels (91.2% even with dropout and weight decay). Random labels make training slower only by a small constant factor.

Interactive figure: loading…

The argument, step by step.

1. If the class fits every labelling of the training inputs, then for every sign pattern $\sigma$ some $h$ has $h(x_i) = \sigma_i$, so the supremum is 1. Then $\mathfrak R_n(\mathcal H) \approx 1$ and $d_{VC} \ge n$. The bound reads $R(h) \le 0 + 1 + \dots$, which says nothing, since always $R \le 1$.
2. More generally, take any bound $R(h) \le \hat R(h) + C(\mathcal H, n, \delta)$ that holds for every data distribution. It must also hold for the random-label distribution. There the trained network has $\hat R = 0$ and $R \approx 0.90$, so $C \ge 0.90$.
3. $C$ does not depend on the labels. So for the true labels, the same bound allows a test error of 0.90, while the real one is 0.14. The bound cannot tell the two cases apart, so it cannot explain the good one.
4. Bounds based on the uniform stability of the algorithm fail for the same reason: they do not look at the labels either.
5. So an explanation must depend on the data and on the solution that the algorithm actually returns (Zhang et al. 2017, section 2.2).

What the experiment does not show. It does not show that learning theory is wrong: the bounds are true, they are only too loose (vacuous) here. It does not show that data do not matter: random labels give chance test accuracy and take longer to fit. And "the network can fit random labels" is not "the network will fit noise": on real labels, gradient descent fits the simple, shared structure first (the fast directions of section 3.4) and memorizes the noise later (Arpit et al. 2017). The 50% curve in the widget shows this.

Partial answers (as of 2026). None of them is a complete explanation for deep networks.

• Implicit bias toward small-norm, max-margin or sparse interpolants (section 3.2). The norm alone does not predict the test error well (Zhang et al. note this).
• Bounds that depend on the solution found: norm-based, margin-based and PAC-Bayes bounds. Dziugaite and Roy (2017) computed non-vacuous PAC-Bayes bounds for small networks on MNIST.
• Flat minima. Solutions in wide, flat valleys tend to generalize. But flatness changes when the network is reparametrized, so it cannot be the full story.
• Benign overfitting. Fitting noisy labels exactly can be almost harmless when many directions with small variance absorb the noise (Bartlett, Long, Lugosi and Tsigler 2020).
• Inductive bias. The architecture (convolutions, attention, residual connections) matches the structure of real data.

Answers to section 1. (1) $p = 100 = n$ is by far the worst: it is the interpolation threshold. Whether $p = 1000$ beats $p = 50$ depends on the problem. In the random-features experiment of section 5 ($n = 40$, inputs in $\mathbb{R}^{10}$), the far over-parameterized fit wins clearly (test error 0.10 against the best U value 0.49). With inputs in $\mathbb{R}^{5}$ it does not quite win (0.142 against 0.125). In the well-specified Gaussian model of section 3.3, the second descent never beats the best small model. (2) Gradient descent returns the minimum-norm interpolant $X^+y$. From $w_0$ it returns $X^+y + (I - X^+X)w_0$: the null-space part of the start is kept. (3) No. The class can fit random labels, so any bound that depends only on the class must allow a test error near chance. An explanation must use the data and the algorithm.

4. Worked example #

A. The theorem by hand. One sample, two features: $X = (1\ \ 2)$, $y = 5$.

• $XX^\top = 1 + 4 = 5$, so $X^+y = X^\top(XX^\top)^{-1}y$ $= (1, 2)\cdot 5/5$ $= (1, 2)$, with $\|X^+y\|^2 = 5$.
• The null space $N(X)$ is spanned by $(2, -1)$ (check: $2 - 2 = 0$). All interpolants are $(1, 2) + c\,(2, -1)$. Their squared norm is $5 + 5c^2$ (the cross term $2c\,(1\cdot 2 + 2\cdot(-1))$ is 0), which is smallest at $c = 0$.
• $s_1^2 = \|(1, 2)\|^2 = 5$, so gradient descent converges for $\eta < 2/5 = 0.4$. Take $\eta = 0.1$. The residual is multiplied by $1 - 0.1 \cdot 5 = 0.5$ at every step:

$k$	$w_k$	residual $r_k = w_1 + 2w_2 - 5$
0	$(0, 0)$	$-5$
1	$(0, 0) - 0.1\cdot(-5)\cdot(1, 2) = (0.5, 1)$	$-2.5$
2	$(0.5, 1) - 0.1\cdot(-2.5)\cdot(1, 2) = (0.75, 1.5)$	$-1.25$
3	$(0.875, 1.75)$	$-0.625$

In general $w_k = (1 - 0.5^k)(1, 2)$: every iterate lies on the ray through $(1, 2)$, the row space.

• Start instead at $w_0 = (1, 0)$. Its null-space part is $\frac{(1, 0)\cdot(2, -1)}{\|(2, -1)\|^2}(2, -1) = \frac{2}{5}(2, -1) = (0.8, -0.4)$. Gradient descent converges to $(1, 2) + (0.8, -0.4) = (1.8, 1.6)$. Check: $1.8 + 2\cdot 1.6 = 5$ ✓, and $\|(1.8, 1.6)\|^2 = 3.24 + 2.56 = 5.8 = 5 + 0.8$: the start added $\|(0.8, -0.4)\|^2 = 0.8$.
• Adam from 0 moves along $(1, 1)$: $w = c(1, 1)$ with $c + 2c = 5$, so $c = 5/3$ and $\|w\|^2 = 50/9 \approx 5.56$.

B. The peak in numbers. Gaussian features, $n = 20$, $\sigma^2 = 0.25$, $\|\beta\|^2 = 1$, with the exact formula of section 3.3:

$p$	5	10	15	18	19–21	22	30	40	80	400
bias²	0	0	0	0		0.091	0.333	0.500	0.750	0.950
variance	0.089	0.278	0.938	4.50	$\infty$	5.00	0.556	0.263	0.085	0.013
risk	0.089	0.278	0.938	4.50	$\infty$	5.09	0.889	0.763	0.835	0.963

• $p = 18$: variance $0.25\cdot 18/(20 - 18 - 1) = 4.5$. $p = 22$: bias $1 - 20/22 = 0.091$, variance $0.25\cdot 20/(22 - 20 - 1) = 5$.
• Past the peak, the variance falls fast and the bias rises slowly. The second descent has a local minimum near $p = 40$. Then the risk climbs toward the null risk $\|\beta\|^2 = 1$.
• The best risk overall is at the smallest $p$ (0.089 at $p = 5$): this is fact 1 of section 3.3.

C. More data can hurt. Fix $p = 40$ and add samples (same formula, now as a function of $n$):

$n$	20	30	37	43	60	100
bias²	0.500	0.250	0.075	0	0	0
variance	0.263	0.833	4.625	5.000	0.526	0.169
risk	0.763	1.083	4.70	5.00	0.526	0.169

Going from 20 to 37 samples multiplies the risk by about 6, because $n$ moves toward $p$. This is sample-wise non-monotonicity. For example, at $n = 37$ the variance is $0.25\cdot 37/(40 - 37 - 1) = 4.625$.

5. In code #

The worked example A with NumPy, including Adam written by hand:

import numpy as np

X = np.array([[1.0, 2.0]])            # one equation: w1 + 2 w2 = 5
y = np.array([5.0])
X_pinv = np.linalg.pinv(X)
print(X_pinv @ y)                      # [1. 2.]  the minimum-norm interpolant X⁺y

def gd(w0, eta=0.1, steps=200):
    w = np.array(w0, dtype=float)
    for _ in range(steps):
        w = w - eta * X.T @ (X @ w - y)  # gradient of ½‖Xw − y‖²
    return w

print(gd([0, 0], steps=3))            # [0.875 1.75 ]  = (1 − 0.5³)(1, 2)
w_lim = gd([1, 0])
P_null = np.eye(2) - X_pinv @ X        # projector onto the null space N(X)
print(w_lim, X_pinv @ y + P_null @ np.array([1.0, 0.0]))   # [1.8 1.6] [1.8 1.6]
print(round(w_lim @ w_lim, 3))        # 5.8 = 5 + 0.8

def adam(w0, lr=0.1, decay=0.995, steps=1500, b1=0.9, b2=0.999, eps=1e-8):
    w, m, v = np.array(w0, dtype=float), np.zeros(2), np.zeros(2)
    for t in range(1, steps + 1):
        g = X.T @ (X @ w - y)
        m, v = b1 * m + (1 - b1) * g, b2 * v + (1 - b2) * g**2
        w = w - lr * decay ** (t - 1) * (m / (1 - b1**t)) / (np.sqrt(v / (1 - b2**t)) + eps)
    return w

w_adam = adam([0, 0])
print(w_adam.round(4), (w_adam @ w_adam).round(3))   # [1.6667 1.6667] 5.556

The theorem with PyTorch's own optimizers, on a random problem with 5 equations and 12 unknowns:

import torch

torch.manual_seed(0)
n, p = 5, 12
X = torch.randn(n, p, dtype=torch.float64)
y = torch.randn(n, dtype=torch.float64)
w_mn = torch.linalg.pinv(X) @ y                          # minimum-norm interpolant X⁺y
P_row = torch.linalg.pinv(X) @ X                         # projector onto the row space
eta = 1.0 / torch.linalg.matrix_norm(X, 2).item() ** 2   # 1/s_1², inside (0, 2/s_1²)

def train(w0, make_opt, steps):
    w = w0.clone().requires_grad_(True)
    opt, sched = make_opt([w])
    for _ in range(steps):
        opt.zero_grad()
        (0.5 * ((X @ w - y) ** 2).sum()).backward()
        opt.step()
        if sched: sched.step()
    return w.detach()

sgd = lambda ps: (torch.optim.SGD(ps, lr=eta), None)     # full-batch gradient descent
zero = torch.zeros(p, dtype=torch.float64)
w_gd = train(zero, sgd, 3000)
print(f"GD from 0:  residual {(X @ w_gd - y).abs().max():.0e}, |w - X⁺y| {(w_gd - w_mn).abs().max():.0e}")

w0 = torch.randn(p, dtype=torch.float64)
w_gd0 = train(w0, sgd, 3000)
print(f"GD from w0: |w - (X⁺y + (I - X⁺X)w0)| {(w_gd0 - (w_mn + w0 - P_row @ w0)).abs().max():.0e}")

def adam(ps):
    opt = torch.optim.Adam(ps, lr=1e-2)
    return opt, torch.optim.lr_scheduler.ExponentialLR(opt, gamma=0.9995)   # decay, so it settles
w_adam = train(zero, adam, 20000)
print(f"Adam from 0: residual {(X @ w_adam - y).abs().max():.0e}, ‖w‖ {w_adam.norm():.3f} "
      f"vs ‖X⁺y‖ {w_mn.norm():.3f}, null-space part {(w_adam - P_row @ w_adam).norm():.3f}")
# GD from 0:  residual 1e-16, |w - X⁺y| 3e-16
# GD from w0: |w - (X⁺y + (I - X⁺X)w0)| 1e-15
# Adam from 0: residual 3e-18, ‖w‖ 0.463 vs ‖X⁺y‖ 0.376, null-space part 0.270

Adam also interpolates, but it leaves the row space and returns a longer solution.

The exact risk formula of section 3.3 against a simulation (the table of part B):

import numpy as np

rng = np.random.default_rng(0)
n, s2 = 20, 0.25                                         # n = 20, σ² = 0.25, ‖β‖² = 1
exact = lambda p: s2 * p / (n - p - 1) if p < n else (1 - n / p) + s2 * n / (p - n - 1)
for p in [10, 15, 30, 40, 80]:
    risk = []
    for _ in range(4000):
        beta = rng.standard_normal(p); beta /= np.linalg.norm(beta)
        X = rng.standard_normal((n, p))
        w = np.linalg.lstsq(X, X @ beta + np.sqrt(s2) * rng.standard_normal(n), rcond=None)[0]
        risk.append(np.sum((w - beta) ** 2))             # excess test risk for x ~ N(0, I)
    print(p, round(np.mean(risk), 3), round(exact(p), 3))
# 10 0.278 0.278
# 15 0.926 0.938
# 30 0.888 0.889
# 40 0.762 0.763
# 80 0.834 0.835

Model-wise double descent with frozen random ReLU features ($n = 40$; medians over 20 draws). np.linalg.lstsq returns the minimum-norm solution when $p > n$:

import numpy as np
from sklearn.linear_model import Ridge

ps = [2, 5, 10, 20, 30, 36, 40, 44, 50, 100, 400, 2000]

def random_features_curve(d, n=40, sigma=0.3, reps=20, seed=0):
    """Median test MSE of min-norm least squares (and ridge) on p frozen random ReLU features."""
    rng = np.random.default_rng(seed)
    beta = rng.standard_normal(d) / np.sqrt(d)
    f = lambda x: x @ beta + 0.3 * np.tanh(x @ beta)              # a mildly nonlinear target
    out = {p: [] for p in ps}
    for _ in range(reps):                                          # new data and features each time
        x_tr, x_te = rng.standard_normal((n, d)), rng.standard_normal((2000, d))
        y_tr = f(x_tr) + sigma * rng.standard_normal(n)
        W, b = rng.standard_normal((max(ps), d)) / np.sqrt(d), rng.uniform(-1, 1, max(ps))
        for p in ps:
            phi = lambda x: np.maximum(x @ W[:p].T + b[:p], 0) / np.sqrt(p)   # p random ReLU features
            F = phi(x_tr)
            w = np.linalg.lstsq(F, y_tr, rcond=None)[0]                     # minimum-norm least squares
            ridge = Ridge(alpha=1e-2, fit_intercept=False).fit(F, y_tr)
            out[p].append([np.mean((phi(x_te) @ w - f(x_te)) ** 2),
                           np.mean((ridge.predict(phi(x_te)) - f(x_te)) ** 2),
                           np.linalg.svd(F, compute_uv=False)[min(n, p) - 1],  # smallest nonzero s_i
                           np.linalg.norm(w)])
    return {p: np.median(out[p], axis=0) for p in ps}

res = random_features_curve(d=10)
print("   p  min-norm   ridge    s_min    ‖w‖")
for p in ps:
    print(f"{p:4d} {res[p][0]:9.3f} {res[p][1]:7.3f} {res[p][2]:8.4f} {res[p][3]:6.1f}")
res5 = random_features_curve(d=5)
best = min((p for p in ps if p < 40), key=lambda p: res5[p][0])
print(f"d = 5: best p < n is p = {best} with {res5[best][0]:.3f}; p = 2000 gives {res5[2000][0]:.3f}")
#    p  min-norm   ridge    s_min    ‖w‖
#    2     0.796   0.796   2.0779    0.7
#    5     0.734   0.730   0.6498    2.5
#   10     0.576   0.571   0.2699    5.0
#   20     0.488   0.394   0.0879   12.6
#   30     0.824   0.291   0.0277   20.3
#   36     2.077   0.252   0.0111   41.2
#   40    14.372   0.224   0.0017  108.7
#   44     1.824   0.217   0.0113   43.1
#   50     0.837   0.176   0.0214   26.4
#  100     0.208   0.111   0.0582   14.1
#  400     0.110   0.072   0.0921   10.5
# 2000     0.099   0.063   0.0983   10.2
# d = 5: best p < n is p = 10 with 0.125; p = 2000 gives 0.142

All three regimes appear: a U with its bottom at $p = 20$, a peak at $p = n = 40$ where the smallest singular value is tiny and the norm of $w$ is huge, and a second descent far below the U. A small ridge penalty removes the peak. With $d = 5$ the second descent does not beat the U.

6. Where it shows up #

• Deep learning practice. Vision models are routinely trained to (almost) zero training error, and away from the critical regime bigger is better. Label noise makes the critical regime dangerous; there, tuned weight decay and ensembling help (Nakkiran et al. 2019 show the effect of ensembling). Language-model pretraining often uses each token about once, so the model does not interpolate its data. The peak matters mainly when a model is trained for many epochs on a fixed, small or noisy dataset: fine-tuning, and most scientific data.
• ML for science. Emulators are often fitted on $10^2$–$10^3$ simulation runs with random features, kernel ridge regression or Gaussian processes. A Gaussian-process posterior mean with zero noise variance is the minimum-norm interpolant in the kernel's function space (the RKHS), the kernel version of the theorem of section 3.2. A nonzero noise variance (a "nugget") gives kernel ridge regression. The $n$ that counts is the number of independent runs, not the number of correlated grid points. When new runs move $n$ toward $p$, re-tune the regularizer.
• Fine-tuning a foundation model (for example a weather or Earth-system model) on a few hundred events is deep in the over-parameterized regime. In the linearized picture, gradient descent returns the solution closest to the pretrained weights (the nonzero-start corollary).
• Early stopping. Because of epoch-wise double descent, the validation curve can rise and then fall again. Keep the best checkpoint of the whole run, not the first local minimum.

7. Common mistakes #

• "Double descent breaks the bias–variance decomposition." The decomposition holds at every $p$. What fails is the idea that the variance always grows with the number of parameters: for the minimum-norm estimator it peaks at $p \approx n$ and then falls.
• "The peak is where the number of parameters equals the number of samples." For least squares with one output it is at $p = n$ features; with $K$ outputs, at about $nK$ parameters. For deep networks it is where the whole training procedure just reaches zero training error (EMC $\approx n$). This point moves with epochs, augmentation and label noise.
• "Past the threshold, bigger always wins." In the well-specified Gaussian model the second descent never beats the best small model, and a tuned ridge beats the ridgeless fit at every size.
• "Gradient descent finds the minimum-norm parameters of any network." The theorem is about linear least squares (and linearized networks) from $w_0 = 0$. The implicit bias changes with the start, the parametrization, the optimizer (Adam leaves the row space) and the loss.
• "The VC bound is just loose by a constant." Here it is vacuous, and any bound that ignores the labels and the algorithm must be vacuous. This does not mean that the data do not matter: random labels give chance test accuracy and are fitted more slowly.
• "More data always helps." Near the threshold, more samples without re-tuning the regularizer can make the test error much worse (part C of section 4).

8. Check yourself #

1. $X = (1\ \ 1\ \ 0)$, $y = 4$. (a) Give the minimum-norm interpolant. (b) Where does gradient descent converge from $w_0 = (1, 0, 3)$?
   Answer

   (a) $XX^\top = 2$, so $X^+y = (1, 1, 0)\cdot 4/2 = (2, 2, 0)$. (b) $N(X)$ is spanned by $(1, -1, 0)$ and $(0, 0, 1)$. The null-space part of $w_0$ is $\tfrac12(1, -1, 0) + (0, 0, 3) = (0.5, -0.5, 3)$, so the limit is $(2, 2, 0) + (0.5, -0.5, 3) = (2.5, 1.5, 3)$. Check: $2.5 + 1.5 = 4$.
2. Same $X$ and $y$. For which steps $\eta$ does gradient descent converge? What happens at $\eta = 1$ from $w_0 = 0$?
   Answer

   $s_1^2 = \|(1, 1, 0)\|^2 = 2$, so $0 < \eta < 1$. At $\eta = 1$ the residual factor is $1 - 2 = -1$: $w$ jumps between $(4, 4, 0)$ and $(0, 0, 0)$, the residual alternates between $+4$ and $-4$, and it never converges. (The average of the two points is the minimum-norm solution.)
3. Gradient descent from 0 with $\eta = 0.01$ runs for $t = 10{,}000$ steps. What fraction of the minimum-norm coefficient is fitted along a direction with $s = 1$, and along one with $s = 0.05$? Which ridge $\lambda$ gives about the same?
   Answer

   The filter factor is $1 - (1 - \eta s^2)^t \approx 1 - e^{-\eta t s^2}$: about 1 for $s = 1$, and $1 - e^{-0.25} = 0.221$ for $s = 0.05$. Ridge with $\lambda \approx 1/(\eta t) = 0.01$ gives $s^2/(s^2 + \lambda) = 0.990$ and $0.200$.
4. In the Gaussian model with $n = 20$, $\sigma^2 = 0.25$ and $\|\beta\|^2 = 1$, compute the risk at $p = 25$. Is it on the rising or the falling side of the peak?
   Answer

   Bias $1 - 20/25 = 0.2$, variance $0.25\cdot 20/(25 - 20 - 1) = 1.25$, risk $1.45$. It is past the peak ($p > n$), on the falling side: from 5.09 at $p = 22$ to 0.889 at $p = 30$.
5. Zhang et al.: true labels give 85.75% test accuracy, random labels 9.78%, both with 100% training accuracy. Show that no bound $R(h) \le \hat R(h) + C(\mathcal H, n, \delta)$ can explain the first number.
   Answer

   The bound must hold for every data distribution, including the random-label one. There $\hat R = 0$ and $R \approx 0.90$, so $C \ge 0.90$. $C$ does not depend on the labels, so for the true labels it allows a test error of 0.90 or more, while the real one is 0.14. Such a bound cannot separate the two cases; an explanation must use the data and the solution the algorithm returns.
6. Does Adam from $w = 0$ on underdetermined least squares return the minimum-norm interpolant?
   Answer

   In general, no. Adam divides each coordinate of the step by its own $\sqrt{\hat v}$, so the step leaves the row space $C(X^\top)$ and the limit has a null-space part. A fixed preconditioner $P$ would give the minimizer of $w^\top P^{-1}w$; Adam's preconditioner changes over time, so there is no closed form. In section 5: norm 0.463 against 0.376, null-space part 0.270. On $w_1 + 2w_2 = 5$ it returns $(5/3, 5/3)$ instead of $(1, 2)$.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the minimum-norm interpolant of $Xw = y$ (with $p > n$, rank $n$)?

   Show answer

2. Why does gradient descent from 0 on $\tfrac12\|Xw - y\|^2$ end at the minimum-norm solution?

   Show answer

3. Where does gradient descent on underdetermined least squares converge from a start $w_0$?

   Show answer

4. What is implicit regularization?

   Show answer

5. What is the interpolation threshold, and why is the test error high there?

   Show answer

6. Bias and variance of the minimum-norm fit for $y = X\beta + \varepsilon$, noise variance $\sigma^2$?

   Show answer

7. What is model-wise double descent?

   Show answer

8. What causes epoch-wise double descent in a linear model?

   Show answer

9. How does ridge change the noise factor $1/s_i$ of the minimum-norm fit?

   Show answer

10. What does the random-label experiment (Zhang et al. 2017) show?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "What is double descent?" Test error against model size: a U, a peak at the interpolation threshold, and a second descent. It also appears over epochs and over the number of samples. It is strongest with label noise, and tuned regularization removes it.
• "Why is there a peak?" At $p \approx n$ there is exactly one interpolant. $X^+$ divides the noise by the singular values, and the smallest one, about $|\sqrt p - \sqrt n|$, is near 0.
• "Prove that gradient descent finds the minimum-norm solution." The iterates stay in $C(X^\top)$; the residual contracts by $I - \eta XX^\top$; the only interpolant in $C(X^\top)$ is $X^+y$; Pythagoras gives minimality.
• "Why do huge networks generalize if VC theory says they should not?" Capacity bounds ignore the data and the algorithm, and the random-label experiment shows they are vacuous here. Current partial answers: implicit bias, norm, margin and PAC-Bayes bounds, benign overfitting, inductive bias of the architecture.
• "Is more data always better?" Not near the threshold without re-tuning the regularizer.
• "Is early stopping like L2 regularization?" For least squares, yes: the gradient-descent filter $1 - (1 - \eta s^2)^t$ is close to the ridge filter $s^2/(s^2 + \lambda)$ with $\lambda \approx 1/(\eta t)$.
• "Does Adam have the same implicit bias as SGD?" No. Adam's per-coordinate scaling leaves the row space, so on least squares it does not return the minimum-norm interpolant in general.

10. Practice #

The notebook for this lesson runs in Colab. Each exercise has a check cell and a folded solution.

1. Gradient descent from 0 and from $w_0$ on underdetermined least squares, compared with np.linalg.pinv and the formula $X^+y + (I - X^+X)w_0$.
2. The double-descent curve with random ReLU features, with and without ridge.
3. The smallest singular value of a random $n \times p$ matrix against $p/n$, compared with $|\sqrt p - \sqrt n|$.
4. A small MLP on the scikit-learn digits with random labels: 100% training accuracy, chance test accuracy.
5. Adam against gradient descent on the same problem: the distance to the minimum-norm solution.

By hand: redo part A of section 4 with $w_0 = (0, 1)$ (answer: $(0.6, 2.2)$), and part B with $\sigma^2 = 1$.

11. Watch and read #

Statistical Learning: 10.7 Interpolation and Double DescentStanford Online · Statistical Learning (ISLR) course · 11 min

A short lecture from the Stanford statistical-learning course: interpolation, minimum-norm fits and the double-descent curve.

What the Books Get Wrong about AI [Double Descent]Welch Labs · 34 min

A visual explanation of double descent and of why the classical bias–variance picture is incomplete.

Mikhail Belkin - From classical bias-variance trade-off to double descentICMNS Digital · Mikhail Belkin · 31 min

The talk by the author who named double descent: interpolation, the peak, and kernel methods.

Benign overfitting- Peter Bartlett, UC BerkleyThe Alan Turing Institute · Peter Bartlett · 42 min

When fitting noisy labels exactly is almost harmless; research level, after section 3.6.

• Simon Prince, Understanding Deep Learning, free PDF: §8.4 (double descent, with chapter notes on effective model complexity and epoch-wise double descent), §9.2 (implicit regularization of gradient descent and SGD, eqs. 9.8–9.9), §20.2 (fitting random data) and §20.4 (factors that determine generalization). The companion notebook 8.3 Double descent trains networks of growing width on MNIST-1D with label noise (it needs a GPU for about 30 minutes).
• Belkin, Hsu, Ma and Mandal, Reconciling modern machine learning practice and the bias-variance trade-off (PNAS, 2019).
• Nakkiran et al., Deep Double Descent: Where Bigger Models and More Data Hurt (2019): §2 (effective model complexity), §5 (model-wise), §6 (epoch-wise), §7 (sample-wise).
• Zhang, Bengio, Hardt, Recht and Vinyals, Understanding deep learning requires rethinking generalization (2017): §2 (random labels and pixels, Table 1), §4 (a two-layer ReLU network with $2n + d$ weights can fit any labels of $n$ points in $\mathbb{R}^d$), §5 (implicit regularization in linear models).
• Hastie, Montanari, Rosset and Tibshirani, Surprises in High-Dimensional Ridgeless Least Squares Interpolation (Annals of Statistics, 2022): the exact risk curves of section 3.3.
• Bartlett, Long, Lugosi and Tsigler, Benign overfitting in linear regression (PNAS, 2020).
• Arpit et al., A Closer Look at Memorization in Deep Networks (ICML, 2017): networks learn simple patterns first and memorize noise later.
• Heckel and Yilmaz, Early Stopping in Deep Networks: Double Descent and How to Eliminate it (ICLR, 2021).
• Soudry et al., The Implicit Bias of Gradient Descent on Separable Data (JMLR, 2018); Woodworth et al., Kernel and Rich Regimes in Overparametrized Models (COLT, 2020).
• Mohri, Rostamizadeh and Talwalkar, Foundations of Machine Learning, 2nd edition, free PDF: chapter 3 (Rademacher complexity and VC dimension).