On the effect of noise on fitting linear regression models

This study aims to improve our understanding of sparsity and overfitting when fitting linear models by exploring both empirically and theoretically the effect of including either irrelevant predictors (increasing $d$) or irrelevant observations (increasing $n$). Specifically, suppose that we arrange the training data in a response vector ${\mathbf{y}}=[y_{1},\ldots,y_{n_{0}}]^{T}\in\mathbb{R}^{n_{0}}$ and a $n_{0}\times d_{0}$ matrix of predictors $\mathbf{X}=[{\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{n_{0}}]^{T}$, such that ${\mathbf{x}}_{i}\sim\mbox{independent}N_{d_{0}}({\mbox{\boldmath$0$}}_{d_{0}},{\mathbf{I}}_{d_{0}})$, where ${\mbox{\boldmath$0$}}_{d_{0}}$ is the $d_{0}$-vector of zeros and ${\mathbf{I}}_{d_{0}}$ is the $d_{0}\times d_{0}$ identity matrix, and ${\mathbf{y}}={\mathbf{X}}{\mbox{\boldmath$\beta$}}_{0}+{\mathbf{e}}$ with $\mathbf{e}\sim N_{n_{0}}({\mbox{\boldmath$0$}}_{n_{0}},\sigma^{2}{\mathbf{I}}_{n_{0}})$, where ${\mbox{\boldmath$\beta$}}_{0}\in\mathbb{R}^{d_{0}}$ is an unknown regression parameter and $\sigma^{2}>0$ is an unknown variance parameter. We consider two sequences of training sets $\mathcal{D}$: I (Adding predictors; $d$ increases with $n=n_{0}$ fixed) for $d\leq d_{0}$, we set $\mathcal{D}^{(d)}=({\mathbf{y}},{\mathbf{X}}^{(d)})$, where ${\mathbf{X}}^{(d)}$ contains the first $d$ columns of ${\mathbf{X}}$, and for $d>d_{0}$, we set $\mathcal{D}^{(d)}=({\mathbf{y}},[{\mathbf{X}},{\mathbf{Z}}])$, where ${\mathbf{Z}}$ is an $n_{0}\times(d-d_{0})$ matrix with independent $N_{d-d_{0}}({\mbox{\boldmath$0$}}_{d-d_{0}},{\mathbf{I}}_{d-d_{0}})$ rows that are independent of ${\mathbf{X}}$ and ${\mathbf{e}}$; II (Adding observations; $n$ increases with $d=d_{0}$ fixed) for $n\leq n_{0}$, we set $\mathcal{D}_{(n)}=({\mathbf{y}}_{(n)},{\mathbf{X}}_{(n)})$, where ${\mathbf{y}}_{(n)}$ contains the first $n$ elements of ${\mathbf{y}}$ and ${\mathbf{X}}_{(n)}$ contains the first $n$ rows of ${\mathbf{X}}$, and for $n>n_{0}$, we set $\mathcal{D}_{(n)}=([{\mathbf{y}}^{T},{\mbox{\boldmath$0$}}_{n-n_{0}}^{T}]^{T},[{\mathbf{X}}^{T},{\mathbf{W}}^{T}]^{T})$, where ${\mathbf{W}}$ is an $(n-n_{0})\times d_{0}$ matrix with independent $N_{d_{0}}({\mbox{\boldmath$0$}}_{d_{0}},{\mathbf{I}}_{d_{0}})$ rows that are independent of ${\mathbf{X}}$ and ${\mathbf{e}}$. We follow Kobak et al. (2020) in setting the responses corresponding to ${\mathbf{W}}$ equal to zero (the marginal mean response); a more complicated alternative would be to set them equal to $N_{n-n_{0}}({\mbox{\boldmath$0$}}_{n-n_{0}},\sigma^{2}{\mathbf{I}}_{n-n_{0}})$ noise that is independent of ${\mathbf{e}}$, ${\mathbf{X}}$ and ${\mathbf{W}}$.

As we are interested in the conflicting consequences of the two types of sequences, we study them separately rather than together (i.e. both $d,n\rightarrow\infty$) as in Bartlett et al. (2020) and Hastie et al. (2022). These authors make all the predictors important ($d=d_{0}$), although Hastie et al. (2022) introduced noise by adding measurement error to the predictors. Nakkiran et al. (2021) increased the number of predictors and the number of observations separately, but made these all important predictors and true observations. Intuitively, insofar as many predictors are “less important” (have non-zero but small coefficients), they have a similar effect to noise predictors. However, clearly separating important and noise predictors as we do is essential for understanding their separate effects. Mitra (2019), albeit studying different estimators from us, included noise predictors but mixed them with important predictors and masked some of the effects we observe. Kobak et al. (2020) also considered both sequences I and II but used specially normalized noise, replacing ${\mathbf{Z}}$ by $(d-d_{0})^{-1/2}\lambda_{z}^{1/2}{\mathbf{Z}}$ and ${\mathbf{W}}$ by $(n-n_{0})^{-1/2}\lambda_{w}^{1/2}{\mathbf{W}}$, where $\lambda_{z},\lambda_{w}>0$. We argue that our results for unnormalized noise predictors are more interesting and relevant because normalized noise predictors have to be constructed and added by the analyst whereas we allow noise predictors to simply occur in the data and/or to be added by the analyst. The noise predictors of Kobak et al. (2020) are also asymptotically zero, so increasingly negligible, and following the general recommendation to standardize the predictors (to have variance one) to reduce heterogeneity would remove the normalization. Muthukumar et al. (2020) considered sequences similar to our sequence I, but in their empirical work they only considered one of the cases we consider, namely $d_{0}<<n$. In this case, we need to include a large number of noise predictors to reach the interpolation point $d=n$, a circumstance that hides some findings we make by considering a wider set of relationships between $d_{0}$ and $n$. Finally, Hellkvist et al. (2023) also considered our sequence I (calling noise predictors fake features) when the unknown parameters are zero-mean random vectors. They did not consider noise observations.

In our empirical and theoretical study of the estimators and their test errors, we establish two important sets of results.

1) Adding noise predictors or observations to the data as in sequences I and II shrinks the estimators to zero as $d\rightarrow\infty$ or $n\rightarrow\infty$. We compute the test errors for the estimators under the two sequences and prove that they both asymptote to the same value as $d\rightarrow 0$ or $\infty$ and $n\rightarrow 0$ or $\infty$. This is a key result for understanding double descent: the convergence of test errors to the same value at both zero and infinity, combined with the ‘condition number anomaly’ (Dax 2022), which is reflected in the test error becoming infinite in Theorems 3.1 and 3.2, collectively drive the double descent phenomenon.

Although there are results relevant to sequence I in the literature (e.g. similar expressions for test errors in Hellkvist et al. (2023)), sequence II has not been considered and our results for sequence II are new. Even for sequence I, our results are not as obvious as they may appear. For example, the estimators do not shrink to zero along the sequences of Kobak et al. (2020) or in the setting without noise predictors of Belkin et al. (2020) who suggest that using as many predictors as possible may be optimal in their setting. For a particular case of sequence I and with a different focus on alias and spuriously correlated predictors, Muthukumar et al. (2020) discussed ‘signal bleeding’ into a large number of alias predictors – which is over-shrinkage – and overfitting of noise under parsimonious selection of noise predictors in scenarios where sparse or limited data makes it challenging to distinguish between true predictors and those that are spuriously correlated with the outcome in the training set. For the test errors when increasing the number of predictors, both Muthukumar et al. (2020) and Hastie et al. (2022) found as we do that the test error of the model approaches that of the null model. We show empirically and theoretically that this holds in greater generality than considered by Muthukumar et al. (2020), when we add noise predictors (rather than predictors with measurement error as in Hastie et al. (2022)), and for sequence II. We suggest that in a limited training set, the less-important predictors in Bartlett et al. (2020) behave similarly to noise, shrinking coefficient estimates toward zero. Importantly, this implies that overfitting noise is not always harmless or benign (Muthukumar et al., 2020; Bartlett et al., 2020; Tsigler & Bartlett, 2023). Finally, for well-specified models, Nakkiran (2019) related the peak in the test error to the condition number of the data matrix. In addition, we use the “condition number anomaly” (Dax, 2022) to obtain additional theoretical insights into why the peak always occurs. This result combined with the fixed asymptotes of the test error establishes that double descent has to occur in our setup (and incidentally, also in that of Hastie et al. (2022)). Showing this for both sequence I and II is important as it highlights the importance of shrinkage induced by noise rather than the specific form of the noise in driving double descent.

2) For both sequences I and II, the estimators are approximately ridge regression estimators. Moreover, for sequence II, the ridge regression estimator applied to the augmented data is also approximately a “double shrunk” ridge regression estimator. The “double shrinking” means that the optimal value of the ridge parameter in this estimator can be negative when $n\rightarrow\infty$, even when the predictors are independent and the number of predictors $d=d_{0}$ is small (i.e. a low-dimensional setting).

The interpretation of least squares estimators applied to noisy data as being like ridge regression estimators applied to “pure” data originated in the measurement error setting (Webb, 1994; Bishop, 1995; Hastie et al., 2022). This is different from adding additional noise predictors or observations. Kobak et al. (2020) added additional normalized noise predictors, but ignored the estimated coefficients of the noise predictors, and noted that a similar interpretation holds when we add normalized noise observations to the data. Our results for the full regression parameter estimator avoid the need to normalize the noise predictors and observations differently from the true versions. The seminal paper on ridge regression (Hoerl & Kennard, 1970) showed that for $d<n$, a ridge estimator with a positive ridge parameter $\lambda$ has a lower mean squared error than the OLS estimator, provided the predictor matrix is of full rank $d$. Subsequent research by Dobriban & Wager (2018) and Hastie et al. (2022) in the high dimensional setting showed that, for predictors with any covariance matrix $\boldsymbol{\Sigma}$, as $d,n\to\infty$ with $d/n=\gamma$, the optimal ridge parameter $\lambda_{\text{opt}}$ is positive when the orientation of $\boldsymbol{\beta}$ is random. Nakkiran et al. (2020) showed that a positive $\lambda$ is optimal regardless of the rank of the predictor matrix or the orientation of $\boldsymbol{\beta}$, as long as the predictors are independent. Recently, Liang & Rakhlin (2020) and Kobak et al. (2020) have shown that $\lambda_{opt}$ can be zero or negative when the predictors are correlated. Furthermore, a recent non-asymptotic analysis by Tsigler & Bartlett (2023) showed that under a specific ‘spiked’ covariance structure for the data in which $\boldsymbol{\Sigma}$ has a few large eigenvalues, a negative $\lambda$ can improve generalization bounds. We find that $\lambda_{opt}$ is positive under sequence I (which is consistent with Nakkiran et al. (2020)) but can be negative under sequence II, even when the predictors are independent and even in the low dimensional setting. This is a completely new result.

Our results under sequence I show double descent occurring in the overparameterized regime, but our results for sequence II show it occurring in the underparameterized regime. This shows that double descent is not driven directly by overparameterization, but rather is driven by shrinkage which can be induced by overparameterization due to including irrelevant noise predictors or by including irrelevant noise observations (or by measurement error in the predictors). The emphasis should be placed on choosing the correct shrinkage, rather than on fitting the exactly correct (possibly sparse) model ($d=d_{0}$) or an overparameterized (usually complex) model ($d>>d_{0}$).

The remainder of the paper is organised as follows. We present numerical results to illustrate our key findings in Section 2 and then present theoretical results that explain our numerical results in Section 3. We present a data application in Section 4 before concluding with a discussion in Section 5. Additional material and proofs of our theorems are given in Supplementary Material.

We construct the two sequences of training data $\mathcal{D}^{(d)}$ and $\mathcal{D}_{(n)}$ as described in the Introduction. Then we fit working linear regression models to the two sequences of training data by estimating the regression parameters using the ordinary least squares (OLS) and minimum norm OLS estimators when $d<n_{0}$ and $d>n_{0}$ for sequence I, and $n>d_{0}$ and $n<d_{0}$ for sequence II, respectively. We write these linear estimators as $\hat{{\mbox{\boldmath$\beta$}}}^{(d)}={\mathbf{M}}^{(d)}{\mathbf{y}}$ for sequence I and $\hat{{\mbox{\boldmath$\beta$}}}_{(n)}={\mathbf{M}}_{(n)}{\mathbf{y}}$ for sequence II, where ${\mathbf{M}}^{(d)}$ is a $d\times n_{0}$ matrix and ${\mathbf{M}}_{(n)}$ is a $d_{0}\times n_{0}$ matrix. The matrix ${\mathbf{M}}^{(d)}$ depends on ${\mathbf{X}}^{(d)}$ or $[{\mathbf{X}},{\mathbf{Z}}]$ and ${\mathbf{M}}_{(n)}$ depends on ${\mathbf{X}}_{(n)}$ or $[{\mathbf{X}}^{T},{\mathbf{W}}^{T}]^{T}$, but neither depends on ${\mathbf{y}}$. For sequence I, when $d\leq\min(n_{0},d_{0})$, ${\mathbf{M}}^{(d)}=({\mathbf{X}}^{(d)T}{\mathbf{X}}^{(d)})^{-1}{\mathbf{X}}^{(d)T}$. If $d_{0}>n_{0}$, we switch for $n_{0}<d\leq d_{0}$ to the minimum norm OLS estimator with ${\mathbf{M}}^{(d)}={\mathbf{X}}^{(d)T}({\mathbf{X}}^{(d)}{\mathbf{X}}^{(d)T})^{-1}$ and then for $d>d_{0}$ ${\mathbf{M}}^{(d)}=[{\mathbf{X}},{\mathbf{Z}}]^{T}({\mathbf{X}}{\mathbf{X}}^{T}+{\mathbf{Z}}{\mathbf{Z}}^{T})^{-1}$; if $d_{0}<n_{0}$, for $d_{0}<d\leq n_{0}$, we use ${\mathbf{M}}^{(d)}=\{[\mathbf{X},\mathbf{Z}]^{T}[\mathbf{X},\mathbf{Z}]\}^{-1}[\mathbf{X},\mathbf{Z}]^{T}$ before we switch to the minimum norm OLS estimator for $d\geq n_{0}$. For sequence II, when $n\leq\min(n_{0},d_{0})$, ${\mathbf{M}}_{(n)}={\mathbf{X}}_{(n)}^{T}({\mathbf{X}}_{(n)}{\mathbf{X}}_{(n)}^{T})^{-1}[{\mathbf{I}}_{n},{\mbox{\boldmath$0$}}_{n,n_{0}-n}]$. If $d_{0}<n_{0}$, we switch to the OLS estimator for $d_{0}<n\leq n_{0}$ so ${\mathbf{M}}_{(n)}=({\mathbf{X}}_{(n)}^{T}{\mathbf{X}}_{(n)})^{-1}{\mathbf{X}}_{(n)}^{T}[{\mathbf{I}}_{n},{\mbox{\boldmath$0$}}_{n,n_{0}-n}]$ and then for $n>n_{0}$ ${\mathbf{M}}_{(n)}=({\mathbf{X}}^{T}{\mathbf{X}}+{\mathbf{W}}^{T}{\mathbf{W}})^{-1}{\mathbf{X}}^{T}$; if $d_{0}>n_{0}$, for $n_{0}<n\leq d_{0}$ ${\mathbf{M}}_{(n)}={\mathbf{X}}^{T}{\mathbf{A}}^{11}+{\mathbf{W}}^{T}{\mathbf{A}}^{21}$, where ${\mathbf{A}}^{11}=\{{\mathbf{X}}{\mathbf{X}}^{T}-{\mathbf{X}}{\mathbf{W}}^{T}({\mathbf{W}}{\mathbf{W}}^{T})^{-1}{\mathbf{W}}{\mathbf{X}}^{T}\}^{-1}$ and ${\mathbf{A}}^{21}=-({\mathbf{W}}{\mathbf{W}}^{T})^{-1}{\mathbf{W}}{\mathbf{X}}^{T}{\mathbf{A}}^{11}$, before we switch to the OLS estimator for $n>d_{0}$. The second expression for ${\mathbf{M}}_{(n)}$ is obtained by letting ${\mathbf{A}}=[{\mathbf{X}}^{T},{\mathbf{W}}^{T}]^{T}[{\mathbf{X}}^{T},{\mathbf{W}}^{T}]$, partitioning it into its natural four submatrices and then using standard formulas for the inverse of a partitioned matrix to obtain the corresponding terms in ${\mathbf{A}}^{-1}$.

We showed in Figures 3 and 5 that adding noise predictors in sequence I or noise observations in sequence II induces shrinkage in $\hat{{\mbox{\boldmath$\beta$}}}^{(d)}$ and $\hat{{\mbox{\boldmath$\beta$}}}_{(n)}$, respectively, and eventually shrinks these estimators to zero. We now confirm this effect theoretically as $d$ or $n\rightarrow\infty$. We let ‘a.s.’ mean ‘almost surely’ and ‘P’ mean ’in probability’. Also, let $\|{\mathbf{x}}\|=({\mathbf{x}}^{T}{\mathbf{x}})^{1/2}$ denote the Euclidean norm of ${\mathbf{x}}$. For sequence I, for $d>\max(d_{0},n_{0})$, write $\hat{\boldsymbol{\beta}}^{(d)}=(\hat{\boldsymbol{\beta}}^{(d)T}_{d_{0}},\hat{\boldsymbol{\beta}}^{(d)T}_{-d_{0}})^{T}$, where $\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}=\mathbf{X}^{T}(\mathbf{X}\mathbf{X}^{T}+{\mathbf{Z}}{\mathbf{Z}}^{T})^{-1}\mathbf{y}$ denotes the first $d_{0}$ elements and $\hat{\boldsymbol{\beta}}^{(d)}_{-d_{0}}=\mathbf{Z}^{T}(\mathbf{X}\mathbf{X}^{T}+{\mathbf{Z}}{\mathbf{Z}}^{T})^{-1}\mathbf{y}$ denotes the remaining $d-d_{0}$ elements of $\hat{\boldsymbol{\beta}}^{(d)}$. By the strong law of large numbers, the $n_{0}\times n_{0}$ matrix ${\mathbf{X}}{\mathbf{X}}^{T}+{\mathbf{Z}}{\mathbf{Z}}^{T}=O(d)$ a.s. as $d\rightarrow\infty$. Since $\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}\in\mathbb{R}^{d_{0}}$ and $\hat{\boldsymbol{\beta}}^{(d)}_{-d_{0}}\in\mathbb{R}^{d-d_{0}}$, we have $\|\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}\|=O(d^{-1}d_{0}^{1/2})$ and $\|\hat{\boldsymbol{\beta}}^{(d)}_{-d_{0}}\|=O(d^{-1}(d-d_{0})^{1/2})$ a.s. as $d\rightarrow\infty$. Similarly, for sequence II, for $n>\max(n_{0},d_{0})$, $\hat{\boldsymbol{\beta}}_{(n)}=(\mathbf{X}^{T}\mathbf{X}+{\mathbf{W}}^{T}{\mathbf{W}})^{-1}\mathbf{X}^{T}\mathbf{y}$, the $d_{0}\times d_{0}$ matrix ${\mathbf{X}}^{T}{\mathbf{X}}+{\mathbf{W}}^{T}{\mathbf{W}}=O(n)$ a.s. as $n\rightarrow\infty$, and $\|\hat{{\mbox{\boldmath$\beta$}}}_{(n)}\|=O(n^{-1}n_{0})$ a.s. as $n\rightarrow\infty$. That is, for $d_{0}$ and $n_{0}$ fixed or increasing sufficiently slowly in sequence I and II, respectively, noise predictors or observations eventually shrink the estimated coefficients to zero.

We examine the performance of the working models fitted along the two training sequences through the test error. In our empirical work, we compute empirical test errors by computing the mean squared error over a large test sample. In theoretical work, we calculate the model-expectation of the squared error at a single observation $(y^{\prime},{\mathbf{x}}^{{}^{\prime}(d)})\in\mathbb{R}^{d+1}$ that is generated independently from the same data generating model as the training data, namely ${\mathbf{x}}^{\prime}\sim N_{d_{0}}({\mbox{\boldmath$0$}}_{d_{0}},{\mathbf{I}}_{d_{0}})$, $y^{\prime}={\mathbf{x}}^{{}^{\prime}T}\boldsymbol{\beta}_{0}+e^{\prime}$ with $e^{\prime}\sim N(0,\sigma^{2})$, and ${\mathbf{x}}^{{}^{\prime}(d)}$ contains the first $d$ elements of ${\mathbf{x}}^{\prime}$ if $d\leq d_{0}$, and ${\mathbf{x}}^{{}^{\prime}(d)}=[{\mathbf{x}}^{{}^{\prime}T},{\mathbf{z}}^{{}^{\prime}T}]$, where ${\mathbf{z}}^{\prime}\sim N_{d-d_{0}}({\mbox{\boldmath$0$}}_{d-d_{0}},{\mathbf{I}}_{d-d_{0}})$ independent of ${\mathbf{x}}^{\prime}$, if $d>d_{0}$. We obtain different test errors, depending on what we condition on. For example, corresponding to the light red curves in our figures, we have the conditional (on the training data) test errors

Note that for sequence II, we compute the test error under the model for the real data $(y^{\prime},\mathbf{x}^{\prime})$ rather than for noise observations $(0,\mathbf{w}^{\prime})$ because this is the prediction we are actually interested in; alternatives such as the test error under the noise observations or a linear combination of these two possibilities (e.g. with weights $n_{0}/n$ and $(n-n_{0})/n$) could also be considered. The dark red curves in our figures are (unconditional) test errors $R^{(d)}$ and $R_{(n)}$ obtained by taking the expectations (over $\mathcal{D}^{(d)}$ or $\mathcal{D}_{(n)}$) of $R^{(d)}(\mathcal{D}^{(d)})$ and $R_{(n)}(\mathcal{D}_{(n)})$.

The test error under sequence I has been obtained by Breiman & Freedman (1983) and Belkin et al. (2019). We state the theorem for completeness.

The test error under sequence II is more complicated to work with than that under sequence I because analytic expressions are not available once we augment the data with noise observations ($n>n_{0}$). Nonetheless, we obtain the following theorem.

We saw empirically in Figures 3 and 5 and theoretically in Section 3.1 that adding noise to the training data shrinks the estimators to zero as $d$ or $n\rightarrow\infty$. This means that, as shown in Figure 1, the test errors should converge as $d$ or $n\rightarrow\infty$ to the test error of the null model with $d=n=0$. That is, both ends of the test error should equal the same value $\operatorname{E}(y^{{}^{\prime}2})=\|{\mbox{\boldmath$\beta$}}_{0}\|^{2}+\sigma^{2}$. This has also been observed by Muthukumar et al. (2020) and, in a different setup, by Hastie et al. (2022) but it has not been proved to occur. That these asymptotes hold under sequence I follows immediately from Theorem 3.1. It is much more challenging to prove that they hold under sequence II (because there is no analytic form for the test error under sequence II). We now prove that this occurs in Theorem 3.3; the proof is given in the Supplementary Material.

The fact that both ends of the test error equal the same value means double descent occurs whenever the test error increases above the initial test error. The test error can increase (e.g. Figure 1c) or decrease (e.g. Figure 1a) from its value at the origin. If the test error increases immediately, it has a local minimum at the origin (which we still call the first descent), while if it decreases and then increases above its initial value, it has a local minimum before the interpolation point ($d=n$). In either case, it will have to descend again (the second descent) to achieve its asymptote. If it decreases below the null model test error, the test error has to increase again to achieve its asymptote. The second descent does not necessarily achieve a global minimum; whether this occurs or not is related to the amount of shrinkage induced by the noise variables or observations.

What makes the test error increase above the null model test error? The estimators and their test errors involve $({\mathbf{D}}^{T}{\mathbf{D}})^{-1}$ or $({\mathbf{D}}{\mathbf{D}}^{T})^{-1}$ with ${\mathbf{D}}$ equal to ${\mathbf{X}}^{(d)}$, $[{\mathbf{X}},{\mathbf{Z}}]$, ${\mathbf{X}}_{(n)}$ or $[{\mathbf{X}}^{T},{\mathbf{W}}^{T}]^{T}$. The underlying driver behind the test error behavior is the generalized condition number $\kappa(\mathbf{D})$ of these matrices. In our examples, the elements of the matrix $\mathbf{D}\in\mathbb{R}^{n\times d}$ are independently sampled from the same probability distribution. The generalized condition number $\kappa(\mathbf{D})$ of a $n\times d$ rectangular matrix $\mathbf{D}$ of rank $r\leq\mbox{min}(n,d)$ is $\kappa(\mathbf{D})=\gamma_{1}(\mathbf{D})/\gamma_{r}(\mathbf{D})$, where $\gamma_{1}(\mathbf{D})\geq\ldots\geq\gamma_{r}(\mathbf{D})$ are the non-zero singular values of $\mathbf{D}$ in descending order (Zielke, 1988). As a direct consequence of the Cauchy Interlacing Theorem (Horn & Johnson, 2012, Theorem 4.3.4), for a fixed $n$ and $d<n$, the condition number of $\mathbf{D}$ increases monotonically as $d\to n$; if $\mathbf{D}$ is of full rank, it reaches its maximum value when $d=n$. The increasing condition number of $\mathbf{D}$ causes the variance and the bias of the estimated coefficients to increase, with the maximum variance and maximum bias and hence maximum test error at the interpolation point. If this maximum is above the asymptote, we observe a second descent in the overparameterized regime. For $d>n$, the condition number of $\mathbf{D}$ is itself likely to decrease as we increase $d$ (the“condition number anomaly”, Dax, 2022). This behavior is illustrated in Figure S2 in the Supplementary Material. The shape of the curves is strikingly similar to that of the variance curves in the second row of Figure 2.

We can obtain additional insight into the estimators by approximating them by the ridge regression estimator $\hat{\boldsymbol{\beta}}_{\lambda}=(\mathbf{X}^{T}\mathbf{X}+\lambda\mathbf{I}_{d_{0}})^{-1}\mathbf{X}^{T}\mathbf{y}={\mathbf{X}}^{T}(\mathbf{X}\mathbf{X}^{T}+\lambda\mathbf{I}_{n_{0}})^{-1}{\mathbf{y}}$. For sequence I, write $\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}-\hat{\boldsymbol{\beta}}_{\lambda}=\mathbf{X}^{T}{\mathbf{L}}^{(d)}(\lambda)\mathbf{y}$, where ${\mathbf{L}}^{(d)}(\lambda)=(\mathbf{X}\mathbf{X}^{T}+{\mathbf{Z}}{\mathbf{Z}}^{T})^{-1}(\lambda\mathbf{I}_{n_{0}}-{\mathbf{Z}}{\mathbf{Z}}^{T})(\mathbf{X}\mathbf{X}^{T}+\lambda\mathbf{I}_{n_{0}})^{-1}$ is an $n_{0}\times n_{0}$ matrix. Let $d^{-1}d_{0}\rightarrow\delta$, with $0\leq\delta<1$, as $d\rightarrow\infty$. Then $(d-d_{0})^{-1}{\mathbf{Z}}{\mathbf{Z}}^{T}-{\mathbf{I}}_{n_{0}}=O_{p}((d-d_{0})^{-1/2})$ as $d\rightarrow\infty$ and hence, ${\mathbf{L}}^{(d)}(d-d_{0})=O_{p}(d^{-3/2})$, so $\|\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}-\hat{\boldsymbol{\beta}}_{d-d_{0}}\|=O_{p}(d^{-1})$ and, from above, $\|\hat{\boldsymbol{\beta}}^{(d)}_{-d_{0}}\|=O(d^{-1/2})$ a.s., so all the elements converge to zero. For sequence II, let $n^{-1}n_{0}\rightarrow\nu$ with $0\leq\nu<1$. Then write $\hat{\boldsymbol{\beta}}_{(n)}-\hat{\boldsymbol{\beta}}_{(n_{0}),n-n_{0}}={\mathbf{L}}_{(n)}(\lambda){\mathbf{X}}^{T}{\mathbf{y}}$, where the $d_{0}\times d_{0}$ matrix ${\mathbf{L}}_{(n)}(\lambda)=(\mathbf{X}^{T}\mathbf{X}+{\mathbf{W}}^{T}{\mathbf{W}})^{-1}\{\lambda\mathbf{I}_{d_{0}}-{\mathbf{W}}^{T}{\mathbf{W}}\}\{{\mathbf{X}}^{T}{\mathbf{X}}+\lambda{\mathbf{I}}_{d_{0}}\}^{-1}$. We have ${\mathbf{L}}_{(n)}(n-n_{0})=O_{p}(n^{-3/2})$ and $\|\hat{\boldsymbol{\beta}}_{(n)}-\hat{\boldsymbol{\beta}}_{n-n_{0}}\|=O_{p}(n^{-1/2})$. We gather the results into a Theorem.

Related results were given by Kobak et al. (2020) for normalized noise variables. Specifically, they normalized the noise predictors and replaced ${\mathbf{Z}}$ in the training set $\mathcal{D}^{(d)}$ by $\tilde{{\mathbf{Z}}}=(d-d_{0})^{-1/2}\lambda_{z}^{1/2}{\mathbf{Z}}$, and replaced ${\mathbf{W}}$ in the training set $\mathcal{D}_{(n)}$ by $\tilde{{\mathbf{W}}}=(n-n_{0})^{-1/2}\lambda_{w}^{1/2}{\mathbf{W}}$. In the first case, ${\mathbf{L}}^{(d)}=O_{p}((d-d_{0})^{-1/2})$, so $\|\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}-\hat{\boldsymbol{\beta}}_{\lambda_{z}}\|=O_{p}((d-d_{0})^{-1/2})$; we can also easily obtain the result of Kobak et al. (2020) that $\|\hat{\boldsymbol{\beta}}^{(d)}_{d_{0}}-\hat{\boldsymbol{\beta}}_{\lambda}\|=o(1)$ a.s., although the rates of convergence give additional information. For the remaining elements, $\|\hat{\boldsymbol{\beta}}^{(d)}_{-d_{0}}\|=O_{p}(1)$ which is not useful. In the second case, ${\mathbf{L}}_{(n)}=O_{p}(n^{-1/2})$ and $\|\hat{\boldsymbol{\beta}}_{(n)}-\hat{\boldsymbol{\beta}}_{\lambda}\|=O_{p}(n^{-1/2})$. As explained in the Introduction, we argue that our results are more interesting and useful than these.

In Figure 6, we showed that the optimal value of the ridge penalty in a ridge regression estimator can be negative. Intuitively, this makes sense because including noise in the test set has the same effect as ridge regularisation so, if we have already overshrunk through including too much noise, the optimal action may be to reduce the shrinkage by using a negative ridge parameter. Theorem 3.5 below which extends Theorem 3.4 for sequence II supports this intuition.

Suppose that for sequence II, we set $n_{0}=n\nu$ with $0<\nu<1$. In Figure 6, this corresponds to making an approximation along the line $n-n_{0}=n_{0}\nu^{-1}(1-\nu)$ which passes through the origin and has slope $\nu^{-1}(1-\nu)$. Let $\hat{{\mbox{\boldmath$\beta$}}}_{(n),\lambda}$ denote the ridge regression estimator computed along sequence II. Then write $\hat{\boldsymbol{\beta}}_{(n),n\lambda}-\hat{\boldsymbol{\beta}}_{n(1-\nu+\lambda)}={\mathbf{L}}_{(n)}{\mathbf{X}}^{T}{\mathbf{y}}$, where the $d_{0}\times d_{0}$ matrix ${\mathbf{L}}_{(n)}=(\mathbf{X}^{T}\mathbf{X}+{\mathbf{W}}^{T}{\mathbf{W}}+n\lambda{\mathbf{I}}_{d_{0}})^{-1}\{(n-n_{0})\mathbf{I}_{d_{0}}-{\mathbf{W}}^{T}{\mathbf{W}}\}\{{\mathbf{X}}^{T}{\mathbf{X}}+n(1-\nu+\lambda){\mathbf{I}}_{d_{0}}\}^{-1}$. We have ${\mathbf{L}}_{(n)}=O_{p}(n^{-3/2})$ and $\|\hat{\boldsymbol{\beta}}_{(n),n\lambda}-\hat{\boldsymbol{\beta}}_{n(1-\nu+\lambda)}\|=O_{p}(n^{-1/2})$, which proves the following Theorem.

There is an analogous result for the normalized predictors used by Kobak et al. (2020) instead of sequence II. Recall that Kobak et al. (2020) keep $n_{0}$ fixed and replace ${\mathbf{W}}$ in the training set $\mathcal{D}_{(n)}$ by $\tilde{{\mathbf{W}}}=(n-n_{0})^{-1/2}\lambda_{w}^{1/2}{\mathbf{W}}$. Then we can show that $\|\hat{\boldsymbol{\beta}}_{(n),\lambda}-\hat{\boldsymbol{\beta}}_{\lambda+\lambda_{w}}\|=O_{p}(n^{-1/2})$. The optimal ridge parameter in the approximating ridge regression estimator is $\lambda_{opt}=\|{\mbox{\boldmath$\beta$}}_{0}\|^{-2}d_{0}\sigma^{2}-\lambda_{w}$ which can also be zero or negative, but does not explicitly depend on $n$. Figure 6 shows empirically that the optimal $\lambda$ should depend on $n$, as shown in our theorem.

Theorem 3.5 does not hold for sequence I because the noise predictors are independent (Nakkiran et al. (2020)).