Surprises in High-Dimensional Ridgeless Least Squares Interpolation

Introduction

Modern deep learning models involve a huge number of parameters. In many applications, current practice suggests that we should design the network to be sufficiently complex so that the model (as trained, typically, by gradient descent) interpolates the data, i.e., achieves zero training error. Indeed, in a thought-provoking experiment, Zhang et al. (2016) showed that state-of-the-art deep neural network architectures are complex enough that they can be trained to interpolate the data even when the actual labels are replaced by entirely random ones.

Despite their enormous complexity, deep neural networks are frequently observed to generalize well in practice. At first sight, this seems to defy conventional statistical wisdom: interpolation (vanishing training error) is commonly taken to be a proxy for overfitting, poor generalization (large gap between training and test error), and hence large test error. In an insightful series of papers, Belkin et al. (2018b, c, a) pointed out that these concepts are in general distinct, and interpolation does not contradict generalization. For example, recent work has investigated interpolation —via kernel ridge regression— in reproducing kernel Hilbert spaces (Liang et al., 2020; Ghorbani et al., 2019). While in low dimension a positive regularization is needed to achieve good interpolation, in certain high dimensional settings interpolation can be nearly optimal.

We study the model (2) in the proportional regime $p\asymp n$, with a special focus on the overparametrized case $p>n$. Our main contribution is to show that, by considering different choices of the features distribution $P_{x}$, we can reproduce a number of statistically interesting phenomena that have emerged in the context of deep learning.

From a technical perspective, our main results are: Theorems 2 and 5, which assume the linear model $x_{i}=\Sigma^{1/2}z_{i}$ with $z_{i}$ a vector with independent coordinates; and Theorem 8, which assumes a nonlinear model $x_{i}=\varphi(Wz_{i})$ with $z_{i}\sim N(0,I_{d})$. While the linear model has already a attracted significant amount of work (see Section 1.3 for an overview), Theorems 2 and 5 provide a more accurate non-asymptotic approximation of the prediction risk, as compared to available results in the literature.

The prediction risk depends on the geometry of the pair $(\Sigma,\beta)$. We consider a few different choices for this geometry, which are broadly motivated by our objective to understand overparametrized models, and specialize our formulas to these special cases:

Isotropic features. This is the simplest case, in which $\Sigma=I_{p}$ and therefore —as we will see— the asymptotic risk depends on $\beta$ only through its norm $\|\beta\|_{2}$. This simple model captures some interesting features of overparametrization, but misses others.

Nonlinear model. In all of the previous cases, the distribution of $x_{i}$ is of the form $x_{i}=\Sigma^{1/2}z_{i}$ where $z_{i}$ is a vector with independent coordinates. In order to test the generality of our results, we consider a model in which $x_{i}$ is obtained by passing $z_{i}\sim{N}(0,I_{d})$ through a one-layer neural net with random first layer weights.

We will summarize our results for these four examples in the next subsection.

A skeptical reader might ask what linear models have to do with neural networks. We emphasize that linear models provide more than a simple analogy, and a recent line of work outlines a concrete connection between the two settings (Jacot et al., 2018; Du et al., 2018b, a; Allen-Zhu et al., 2018; Chizat and Bach, 2018b). We will discuss this connection in Section 1.2.

We denote by $\gamma:=p/n\in(0,\infty)$ the overparametrization ratio. When $\gamma<1$, we call the problem underparametrized, and when $\gamma>1$, we call it overparametrized. Our most general results for the linear model (Theorem 2 and 5) apply to a non-asymptotic setting in which $n,p$ are finite, and provide a deterministic approximation of the risk with error bounds that are uniform in the distribution of the data. We will occasionally consider a simpler asymptotic scenario in which both $p$ and $n$ diverge with $p/n\to\gamma$.

Our main results are twofold: $(i)$ We show that by taking $\Sigma\neq I_{p}$, we can easily construct scenarios in which the minimum of the risk is achieved in the overparamertized regime $p>n$; $(ii)$ We show that these findings are robust to the details of the distribution of $(y_{i},x_{i})$.

As a preliminary remark, note that in the underparametrized regime ($\gamma<1$), the min-norm estimator coincides with the standard least squares estimator. Its risk is purely variance (there is no bias), and does not depend on $\beta,\Sigma$ (see Proposition 2). Interestingly, the asymptotic risk diverges as we approach the interpolation boundary (as $\gamma\to 1$).

In contrast, in the overparametrized regime ($\gamma>1$), the risk is composed of both bias and varianceNote that in the overparametrized regime the bias is non-vanishing even in the interpolation limit $\lambda\to 0$. The reason is that the set of interpolators is an affine space of dimension $p-n$, and the min-norm criterion selects one specific interpolator which has norm smaller than $\beta$., and generally depends on $\beta,\Sigma$ (see Theorem 2).

We next highlight some concrete results for the four models discussed in the previous section (unless explicitly said, we refer to the min-norm estimator).

In either case, the risk approaches the null risk as $\gamma\to\infty$, and achieves its global minimum in the underparametrized regime (see Section 3.2).

Optimal tuning of the ridge penalty can be achieved by leave-one-out cross-validation (see Theorem 7).

In this case $\Sigma\neq{{I}}$ and the risk depends on the geometry of $(\Sigma,\beta)$, and in particular on how $\beta$ aligns with the eigenvectors of $\Sigma$.

If the coefficients vector is equidistributed along the eigenvectors of $\Sigma$, the behavior is qualitatively similar to the isotropic case. This situation arises, for instance, if $\beta$ is itself random with a spherical prior.

If $\beta$ is aligned with the top eigenvectors of $\Sigma$, the situation is qualitatively different. As an example we obtain an explicit formula for the asymptotic risk in the latent space model discussed above, see Figure 5 for an illustration. We find that, for natural choices of the model parameters, the risk is monotone decreasing in the overparametrized regime, and reaches its global minimum as $\gamma\to\infty$. This qualitative behavior matches the one observed for neural networks (see Section 5.4).

For the latent space model, we observe that, at large overparametrization, the minimum error is achieved as $\lambda\to 0$, i.e. by min-norm interpolators (see Section 6.2, and Section 1.3 for related work).

From a technical viewpoint, analysis of the isotropic covariates model is straightforward and relies on standard random matrix theory results. However, we believe it provides useful insights.

In contrast, the results for general covariance and coefficients structure $(\Sigma,\beta)$ is technically novel. We discuss related work in Section 1.3. Our results for the nonlinear model are also novel. In this setting, we derive a new asymptotic result on resolvents of certain block matrices, which may be of independent interest (see Lemma 3).

We next discuss some intuition behind and implications of our results.

The shape of the asymptotic risk curve for min-norm least squares is, of course, controlled by its components: bias and variance. For fully specified models, the bias increases with $\gamma$ in the overparametrized regime, which is intuitive. When $p>n$, the min-norm least squares estimate of $\beta$ is constrained to lie the row space of $X$, the training feature matrix. This is a subspace of dimension $n$ lying in a feature space of dimension $p$. Thus as $p$ increases, so does the bias, since this row space accounts for less and less of the ambient $p$-dimensional feature space.

Recently, Belkin et al. (2018a) pointed out a fascinating empirical trend where, for popular methods like neural networks and random forests, we can see a second bias-variance tradeoff in the out-of-sample prediction risk beyond the interpolation limit. The risk curve here resembles a traditional U-shape curve before the interpolation limit, and then descends again beyond the interpolation limit, which these authors call “double descent”. A closely related phenomenon was found earlier by Spigler et al. (2018), who studied the “jamming transition” from underparametrized to overparametrized neural networks. Our results formally verify that this double descent phenomenon occurs even in the simple and fundamental case of least squares regression. The appearance of the second descent in the risk, past the interpolation boundary ($\gamma=1$), is explained by the fact that the variance decreases as $\gamma$ grows, as discussed above.

In the misspecified case, the variance still decreases with $\gamma$ (for the same reasons), but interestingly, the bias can now also decrease with $\gamma$, provided $\gamma$ is not too large (not too far past the interpolation boundary). The intuition here is that in a misspecified model, some part of the true regression function is always unaccounted for, and adding features generally improves our approximation capacity. As a consequence, the double descent phenomenon can be even more pronounced in the misspecified case (depending on the strength of the approximation bias), and that the risk can attain its global minimum past the interpolation limit.

The min-norm least squares estimator can be seen as the limit of ridge regression as the tuning parameter tends to zero. It is also the convergence point of gradient descent run on the least squares loss. We would not in general expect the best-predicting ridge solution to be at the end of its regularization path. Our results, comparing min-norm least squares to optimally-tuned ridge regression, show that (asymptotically) this is never the case, when $\beta$ is incoherent with respect to the eigenvectors of $\Sigma$. This is for instance the case when $\Sigma=I_{p}$, or $\beta$ is distributed according to a spherically symmetric prior. In contrast, Kobak et al. (2020) recently pointed out that —when $\beta$ is aligned with the leading eigenvectors of $\Sigma$— min-norm regression can have optimal risk (i.e. the optimal regularization vanishes). We show that this is indeed the case in the latent space model mentioned above: this provides indeed an extremely simple example of a phenomenon that has been observed in the past for kernel methods Liang et al. (2020).

As mentioned above, early-stopped gradient descent is known to be closely connected to ridge regularization, see, e.g., Ali et al. (2019) which proves a tight coupling between the two (see Section 1.3 for further related work).

In practice, of course, we would not have access to the optimal tuning parameter for ridge (optimal stopping for gradient descent), and we would rely on, e.g., cross-validation (CV). Our theory shows that for ridge regression, CV tuning is asymptotically equivalent to optimal tuning (and we would expect the same results to carry over to gradient descent, but have not pursued this formally).

2 Connection to neural networks

As mentioned above, recent literature has pointed out a direct connection between linear models and more complex models such as neural networks (Jacot et al., 2018; Du et al., 2018b, a; Allen-Zhu et al., 2018; Chizat and Bach, 2018b). Our work contributes to this line of work and in particular points at the important role played by universality, a core concept in random matrix theory Tao (2012).

This model is still nonlinear in the input $z$, but is linear in the parameters $\beta$. We are therefore led to consider a linear regression problem, with random features $x_{i}=\nabla_{\theta}f(z_{i};\theta_{0})$, $i=1,\ldots,n$, of high-dimensionality ($p$ much greater than $n$). Notice that the features are random because of the initialization $\theta_{0}$. Further, since $p>n$, many vectors $\beta$ give rise to a model that interpolates the data.

A different approach relies on a universality hypothesis. Recall that $\{(y_{i},z_{i})\}_{i\leq n}$ have joint distribution $P_{y,z}$, and $x_{i}=\nabla_{\theta}f(z_{i};\theta_{0})$ and assume these quantities to have zero meanWe further note that the assumption of zero mean is mainly for convenience, as a non-zero mean of the covariates can be effectively compensated by an intercept.. Universality posits that the asymptotic risk of ridge regression (or min-norm interpolation) does not change if we replace the joint distribution $(x_{i},z_{i})$ by a Gaussian with the same covariance.

While universality is not expected to hold for any distribution of the data $(y_{i},z_{i})$, and for any function $f$, it has been confirmed in a number of examples of interest (see Section 1.3). In fact, our analysis of the nonlinear model in Section 8 provides a concrete example in which universality can be confirmed rigorously.

Under the universality hypothesis, we can assume $(z_{i},x_{i})$ to be jointly Gaussian. If the true joint distribution of the data $(y_{i},z_{i})$ is also Gaussian, it is sufficient to study the data distribution

where $\Sigma,\sigma$ are chosen to match the second order statistics of the original model. This is a special case of our ‘linear model’ (2). Summarizing, replacing the original nonlinear model $y_{i}=f(z_{i};\theta)$, with the linear model rests on two key steps. The first step is the linearization in Eq. (3): this has been proved to be a good approximation under certain ‘lazy training’ schemes Chizat and Bach (2018b). The second step is the universality hypothesis, which is quite common (albeit for simpler models) in random matrix theory Tao (2012).

Finally, although we believe our our results are relevant to understanding overparametrized neural networks, the concrete correspondence outlined above only holds in a certain ‘lazy training’ regime, in which network weights do not change much during training, More generally, in a neural network, the feature representation and the regression function or classifier are learned simultaneously. In both our linear and nonlinear model settings, the features $X$ are not learned, but observed. Learning $X$ could significantly change some aspects of the behavior of an interpolator. (See for instance Chapter 9 of Goodfellow et al. (2016), and also Chizat and Bach (2018b); Zhang et al. (2019), which emphasize the importance of learning the representation.)

3 Related work

Ridge regression with random features has been studied in the past. Dicker (2016) considers a model in which the covariates are isotropic Gaussian $x_{i}\sim N(0,I_{p})$ and computes the asymptotic risk of ridge regression in the proportional asymptotics $p,n\to\infty$, with $p/n\to\gamma\in(0,\infty)$. Dobriban and Wager (2018) generalize these results to $x_{i}=\Sigma^{1/2}z_{i}$, where $z_{i}$ has independent entries with bounded $12$-th moment.

Recently, Advani and Saxe (2017) study the effect of early stopping and ridge regularization in a model with isotropic Gaussian covariates $x_{i}\sim N(0,I_{p})$, again focusing on the proportional asymptotics $p,n\to\infty$, with $p/n\to\gamma\in(0,\infty)$. They show that this simple model reproduces several phenomena observed in neural networks training. The same model is reconsidered in concurrent work by Belkin et al. (2019), who obtain exact results for the expected risk of min-norm regression, relying on the jointly Gaussian distribution of $(y_{i},x_{i})$. We contribute to this line of work by extending the analysis to general covariance structures and to misspecified models. As we will see, these generalizations allow to produce examples for which the global minimum of the risk is achieved in the overparametrized regime $\gamma>1$.

The importance of the relation between the coefficient vector $\beta$ and the eigenvectors of $\Sigma$ was emphasized by Kobak et al. (2020) and Bartlett et al. (2020). These papers point out —under different asymptotic settings— that $\lambda=0+$ (i.e., min-norm regression) can be optimal or nearly optimal. After a preprint of this paper appeared, Wu and Xu (2020) and Richards et al. (2020) generalized our earlier results to cover the case in which $\beta$ is potentially aligned with $\Sigma$. We review in further detail these important generalizations in Section 4. We contribute to this line of work by obtaining non-asymptotic approximations for the risk, with explicit and nearly optimal error bounds. These hold under weaker assumptions on the geometry of $(\Sigma,\beta)$ than the results of Wu and Xu (2020); Richards et al. (2020).

For the nonlinear model, the random matrix theory literature is much sparser, and focuses on the related model of kernel random matrices, namely, symmetric matrices of the form $K_{ij}=\varphi(z_{i}^{T}z_{j})$. El Karoui (2010) studied the spectrum of such matrices in a regime in which $\varphi$ can be approximated by a linear function (for $i\neq j$) and hence the spectrum converges to a rescaled Marchenko-Pastur law. This approximation does not hold for the regime of interest here, which was studied instead by Cheng and Singer (2013) (who determined the limiting spectral distribution) and Fan and Montanari (2015) (who characterized the extreme eigenvalues). The resulting eigenvalue distribution is the free convolution of a semicircle law and a Marchenko-Pastur law. In the current paper, we must consider asymmetric (rectangular) matrices $x_{ij}=\varphi(w_{j}^{T}z_{i})$, whose singular value distribution was recently computed by Pennington and Worah (2017), using the moment method. Unfortunately, the prediction variance depends on both the singular values and vectors of this matrix. In order to address this issue, we apply the leave-one out method of Cheng and Singer (2013) to compute the asymptotics of the resolvent of a suitably extended matrix. We then extract the information of interest from this matrix. After appearance of a preprint of this paper, Mei and Montanari (2019) extended the results presented here, to obtain a complete characterization of the risk for the non-linear random features model.

Let us finally mention that the universality (or ‘invariance’) phenomenon is quite common in random matrix theory Tao (2012). In the context of kernel inner product random matrices, it appears (somewhat implicitly) in Cheng and Singer (2013) and (more explicitly) in Fan and Montanari (2015). After a first appearance of this manuscript, universality has been investigated in the context of neural networks in several papers Mei and Montanari (2019); Montanari et al. (2019); Gerace et al. (2020); Goldt et al. (2020); Hu and Lu (2020); Adlam and Pennington (2020).

4 Outline

Section 2 provides important background. Sections 3–7 consider the linear model case, focusing on isotropic features, correlated features, misspecified models, ridge regularization, and cross-validation, respectively. Section 8 covers the nonlinear model case. Nearly all proofs are deferred until the appendix.

Preliminaries

We describe our setup and gather a number of important preliminary results.

Consider a test point $x_{0}\sim P_{x}$, independent of the training data. For an estimator $\hat{\beta}$ (a function of the training data $X,y$), we define its out-of-sample prediction risk (or simply, risk) as

where $\|x\|_{\Sigma}^{2}=x^{T}\Sigma x$. Note that our definition of risk is conditional on $X$ (as emphasized by our notation $R_{X}$). Note also that we have the bias-variance decomposition

2 Ridgeless least squares

This can be equivalently written as $\hat{\beta}=(X^{T}X)^{+}X^{T}y$, where $(X^{T}X)^{+}$ is the pseudoinverse of $X^{T}X$. An alternative name for (4) is the “ridgeless” least squares estimator, motivated by the fact that $\hat{\beta}=\lim_{\lambda\to 0^{+}}\hat{\beta}_{\lambda}$, where $\hat{\beta}_{\lambda}$ denotes the ridge regression estimator:

or, equivalently, $\hat{\beta}_{\lambda}=(X^{T}X+n\lambda I)^{-1}X^{T}y$.

When $X$ has full column rank the min-norm least squares estimator reduces to $\hat{\beta}=(X^{T}X)^{-1}X^{T}y$, the usual least squares estimator. When $X$ has rank $n$, importantly, this estimator interpolates the training data: $y_{i}=x_{i}^{T}\hat{\beta}$, for $i=1,\ldots,n$.

Lastly, the following is a well-known fact that connects the min-norm least squares solution to gradient descent (as referenced in the introduction).

Initialize $\beta^{(0)}=0$, and consider running gradient descent on the least squares loss, yielding iterates

where we take $0<t\leq 1/\lambda_{\max}(X^{T}X)$ (and $\lambda_{\max}(X^{T}X)$ is the largest eigenvalue of $X^{T}X$). Then $\lim_{k\to\infty}\beta^{(k)}=\hat{\beta}$, the min-norm least squares solution in (4).

3 Bias and variance

We recall expressions for the bias and variance of the min-norm least squares estimator, which are standard.

Under the model (1), (2), the min-norm least squares estimator (4) has bias and variance

where $\hat{\Sigma}=X^{T}X/n$ is the (uncentered) sample covariance of $X$, and $\Pi=I-\hat{\Sigma}^{+}\hat{\Sigma}$ is the projection onto the null space of $X$.

4 Underparametrized asymptotics

We consider an asymptotic setup where $n,p\to\infty$, in such a way that $p/n\to\gamma\in(0,\infty)$. Recall that when $\gamma<1$, we call the problem underparametrized; when $\gamma>1$, we call it overparametrized. Here, we recall the risk of the min-norm least squares estimator in the underparametrized case. The rest of this paper focuses on the overparametrized case.

The following is a known result in random matrix theory, and can be found in Chapter 6 of Serdobolskii (2007). It can also be found in the wireless communications literature, see Chapter 4 of Tulino and Verdu (2004).

Assume the model (1), (2), and assume $x\sim P_{x}$ is of the form $x=\Sigma^{1/2}z$, where $z$ is a random vector with i.i.d. entries that have zero mean, unit variance, and a finite 4th moment, and $\Sigma$ is a (sequence of) deterministic positive definite matrix, such that $\lambda_{\min}(\Sigma)\geq c>0$, for all $n,p$ and a constant $c$ (here $\lambda_{\min}(\Sigma)$ is the smallest eigenvalue of $\Sigma$). Then as $n,p\to\infty$, such that $p/n\to\gamma<1$, the risk of the least squares estimator (4) satisfies, almost surely,

As it can be seen from the last proposition, in the underparametrized case the risk is just variance. In contrast, in the overparametrized case, the bias $B_{X}(\hat{\beta};\beta)=\beta^{T}\Pi\Sigma\Pi\beta$ is nonzero, because $\Pi$ is. This will be the focus of the next sections.

Isotropic features

We begin by considering the simpler case in which $\Sigma=I$. In this case the limiting bias is relatively straightforward to compute and depends on $\beta$ only through $\|\beta\|_{2}^{2}$. In Section 4, we generalize our analysis and study the dependence of the bias on the geometry of $\Sigma$ and $\beta$.

Choosing $U$ so that $U\beta=re_{i}$, the $i$th standard basis vector, then averaging over $i=1,\ldots,p$, yields

It is possible to show that, $B_{X}(\hat{\beta};\beta)$ concentrates around its expectation ant therefore $B_{X}(\hat{\beta};\beta)\to r^{2}(1-1/\gamma)$, almost surely. This is stated formally in the next section.

2 Limiting risk

As the next result shows, the independence of the risk on $\beta$ is still true outside of the Gaussian case, provided the features are isotropic. The next result can be proved as a corollary of the more general Theorem 3 below. We give a simpler self-contained proof using a theorem of Rubio and Mestre (2011) in Appendix C.2.

Assume the model (1), (2), where $x\sim P_{x}$ has i.i.d. entries with zero mean, unit variance, and a finite moment of order $4+\eta$, for some $\eta>0$. Also assume that $\|\beta\|_{2}^{2}=r^{2}$ for all $n,p$. Then for the min-norm least squares estimator $\hat{\beta}$ in (4), as $n,p\to\infty$, such that $p/n\to\gamma\in(0,\infty)$, it holds almost surely that

Hence, summarizing with Proposition 2, we have

We can see that the limiting norm, as a function of $\gamma$, has a somewhat similar profile to the limiting risk in (8): it is monotonically increasing on $(0,1)$, diverges at the interpolation boundary, and is monotonically decreasing on $(1,\infty)$.

Correlated features

We broaden the scope of our analysis from the last section, where we examined isotropic features. In this section, we take $x\sim P_{x}$ to be of the form $x=\Sigma^{1/2}z$, where $z$ is a random vector with independent entries that have zero mean and unit variance, and $\Sigma$ is arbitrary (but still deterministic and positive definite).

We next state our assumptions about the data distribution: our results will be uniform with respect to the (large) constants $M$, $\{C_{k}\}_{k\geq 2}$ appearing in this assumption.

The covariates vector $x\sim P_{x}$ is of the form $x=\Sigma^{1/2}z$, where, defining $\widehat{H}_{n}$ as per Eq. (9), we have

$s_{1}=\|\Sigma\|_{op}\leq M$, $\int s^{-1}d\widehat{H}_{n}(s)<M$.

$|1-(p/n)|\geq 1/M$, $1/M\leq p/n\leq M$.

Condition $(a)$ bounds the tail probabilities on the covariates. Requiring finite moment of all order is useful to get strong bounds on the deviations of the risk from its predicted value. As discussed below, bounds on the first few moments are sufficient if we are satisfied in weaker probability bounds.

Conditions $(b)$ requires the eigenvalues of $\Sigma$ to be bounded, and not to accumulateThe latter assumption could have been further relaxed, by requiring only $p_{+}$ of the $p$ eigenvalues to be non-vanishing and to satisfy the other conditions. This would require to redefine $\gamma$ as $p_{+}/n$. near . For the analysis of min-norm interpolation, we will add the additional assumption that the minimum eigenvalue of $\Sigma$ is bounded away from zero. However condition $(b)$ is sufficient for the analysis of ridge regression in Section 6.

Finally, as our statements are non-asymptotic, we do not assume $p/n$ to converge to to a value. However condition $(c)$ requires $p/n$ to be bounded and bounded away from the interpolation threshold $p/n=1$.

We then define the predicted bias and variance by

Note that evaluating $\mathscrsfs{B}(H,G,\gamma)$, $\mathscrsfs{V}(H,\gamma)$ numerically is relatively straightforward, with the most complex part being the solution of Eq. (10). The next theorem establishes that —under suitable technical conditions— the functions $\mathscrsfs{B}$, $\mathscrsfs{V}$ characterize the test error. Similar theorems were proved in Wu and Xu (2020); Richards et al. (2020), which generalized an earlier version of this manuscript to account for the geometry of $(\Sigma,\beta)$.

Assume the data model (1), (2) and that the covariates distribution satisfies Assumption 1. Further assume $s_{p}=\lambda_{\min}(\Sigma)>1/M$. Define $\gamma=p/n$ and let $\hat{\beta}$ be the min-norm least squares estimator in Eq. (4).

Then for any constants $D>0$ (arbitrarily large) there exist $C=C(M,D)$ such that, with probability at least $1-Cn^{-D}$ the following hold

where $\mathscrsfs{B}$ and $\mathscrsfs{V}$ are given in Definition 2.

The proof of this theorem is deferred to Section B.

The order of the error bound in Eqs. (14), (15) is not optimal: a central limit theorem heuristics suggests the deterministic approximation to be accurate up to an error of order $n^{-1/2}$. Indeed, we are able to establish the correct order in the case of ridge regression, see Theorem 5.

Note that Theorem 2 establishes deterministic approximations for the bias and variance, that are valid at finite $n,p$. The overparametrization ratio $\gamma=p/n$ is a non-asymptotic quantity, and the error bounds are uniform, i.e. depend uniquely on the constant $M$. This is to be contrasted with the asymptotic setting of Richards et al. (2020); Wu and Xu (2020). Both of these papers assume a sequence of regression problems with $n,p\to\infty$, $p/n\to\gamma$, and obtain an asymptotically exact expression for the risk.

Technically, Richards et al. (2020); Wu and Xu (2020) apply asymptotic random matrix theory results, such as Ledoit and Peche (2011), while we have to take a longer detour to exploit non-asymptotic results established inKnowles and Yin (2017). We believe that the non-asymptotic approach provides more concrete and accurate statements.

Theorem 2 also implies asymptotic predictions under minimal assumptions. In particular, if the two probability measures $\widehat{H}_{n}$, $\widehat{G}_{n}$ converge weakly to probability measures $H$, $G$ on $[0,\infty)$, then we obtain $B_{X}(\hat{\beta};\beta)/\|\beta\|^{2}\to\mathscrsfs{B}(H,G,\gamma)$, $V_{X}(\hat{\beta};\beta)\to\mathscrsfs{V}(H,\gamma)$ almost surely.

As pointed out above the condition $\widehat{H}_{n}\Rightarrow H$, $\widehat{G}_{n}\Rightarrow G$ (here $\Rightarrow$ denotes weak convergence) is strictly weaker than the condition assumed in Wu and Xu (2020) to establish asymptotic results. Further, we require weaker moment conditions.

In the next sections we illustrate the role of the geometry of $\beta,\Sigma$ by considering two models for which $\widehat{H}_{n}\Rightarrow H$, $\widehat{G}_{n}\Rightarrow G$ as $n,p\to\infty$. First we consider the case $G=H$, which we refer to as ‘equidistributed:’ the components of $\beta$ are roughly equally distributed along the eigenvectors of $\Sigma$. In this case there is no special relation between $\beta$ and $\Sigma$.

As a further application, we consider a latent space model in which $\beta$ is aligned with the top eigenvectors of $\Sigma$. This can be regarded as a misspecified model, and is therefore presented in Section 5.4 below.

2 Equidistributed coefficients

In this section we assume $G=H$, $\|\beta\|_{2}\to r$, and $p/n\to\gamma$. One way to generate $\beta$ satisfying this condition is to draw it uniformly at random on the $p$-dimensional sphere of radius $\|\beta\|_{2}=r$. In this case the conditions of Theorem 2 (or Theorem 3), hold with $\widehat{G}_{n}\to G=H$.

Under the assumptions of Theorem 3, further assume $G=H$, $\|\beta\|^{2}\to r^{2}$. Then for $n,p\to\infty$, with $p/n\to\gamma>1$, almost surely

As a special case, we can revisit the isotropic case $\Sigma=I$, which results in $dH=\delta_{1}$. In this case $c_{0}(H,\gamma)=\gamma(\gamma-1)$ yielding immediately $\mathscrsfs{B}_{\mbox{\tiny\rm equi}}(H,\gamma)=1-\gamma^{-1}$ and $\mathscrsfs{V}_{\mbox{\tiny\rm equi}}(H,\gamma)=1/(\gamma-1)$.

Under the assumptions of Theorem 3, further assume $\|\beta\|^{2}\to r^{2}$, and let $c_{0}=c_{0}(H,\gamma)$ be defined as there. Then as $n,p\to\infty$, such that $p/n\to\gamma$, the min-norm least squares estimator (4) satisfies, almost surely,

Misspecified model

In this section, we consider a misspecified model, in which the regression function is still linear, but we observe only a subset of the features. Such a setting provides another potential motivation for interpolation: in many problems, we do not know the form of the regression function, and we generate features in order to improve our approximation capacity. Increasing the number of features past the point of interpolation (increasing $\gamma$ past 1) can now decrease both bias and variance.

As such, the misspecified model setting also yields further interesting asymptotic comparisons between the $\gamma<1$ and $\gamma>1$ regimes. Recall the isotropic features model of Section 3.2: the risk function in (8) is globally minimized at $\gamma=0$. This is a consequence of the fact that, in a well-specified linear model at $\gamma=0$ there is no bias and no variance, and hence no risk. In a misspecified model, we will see that the story can be quite different, and the asymptotic risk can actually attain its global minimum on $(1,\infty)$.

Consider, instead of (1), (2), a data model

Given a test point $(x_{0},w_{0})\sim P_{x,w}$, and an estimator $\hat{\beta}$ (fit using $X,y$ only, and not $W$), we define its out-of-sample prediction risk as

Note that this definition is conditional on $X$, and we are integrating over the randomness not only in $\epsilon$ (the training errors), but in the unobserved features $W$, as well. The next lemma decomposes this notion of risk in a useful way.

Under the misspecified model (19), (20), for any estimator $\hat{\beta}$, we have

The first term $R^{*}_{X}(\hat{\beta};\beta,\theta)$ in the decomposition in Lemma 2 is precisely the risk that we studied previously in the well-specified case, except that the response distribution has changed (due to the presence of the middle term in (20)). We call the second term $M(\beta,\theta)$ in Lemma 2 the misspecification bias.

In order to discuss some qualitative features, we focus on the simplest possible model by assuming independent covariates.

2 Isotropic features

Here, we make the additional simplifying assumption that $(x,w)\sim P_{x,w}$ has i.i.d. entries with unit variance, which implies that $\Sigma=I$. (The case of independent features but general covariances $\Sigma_{x},\Sigma_{w}$ is similar, and we omit the details.) Therefore, we may write the response distribution in (20) as

where $\delta_{i}$ is independent of $x_{i}$, having mean zero and variance $\sigma^{2}+\|\theta\|_{2}^{2}$, for $i=1,\ldots,n$. Denote the total signal by $r^{2}=\|\beta\|_{2}^{2}+\|\theta\|_{2}^{2}$, and the fraction of the signal captured by the observed features by $\kappa=\|\beta\|_{2}^{2}/r^{2}$. Then $R_{X}^{*}(\hat{\beta};\beta,\theta)$ behaves exactly as we computed previously, for isotropic features in the well-specified setting (Theorem 2 for $\gamma<1$, and Theorem 1 for $\gamma>1$), after we make the substitutions:

Furthermore, we can easily calculate the misspecification bias:

Putting these results together leads to the next conclusion.

Assume the misspecified model (19), (20), and assume $(x,w)\sim P_{x,w}$ has i.i.d. entries with zero mean, unit variance, and a finite moment of order $8+\eta$, for some $\eta>0$. Also assume that $\|\beta\|_{2}^{2}+\|\theta\|_{2}^{2}=r^{2}$ and $\|\beta\|_{2}^{2}/r^{2}=\kappa$ for all $n,p$. Then for the min-norm least squares estimator $\hat{\beta}$ in (4), as $n,p\to\infty$, with $p/n\to\gamma$, it holds almost surely that

We remark that, in the independence setting considered in Theorem 4, the dimension $q$ of the unobserved feature space does not play any role: we may equally well take $q=\infty$ for all $n,p$ (i.e., infinitely many unobserved features).

3 Polynomial approximation bias

Since adding features should generally improve our approximation capacity, it is reasonable to model $\kappa=\kappa(\gamma)$ as an increasing function of $\gamma$. To get an idea of the possible shapes taken by the asymptotic risk curve from Theorem 4, we consider the example of a polynomial decay for the approximation bias,

for some $a>0$. In this case, the limiting risk in the isotropic setting, from Theorem 4, becomes

4 A latent space model

Here $(\xi_{i})_{i\leq n}$, $(u_{ij})_{i\leq n,j\leq p}$ are noise variables that are mutually independent, and independent of $z_{i}$, with $\xi_{i}\sim{N}(0,\sigma_{\xi}^{2})$, $u_{ij}\sim{N}(0,1)$. The features matrix takes the form $X=ZW^{T}+U$ and therefore, for $p>n$, min-norm regression amounts to

Notice that the signal-to-noise ratio in the features $x_{ij}$ is equal to $\|w_{i}\|_{2}^{2}$. To be definite, we’ll fix $\|w_{i}\|_{2}=1$, which implies in particular $\|W\|_{F}^{2}=p$. In order to simplify our calculations, we assume all the non-zero singular values of $W$ to be equal, whence $W^{T}W=(p/d)I_{d}$. It is not hard to check that this model satisfies the assumptions of Theorem 2, with

Using Theorem 2, we get the following explicit expressions.

Consider the latent space model described above, and assume $d/p\to\psi\in(1,\infty)$, $p/n\to\gamma\in(0,\infty)$. Then, almost surely

where $\sigma^{2}=\sigma_{\xi}^{2}+\psi r_{\theta}^{2}/(1+\psi)$, and $c_{0}=c_{0}(\psi,\gamma)\geq 0$ is the unique non-negative solution of the following second order equation

Figures 5 and 6 illustrate this corollary by comparing analytical predictions to numerical simulations. We observe that the prediction risk is monotone decreasing in the overparametrization ratio for $\gamma>1$, and reaches its global minimum asymptotically as $\gamma\to\infty$ (after $p,n,d\to\infty$). To understand why this happens, notice that each feature vector $x_{i}$ can be viewed as a noisy measurement of the latent covariates $z_{i}$. If the noise $u_{ij}$ was absent, then performing min-norm regression with respect to $(x_{i})_{i\leq n}$ amounts to minimize would be equivalent to min-norm regression with respect to $(z_{i})_{i\leq n}$. To see this, consider again Eq. (25). If we drop the noise $U$, we are minimizing $\|b\|_{2}$ subject to $Z(W^{T}b)=0$, and the regression function is $\hat{f}(z)=x^{T}\hat{\beta}=z^{T}(W^{T}\hat{\beta})$. Since $W$ is orthogonal, this is equivalent to computing $\hat{\theta}=\arg\min\{\|t\|_{2}:$ subject to $Zt=0\}$.

In presence of noise $u_{ij}$, the latent features cannot be estimated exactly. However as $p$ gets larger, the noise is effectively ‘averaged out’ and we approach the idealized situation in which the $z_{i}$’s are observed.

Ridge regularization

We generalize the formulas of Section 4 to non-vanishing ridge regularization. We work under the same assumptions of that section. In particular, recall that $\widehat{H}_{n}(s)=p^{-1}\sum_{i=1}^{p}1_{\{s\geq s_{i}\}}$ is the empirical distribution of eigenvalue of $\Sigma$, and $\widehat{G}_{n}(s)=\sum_{i=1}^{p}\langle\beta,v_{i}\rangle^{2}1_{\{s\geq s_{i}\}}/\|\beta\|^{2}$ the same empirical distribution, reweighted by the projection of $\beta$ onto the eigenvector $v_{i}$ of the covariance $\Sigma$. (Recall the eigenvalue decomposition $\Sigma=\sum_{i=1}^{p}s_{i}v_{i}v_{i}^{T}$.)

Further define $m_{n,1}(z)=m_{n,1}(z;\widehat{H}_{n},\gamma)$ via

These definitions are extended analytically to ${\operatorname{Im}}(z)=0$ whenever possible. We then define the predicted bias and variance by

We next state our deterministic approximation of the risk.

Let $M^{-1}\leq p/n\leq M$, and Assumption 1 hold. Further assume $\lambda\vee\lambda_{\min}(\Sigma)>1/M$ and $n^{-2/3+1/M}<\lambda<M$. Let $\hat{\beta}_{\lambda}$ be the ridge estimator of Eq. (5).

Then for any constants $D>0$ (arbitrarily large) and ${\varepsilon}>0$ (arbitrarily small), there exist $C=C(M,D)$ such that, with probability at least $1-Cn^{-D}$ the following hold

where $\mathscrsfs{B}$ and $\mathscrsfs{V}$ are given in Definition 2.

The proof of this theorem is deferred to Appendix A. As for theorem 2, similar results were proved in Wu and Xu (2020); Richards et al. (2020), subsequently to a first version of this manuscript that only focused on random $\beta$. The same comparison of Remark 1 applies here.

In particular, Theorem 5 establishes non-asymptotic deterministic approximations for the bias $B_{X}(\hat{\beta}_{\lambda};\beta)$ and variance $V_{X}(\hat{\beta};\beta)$. The error terms are uniform over the covariance matrix, and have nearly optimal dependence upon the sample size $n$. Indeed, a central-limit theorem heuristics suggests indeed fluctuations of order $n^{-1/2}$.

As for the case of min-norm regression, it is easy to extract asymptotic prediction under minimal assumptions.

As a special case, we can consider the simple isotropic model that was already studied in Section 3. Very similar (though not identical) results can be found in Dicker (2016); Dobriban and Wager (2018).

Assume the conditions of Theorem 1 (well-specified model, isotropic features). Then for ridge regression estimator in (5) as $n,p\to\infty$, such that $p/n\to\gamma\in(0,\infty)$, it holds almost surely that

Here $m(z)$ is given by Eq. (36) which in this case has the explicit solution $m(z)=[1-\gamma-z-\sqrt{(1-\gamma-z)^{2}-4\gamma z}]/(2\gamma z)$.

Furthermore, the limiting ridge risk is minimized at $\lambda^{*}=\sigma^{2}\gamma/r^{2}$, in which case we have the simpler expression $R_{X}(\hat{\beta}_{\lambda};\beta)\to\sigma^{2}\gamma m(-\lambda^{*})$.

Under the conditions of Theorem 4 (misspecified model, isotropic features), the limiting risk of ridge regression is as in the first two displays, and the optimal limiting risk is as in the third, after we make the substitutions in (21) and add $r^{2}(1-\kappa)$, to each expression.

It is easy to recover the formulas in Theorem 1 as a limiting case of Eq. (43), by using the $z\to 0$ asymptotics $m(z)=(1-\gamma)^{-1}+O(z)$ for $\gamma<1$ and $m(z)=(1-\gamma)^{-1}+O(z)$ for $\gamma<1$ and $m(z)=-(\gamma-1)/(\gamma z)+[(\gamma-1)\gamma]^{-1}+O(z)$ for $\gamma>1$.

Figures 8 and 8 compare the risk curves of min-norm least squares to those from optimally-tuned ridge regression, in the well-specified and misspecified settings, respectively. There are two important points to make. The first is that optimally-tuned ridge regression is seen to have strictly better asymptotic risk throughout, regardless of $r^{2}$, $\gamma$, $\kappa$. This should not be a surprise, as by definition optimal tuning should yield better risk than min-norm least squares, which is the special case given by $\lambda\to 0^{+}$.

As we will see, the behavior is rather different in the latent space model.

2 Latent space model

As a special application, we consider the latent space model of Section 5.4. It is immediate to specialize Eqs. (38) and (39) to this case. We omit giving giving explicit formulas for brevity, and instead plot the resulting curves for the prediction risk (test error).

In Figure 10 we plot the risk as a function of the overparametrization ration $\gamma=p/n$ for several values of the regularization parameter $\lambda$ (included the ridgeless limit $\lambda\to 0$). The setting here is analogous to the one of Fig. 5. We observe several interesting phenomena:

Independently of $\lambda$ in the probed range, the risk is minimized at large overparametrization $\gamma\gg 1$.

As expected, the divergence of the risk at the interpolation threshold $\gamma=1$ is smoothed out by regularization, and the risk becomes a monotone decreasing function of $\gamma$ when $\lambda$ is large enough. Crucially, the optimal amount of regularization (corresponding to the lower envelope of these curves) results in a monotonically decreasing risk.

At large overparametrization, the optimal value of the regularization parameter is $\lambda\to 0$.

We confirm the last finding in Fig. 10 which plots the risk as a function of $\lambda$: the optimal regularization is $\lambda\to 0$. Notice that this is the case despite the fact that the observations are noisy, namely $\sigma_{\xi}>0$ strictly.

The optimality of $\lambda\to 0$ has a known intuitive explanation that is worth recalling here. Recall that ridge predictor at point $x_{0}$ takes the form

where we introduced the kernel matrix $K(X,X)=XX^{T}/n$, and the vector $K(x_{0},X)=x_{0}^{T}X/n$. In the present case $X$ can ve viewed as a noisy version of $\overline{X}=ZW^{T}$, namely Consider the case in which the covariates $X$ contain noise, as is our case, $X=\overline{X}+\eta U$ where $u_{ij}\sim N(0,1)$. (While we are considering $\eta=1$, it is instructive to regard the noise standard deviation as a parameter.) Then, we might expect $K(X,X)\approx K(\overline{X},\overline{X})+\eta^{2}UU^{T}\approx K(\overline{X},\overline{X})+\eta^{2}I_{n}$. If this approximations hold, the noise acts as an extra ridge term, which can be sufficient to regularize the problem.

To the best of our knowledge, this argument was first presented by Webb (1994) and, more explicitly, by Bishop (1995). Recently, Kobak et al. (2020) elucidated its role in linear regression, establishing several of its consequences. In a parallel line of work, a closely related idea has recently emerged in the analysis of kernel methods in high dimension El Karoui (2010); Liang et al. (2020).

Cross-validation

We analyze the effect of using cross-validation to choose the tuning parameter in ridge regression. In short, we find that choosing the ridge tuning parameter to minimize the leave-one-out cross-validation error leads to the same asymptotic risk as the optimally-tuned ridge estimator. The next subsection gives the details; the following subsection presents a new “shortcut formula” for leave-one-out cross-validation in the overparametrized regime, for min-norm least squares, akin to the well-known formula for underparametrized least squares and ridge regression. We refer to Craven and Wahba (1979); Golub et al. (1979a) for background on CV and GCV, and to Arlot and Celisse (2010) for a more recent review.

Recomputing the leave-one-out predictors $\hat{f}^{-i}_{\lambda}$, $i=1,\ldots,n$ can be burdensome, especially for large $n$. Importantly, there is a well-known “shortcut formula” that allows us to express the leave-one-out CV error (45) as a weighted average of the training errors,

where $S_{\lambda}=X(X^{T}X+n\lambda)^{-1}X^{T}$ is the ridge smoother matrix. This identity is an immediate consequence of the Sherman-Morrison-Woodbury formula. In the next subsection, we give an extension to the case $\lambda=0$ and ${\rm rank}(X)=n$, i.e., to min-norm least squares.

The next result shows that, for isotropic features, the CV error of a ridge estimator converges almost surely to its prediction error. The focus on isotropic features is only for simplicity: a more general analysis is possible but is not pursued here. The proof, given in Appendix C.6, relies on the shortcut formula (46). In the proof, we actually first analyze generalized cross-validation (GCV), which turns out to be somewhat of an easier calculation (see the proof for details on the precise form of GCV), and then relate leave-one-out CV to GCV.

with the right-hand side above being the asymptotic risk of optimally-tuned ridge regression. Further, the exact same set of results holds for GCV.

Similar results were obtained for various linear smoothers in Li (1986, 1987), for the lasso in the high-dimensional (proportional) asymptotics in Miolane and Montanari (2018), and for general smooth penalized estimators in Xu et al. (2019). The latter paper covers ridge regression as a special case, and gives more precise results (convergence rates), but assumes more restrictive conditions.

The key implication of Theorem 7, in the context of the current paper and its central focus, is that the CV-tuned or GCV-tuned ridge estimator has the same asymptotic performance as the optimally-tuned ridge estimator. In other words, the ridge curves in Figures 1, 8, and 8 can be alternatively viewed as the asymptotic risk of ridge under CV tuning.

2 Shortcut formula for ridgeless CV

We extend the leave-one-out CV shortcut formula (46) to work when $p>n$ and $\lambda=0+$, i.e., for min-norm least squares. In this case, both the numerator and denominator are zero in each summand of (46). To circumvent this, we can use the so-called “kernel trick” to rewrite the ridge regression solution (5) with $\lambda>0$ as

Using this expression, the shortcut formula for leave-one-out CV in (46) can be rewritten as

Taking $\lambda\to 0^{+}$ yields the shortcut formula for leave-one-out CV in min-norm least squares (assuming without a loss of generality that ${\rm rank}(X)=n$),

In fact, the exact same arguments given here still apply when we replace $XX^{T}$ by a positive definite kernel matrix $K$ (i.e., $K_{ij}=k(x_{i},x_{j})$ for each $i,j=1,\ldots,n$, where $k$ is a positive definite kernel function), in which case (48) gives a shortcut formula for leave-one-out CV in kernel ridgeless regression (the limit in kernel ridge regression as $\lambda\to 0^{+}$). We also remark that, when we include an unpenalized intercept in the model, in either the linear or kernelized setting, the shortcut formula (48) still applies with $XX^{T}$ or $K$ replaced by their doubly-centered (row- and column-centered) versions, and the matrix inverses replaced by pseudoinverses.

Nonlinear model

Our analysis in the previous sections assumed $x_{i}=\Sigma^{1/2}z_{i}$, with $z_{i}$ a vector with i.i.d. entries. As discussed in the introduction, we expect our results to have implications on certain neural networks in the lazy training regime, via the universality hypothesis. Roughly speaking, if covariates are obtained by applying the featurization map $z_{i}\mapsto x_{i}=\nabla f(z_{i};\theta_{0})$ to data $z_{i}$, the asymptotic behavior of ridge regression is expected to be close to the one of an equivalent Gaussian model, whose second-order statistics match the ones of $\nabla f(z_{i};\theta_{0})$.

We focus on a simple case in which the second order statistics match the ones of the isotropic model.

Notice that, conditionally on $W$, the vectors $x_{i}=\varphi(Wz_{i})$, $i\leq n$ are independent. However, they do not have independent coordinates: intuitively if $d$ is significantly smaller than $p$, $x_{i}$ contains much less randomness than a vector of $p$ independent entriesFor instance, if $\varphi(t)=at^{2}+b$, we can reconstruct $z_{i}$ from the first $2d$ coordinates of $x_{i}$, and therefore the remaining $p-2d$ coordinates of $x_{i}$ are a function of the first $2d$..

Then for $\gamma>1$, the variance satisfies, almost surely,

which is again as in the case of linear isotropic features, recall Theorem 1.

The proof of Theorem 8 is lengthy and will be sketched shortly.

where $G_{1},G_{2}$ are jointly Gaussian with unit variance and covariance $w_{i}^{T}w_{j}$. Denoting by $\varphi(x)=\sum_{k\geq 0}\lambda_{k}(\varphi)\,h_{k}(x)$ the decomposition of $\varphi$ into orthonormal Hermite polynomials, we thus obtain

Despite the nonlinear model $x_{i}=\varphi(Wz_{i})$ matches the second-order population statistics of the isotropic model, it is far from obvious that the resulting least norm regression risk is the same. Indeed the coordinates of vector $x_{i}$ are highly dependent, as it can be seen by noting that they are a function of only $d\ll p$ independent random variables. Theorem 8 confirms that —despite dependence— the risk is asymptotically the same, thus providing a concrete example of the general universality phenomenon.

2 Proof outline for Theorem 8

We introduce the following resolvents (as usual, these are defined for ${\operatorname{Im}}(\xi)>0$ and by analytic continuation, whenever possible, for ${\operatorname{Im}}(\xi)=0$):

Here and henceforth, we write $[i,j]=\{i+1,\ldots,i+j\}$ for integers $i,j$. We also write $M^{-1}_{ij}=(M^{-1})_{ij}$ for a matrix $M$, and ${\sf Tr}_{S}(M)=\sum_{i\in S}M_{ii}$ for a subset $S$. The equalities in the first and third lines above follow by invariance of the distribution of ${A}(s)$ under permutations of $[1,p]$ and $[p+1,p+n]$. Whenever clear from the context, we will omit the arguments from block matrix and resolvents, and write ${A}={A}(s)$, $m_{1,n}=m_{1,n}(\xi,s)$, and $m_{2,n}=m_{2,n}(\xi,s)$.

The next lemma characterizes the asymptotics of $m_{1,n}$, $m_{2,n}$.

Assume the conditions of Theorem 8. Consider ${\operatorname{Im}}(\xi)>0$ or ${\operatorname{Im}}(\xi)=0$, ${\operatorname{Re}}(\xi)<0$, with $s\geq t\geq 0$. Let $m_{1}$ and $m_{2}$ be the unique solutions of the following quadratic equations:

subject to the condition of being analytic functions for ${\operatorname{Im}}(z)>0$, and satisfying $|m_{1}(z,s)|,|m_{2}(z,s)|\leq 1/{\operatorname{Im}}(z)$ for ${\operatorname{Im}}(z)>C$ (with $C$ a sufficiently large constant). Then, as $n,p,d\to\infty$, such that $p/n\to\gamma$ and $d/p\to\psi$, we have almost surely (and in $L^{1}$),

The proof of this lemma is given in Appendix D.2. As a corollary of the above, we obtain that the asymptotical empirical spectral distribution of the empirical covariance $\hat{\Sigma}=X^{T}X/n$ matches the one for the independent entries model, and is hence given by the Marchenko-Pastur law (a result already obtained in Pennington and Worah (2017)). We state this formally using the Stieltjes transform

Assume the conditions of Theorem 8. Consider ${\operatorname{Im}}(z)>0$. As $n,p,d\to\infty$, with $p/n\to\gamma$ and $d/p\to\psi$, we have (almost surely and in $L^{1}$) $R_{n}(\xi)\to r(\xi)$ where $r$ is nonrandom and coincides with the Stieltjes transform of the Marchenko-Pastur law, namely

We refer to Appendix D.4 for a proof of this corollary. The next lemma connects the above resolvents computed in Lemma 3 to the variance of min-norm least squares, hence finishing our proof outline.

Assume the conditions of Theorem 8. Let $m_{1}$, $m_{2}$ be the asymptotic resolvents given in Lemma 3. Define

Then for $\gamma\neq 1$, \partial_{x}m(\xi,x)\big{|}_{x=0} as a simple pole at $\xi=0$, and hence admits a Taylor-Laurent expansion around $\xi=0$, whose coefficients will be denoted by $D_{-1}$, $D_{0}$

Here each coefficient is a function of $\gamma,\psi$: $D_{-1}=D_{-1}(\gamma,\psi)$, $D_{0}=D_{0}(\gamma,\psi)$. Furthermore, for the ridge regression estimator $\hat{\beta}_{\lambda}$ in (5), as $n,p,d\to\infty$, such that $p/n\to\gamma\in(0,\infty)$, $d/p\to\psi\in(0,1)$, the following ridgeless limit holds almost surely:

The proof of this lemma can be found in Appendix D.3. Theorem 8 follows by evaluating the formula in Lemma 4, by using the result of Lemma 3. We refer to the appendices material for details.

Acknowledgements

The authors are grateful to Brad Efron, Rob Tibshirani, and Larry Wasserman for inspiring us to work on this in the first place. RT sincerely thanks Edgar Dobriban for many helpful conversations about the random matrix theory literature, in particular, the literature trail leading up to Proposition 2, and the reference to Rubio and Mestre (2011). We are also grateful to Dmitry Kobak for bringing to our attention the Kobak et al. (2020) and clarifying the implications of this work on our setting.

TH was partially supported by grants NSF DMS-1407548, NSF IIS-1837931, and NIH 5R01-EB001988-21. AM was partially supported by NSF DMS-1613091, NSF CCF-1714305, NSF IIS-1741162, and ONR N00014-18-1-2729. RT was partially supported by NSF DMS-1554123.

Appendix A Proofs for the linear model: Risk of ridge regression

Before passing to the proof of Theorem 5, we state and prove a useful lemma in analysis.

Without loss of generality we can assume $x=0$. Further we set $h=f-g$. We consider a finite difference approximation of order $2k$ of the first order derivative of $h$ at . One such approximation is given by (Trefethen, 1996, Chapter 3)

where $u_{x}$ is a point in the interval $[0,x]$. Using the properties of $D_{2k}$ listed above, we obtain

We are left wih the task of proving claims $(i)$ and $(ii)$ about the operator $D_{2k}$. Substituting, these are equivalent to the claims

We further state and proof a lemma about a modified Stieltjes transform.

Let $s_{1}\geq\dots\geq s_{p}>0$ be such that $s_{1}\leq M$ and $\#\{i:\,s_{i}\geq 1/M\}\geq p/M$ for some large constant $M>1$, and $M^{-1}\leq p/n\leq M$. Further assume $|1-(p/n)|\geq 1/M$ and $n^{-1}\sum_{i=1}^{p}s_{i}^{-1}1_{s_{i}>0}<M$. For $(\lambda,\eta)\in(0,M]\times(-1/(10M^{2}),1/(10M^{2}))$ let $u=u(\lambda,\eta)$ be the only solution in $(0,\infty)$ of the equation

$\lambda\mapsto u(\lambda,\eta)$ is analytic and non-decreasing with $0<\partial_{\lambda}u(\lambda,\eta)\leq C(M)$ for a constant $C(M)$ uniquely dependent on $M$.

If further $\lambda\vee s_{p}\geq 1/M$, then $|\partial^{k}_{\eta}u(\lambda,\eta)|\leq C_{k}(M)$ for all $k$.

We begin from poinr $(a)$. Since $\eta$ is fixed here, we write $u(\lambda)=u(\lambda,\eta)$, dropping the dependence on $\eta$. Note that, for $\eta>0$, any solution $u$ must be such that $u\in(0,1/\eta)$ because otherwise the left hand side is negative and the right-hand side is positive. Call $x=u/(1-\eta u)$ whence the above equation is equivalent to the following equation on $x\in(0,\infty)$:

where $s_{i}(\eta)=s_{i}/(1+\eta s_{i})$. For $\eta\in(-1/(10M^{2}),1/(10M^{2}))$, we have $(1-(10M)^{-1})s_{i}<s_{i}(\eta)<(1+(10M)^{-1})s_{i}$, and therefore the sequence $(s_{i}(\eta))_{i\leq p}$ satisfies the same assumptions as $(s_{i})_{i\leq p}$, with $M$ replaced by $2M$. We will therefore drop the argument $\eta$ and study the fixed point equation (61) (which is equivalent to the case $\eta=0$ of Eq. (60)).

Let $g(x):=n^{-1}\sum_{i=1}^{p}s_{i}/(s_{i}+x)$, and note that this function is mononone decreasing on $[0,\infty)$ with $g(0)\leq C_{0}(M)$, and $g(x)\leq M^{2}/x$. Equation (61) can be rewritten as

Since $x\mapsto\hat{g}(x):=(\lambda/x)+g(x)$ is monotone decreasing with $\hat{g}(0+)=+\infty$ and $\hat{g}(x)\to 0$ as $x\to\infty$, this equation has a unique solution, denoted by $x(\lambda)$. By the implicit function theorem, we have

The second inequality follows for instance by noting that $n^{-1}\sum_{i=1}^{p}s_{i}^{2}/(s_{i}+x(\lambda))^{2}<g(x(\lambda))=1-\lambda/x(\lambda)<1$. Notice that $g(0)=p/n$. By assumption $|g(0)-1|\geq 1/M$. We will distinguish the two cases $g(0)<1$ and $g(0)>1$.

Case $g(0)<1$. Since $g(x)$ is decreasing, (62) implies $1\leq\lambda x^{-1}+g(0)$ whence

and therefore, by Eq. (63), we have $\partial_{\lambda}x(\lambda)\leq M$ for all $\lambda\in(0,\infty)$.

Case $g(0)>1$. Let $x_{0}=g^{-1}(1)>0$. Since $x\mapsto g(x)$ is convex, $x_{0}\geq(1-g(0))/|g^{\prime}(0)|$. Note that $|g^{\prime}(0)|=n^{-1}\sum_{i=1}^{p}s_{i}^{-1}1_{s_{i}>0}\leq M$. Therefore $x_{0}\geq 1/C(M)$ for some finite constant $C(M)$.

Since $\hat{g}(x)=\lambda/x+g(x)>g(x)$ on $(0,\infty)$, we have $x(\lambda)>x(0+)=x_{0}\geq 1/C(M)$ for all $\lambda\in(0,\infty)$. Further,

Therefore by Eq. (63), $0<\partial_{\lambda}x(\lambda)\leq C(M)$ also in this case.

and therefore $7/10\leq(1+\eta x(\lambda))\leq 13/10$, whence

Claim $(b)$ simply follows since $u/(1-\eta u)\geq\lambda$ by Eq. (60). Since $(1-\eta u)=1/(1+\eta x)\geq 1/2$ as just shown, this implies $u\geq\lambda/2$.

Finally, for claim $(c)$, notice that $u(\lambda,\eta)$ is a solution of the equation

By the implicit function theorem, it is sufficient to prove that \big{|}\frac{\partial^{k+l}f}{\partial u^{k}\partial\eta^{l}}\big{|}\leq C_{k,l}(M) for all $l\geq 1,k\geq 0$ and $l=0$, $k\geq 2$, and that \big{|}\frac{\partial f}{\partial u}\big{|}\geq 1/C(M).

In order to bound the first derivative away from , note that, again by the implicit function theorem,

where the last inequality follows from point $(a)$ above. In order to bound higher order derivatives, write $f(u,\eta)=f_{1}(\eta u)/\eta-\lambda-f_{2}(u)$ where $f_{1}(v):=v/(1-v)$. Recall that $x=u/(1-\eta u)$, and $x\leq 3M^{2}$ by Eq. (66). Therefore

where $C_{0}$, $c_{1}(m_{1},m_{2},m_{3})$ are coefficients depending uniquly on $k,l$. This sum is bounded since each term is bounded, by the above estimate on $\partial^{k}f_{1}$, and further using the fact that $u=x/(1+\eta x)\leq 10M^{2}$ since $x\leq 3M^{2}$ by Eq. (66). ∎

For the proof of Theorem 5, it is convenient to introduce the notations $S_{Z}:=Z^{T}Z/n$ and $S_{X}:=X^{T}X/n=\Sigma^{1/2}S_{Z}\Sigma^{1/2}$. We begin by recalling the expressions for bias and variance of ridge regression at regularization parameter $\lambda$: