The distribution of Ridgeless least squares interpolators

Consider the standard linear regression model

where we observe i.i.d. feature vectors $X_{i}\in\mathbb{R}^{n}$ and responses $Y_{i}\in\mathbb{R}$, and $\xi_{i}$’s are unobservable errors. For notational simplicity, we write $X=[X_{1}\cdots X_{m}]^{\top}\in\mathbb{R}^{m\times n}$ as the design matrix that collects all the feature vectors, and $Y=(Y_{1},\ldots,Y_{m})^{\top}\in\mathbb{R}^{m}$ as the response vector. The feature vectors $X_{i}$’s are assumed to satisfy $\operatorname{\mathbb{E}}X_{1}=0$ and $\operatorname{Cov}(X_{1})=\Sigma$, and the errors satisfy $\operatorname{\mathbb{E}}\xi_{1}=0$ and $\operatorname{Var}(\xi_{1})=\sigma_{\xi}^{2}$.

Throughout this paper, we reserve $m$ for the sample size, and $n$ for the signal dimension. The aspect ratio $\phi$, i.e., the number of samples per dimension, is then defined as $\phi\equiv m/n$. Accordingly, we refer to $\phi^{-1}>1$ as the overparametrized regime, and $\phi^{-1}<1$ as the underparametrized regime.

Within the linear model (1.1), the main object of interest is to recover/estimate the unknown signal $\mu_{0}\in\mathbb{R}^{n}$. While a large class of regression techniques can be used for the purpose of signal recovery under various structural assumptions on $\mu_{0}$, here we will focus our attention on one widely used class of regression estimators, namely, the Ridge estimator (cf. [HK70]) with regularization $\eta>0$,

and the Ridgeless estimator (also known as the minimum-norm interpolator),

which is almost surely (a.s.) well-defined in the overparametrized regime $\phi^{-1}>1$. Here $A^{-}$ is the Moore-Penrose pseudo-inverse of $A$. The notation $\widehat{\mu}_{0}$ is justified since for $\phi^{-1}>1$, $\widehat{\mu}_{\eta}\to\widehat{\mu}_{0}$ a.s. as $\eta\downarrow 0$.

From a conventional statistical point of view, the Ridgeless estimator seems far from an obviously good choice: As $\widehat{\mu}_{0}$ perfectly interpolates the data, it is susceptible to high variability due to the widely recognized bias-variance tradeoff inherent in ‘optimal’ statistical estimators [JWHT21, DSH23]. On the other hand, as the Ridgeless estimator $\widehat{\mu}_{0}$ is the limit point of the gradient descent algorithm run on the squared loss in the overparametrized regime $\phi^{-1}>1$, it provides a simple yet informative test case for understanding one major enigma of modern machine learning methods: these methods typically interpolate training data perfectly; still, they enjoy good generalization properties [JGH18, DZPS18, AZLS19, BHMM19, COB19, ZBH⁺21].

Inspired by this connection, recent years have witnessed a surge of interest in understanding the behavior of the Ridgeless estimator $\widehat{\mu}_{0}$ and its closely related Ridge estimator $\widehat{\mu}_{\eta}$, with a particular focus on their prediction risks, cf. [TV⁺04, EK13, HKZ14, Dic16, DW18, EK18, ASS20, BHX20, WX20, BLLT20, BMR21, RMR21, HMRT22, TB22, CM22]. An emerging picture from these works is that, the norm minimizing property of the Ridgeless estimator $\widehat{\mu}_{0}$ inherited from the gradient descent algorithm, induces certain ‘implicit regularization’ that helps prediction accuracy despite interpolation. More interestingly, in certain scenarios of $(\Sigma,\mu_{0})$, the Ridgeless estimator $\widehat{\mu}_{0}$ has a smaller (limiting) prediction risk compared to any Ridge estimator $\widehat{\mu}_{\eta}$ with $\eta>0$, and therefore in such scenarios interpolation becomes optimal for the task of prediction, cf. [KLS20, HMRT22, TB22]. As a result, the Ridgeless estimator $\widehat{\mu}_{0}$ can be viewed, at least to some extent, as a (substantially) simplified yet theoretically tractable proxy that captures some important features of modern machine learning methods. Moreover, understanding for the prediction risk of $\widehat{\mu}_{0}$ serves as a basis for more complicated interpolation methods in e.g., kernel ridge regression [LR20, BMR21], random features model [MM22, MMM22], neural tangent model [MZ22], among others.

Despite these encouraging progress, there remains limited understanding of the behavior of the Ridgeless estimator $\widehat{\mu}_{0}$ as a statistical estimator. From a statistical perspective, a ‘good’ estimator usually possesses multiple desirable properties, and therefore it is natural to ask whether the Ridgeless estimator $\widehat{\mu}_{0}$ is also ‘good’ in other senses beyond the prediction accuracy. This is particularly relevant, if we aim to consider $\widehat{\mu}_{0}$ also as a ‘good’ estimator that can be applied in broader statistical contexts, rather than viewing it solely as a theoretical proxy designed to provide insights into the mechanisms of modern machine learning methods.

From a different angle, while much of the aforementioned research has focused on identifying favorable scenarios of $(\Sigma,\mu_{0})$ in which prediction via $\widehat{\mu}_{0}$ can be accurate or even optimal, this line of theory does not automatically imply that in the practical statistical applications of overparametrized linear regression, the Ridgeless estimator $\widehat{\mu}_{0}$ remains the optimal choice for the task of prediction for ‘typical’ scenarios of $(\Sigma,\mu_{0})$. The prevailing conventional statistical wisdom, which strongly advocates the use of regularization to strike a balance between bias and variance, may be still at work. From both conceptual and practical standpoints, it is therefore essential to understand the extent to which the optimality phenomenon of interpolation, as observed in certain scenarios of $(\Sigma,\mu_{0})$ in the literature, can be considered ‘generic’ within the context of Ridge regression.

The main goal of this paper is to make a further step in understanding the precise stochastic behavior of the Ridge(less) estimator $\widehat{\mu}_{\eta}$, by developing a high-dimensional distributional characterization in the so-called proportional regime, where $m$ and $n$ is of the same order. As will be clear from below, the distributional characterization of the Ridge(less) estimator $\widehat{\mu}_{\eta}$ offers a detailed understanding of its statistical properties across various statistical tasks, extending beyond the prediction accuracy. Notable examples include a characterization for a general class of weighted $\ell_{q}$ risks ($0<q<\infty$) associated with $\widehat{\mu}_{\eta}$, as well as its capacity for statistical inference in terms of constructing confidence intervals.

In addition, the distributional characterization also leads to new insights on the (sub-)optimality of the Ridgeless estimator $\widehat{\mu}_{0}$, which, interestingly, further reconciles with the conventional statistical wisdom: given any covariance structure $\Sigma$, except for an exponentially small proportion of signal $\mu_{0}$’s within the $\ell_{2}$ norm ball, interpolation is optimal in prediction, estimation and inference, if and only if the noise level $\sigma_{\xi}^{2}=0$. In other words, in the common statistical setting where the noise level is nontrivial $\sigma_{\xi}^{2}>0$, a certain amount of regularization in Ridge regression remains beneficial for ‘most’ signal $\mu_{0}$’s across a number of different statistical tasks in the linear model (1.1)—including the most intensively studied prediction task in the literature.

Of course, in practice, the optimal regularization is unknown and may differ significantly for each statistical task. Surprisingly, our distributional characterization reveals that two widely used adaptive tuning methods, the generalized cross-validation scheme [CW79, GHW79] and the $k$-fold cross-validation scheme [GKKW02, JWHT21]—both of which designed for tuning the prediction risk—are actually simultaneously optimal for all the aforementioned statistical tasks, at least for ‘most’ signal $\mu_{0}$’s. In particular, for ‘most’ signal $\mu_{0}$’s, these two popular adaptive tuning methods automatically lead to optimal prediction, estimation and in-sample risks. Moreover, when combined with the debiased Ridge estimator, they produce the shortest confidence intervals for the coordinates of $\mu_{0}$ with asymptotically valid coverage.

For simplicity of discussion, we will focus here on the overparametrized regime $\phi^{-1}>1$ that is of our main interest.

While our theory (1.6) is strong enough to characterize all $\ell_{q}$ risks of the Ridge(less) estimator $\widehat{\mu}_{\eta}$ in the sense of (1.7), the uniform nature of (1.6) also illuminates novel insights into certain global, qualitative behavior of the most commonly studied $\ell_{2}$ risks for finite samples. To fix notation, we define

• •

  (prediction risk) $R^{\operatorname{\mathsf{pred}}}_{(\Sigma,\mu_{0})}(\eta)\equiv\lVert\Sigma^{1/2}(\widehat{\mu}_{\eta}-\mu_{0})\rVert^{2}$,

• •

  (estimation risk) $R^{\operatorname{\mathsf{est}}}_{(\Sigma,\mu_{0})}(\eta)\equiv\lVert\widehat{\mu}_{\eta}-\mu_{0}\rVert^{2}$,

• •

  (in-sample risk) $R^{\operatorname{\mathsf{in}}}_{(\Sigma,\mu_{0})}(\eta)\equiv n^{-1}\lVert X(\widehat{\mu}_{\eta}-\mu_{0})\rVert^{2}$.

Using our uniform distributional characterization in (1.6), we show in Theorem 3.4 that for ‘most’ $\mu_{0}$’s and all $\#\in\{\operatorname{\mathsf{pred}},\operatorname{\mathsf{est}},\operatorname{\mathsf{in}}\}$, with high probability,

In fact, we prove in Theorem 3.4 a much stronger statement: for ‘most’ $\mu_{0}$’s, the (random) global optimum of $\eta\mapsto R^{\#}_{(\Sigma,\mu_{0})}(\eta)$ for all $\#\in\{\operatorname{\mathsf{pred}},\operatorname{\mathsf{est}},\operatorname{\mathsf{in}}\}$ will be achieved approximately at $\eta_{\ast}=\operatorname{\mathsf{SNR}}_{\mu_{0}}^{-1}$ with high probability. Here $\operatorname{\mathsf{SNR}}_{\mu_{0}}=\lVert\mu_{0}\rVert^{2}/\sigma_{\xi}^{2}$ is the usual notion of signal-to-noise ratio; when $\mu_{0}\neq 0$ and $\sigma_{\xi}^{2}=0$, we shall interpret $\operatorname{\mathsf{SNR}}_{\mu_{0}}^{-1}=0$.

It must be stressed that, for different $\#\in\{\operatorname{\mathsf{pred}},\operatorname{\mathsf{est}},\operatorname{\mathsf{in}}\}$, the empirical risk curves $\eta\mapsto R^{\#}_{(\Sigma,\mu_{0})}(\eta)$ concentrate on genuinely different deterministic counterparts $\eta\mapsto\bar{R}^{\#}_{(\Sigma,\mu_{0})}(\eta)$ with different mathematical expressions (cf. Theorem 3.1). As such, there are no apriori reasons to expect that these risk curves share approximately the same global minimum. Remarkably, as a consequence of the approximate formulae for the deterministic risk curves $\eta\mapsto\bar{R}^{\#}_{(\Sigma,\mu_{0})}(\eta)$ (cf. Theorem 3.2), we show that the curves $\eta\mapsto\bar{R}^{\#}_{(\Sigma,\mu_{0})}(\eta)$ are qualitatively similar, in that they approximately behave locally like a quadratic function centered around $\eta_{\ast}=\operatorname{\mathsf{SNR}}_{\mu_{0}}^{-1}$ (cf. Proposition 3.3), at least for ‘most’ signal $\mu_{0}$’s.

From a broader perspective, the phase transition (1.8) aligns closely with the conventional statistical wisdom: for ‘most’ signal $\mu_{0}$’s, certain amount of regularization is necessary to achieve optimal performance for all the prediction, estimation and in-sample risks, as long as the noise level $\sigma_{\xi}^{2}>0$.

The phase transition in (1.8) naturally raises the question of how one can choose the optimal regularization in a data-driven manner. Here we study two widely used adaptive tuning methods, namely,

1. (1)

   the generalized cross-validation scheme $\widehat{\eta}^{\operatorname{\mathsf{GCV}}}$, and

2. (2)

   the $k$-fold cross-validation scheme $\widehat{\eta}^{\operatorname{\mathsf{CV}}}$.

The readers are referred to (4.3) and (4.5) for precise definitions of $\widehat{\eta}^{\operatorname{\mathsf{GCV}}},\widehat{\eta}^{\operatorname{\mathsf{CV}}}$ in the context of Ridge regression. Both methods have a long history in the literature; see, e.g., [Sto74, Sto77, CW79, GHW79, Li85, Li86, Li87, DvdL05] for some historical references. In essence, both methods $\widehat{\eta}^{\operatorname{\mathsf{GCV}}},\widehat{\eta}^{\operatorname{\mathsf{CV}}}$ are designed to estimate the prediction risk, so it is natural to expect that they perform well for the task of prediction. Rigorous theoretical justifications along this line in Ridge regression in high dimensional settings can be found in, e.g., [LD19, PWRT21, HMRT22].

Interestingly, the insight from (1.8) suggests a far broader utility of these adaptive tuning methods. As all the empirical risk curves $\eta\mapsto R^{\#}_{(\Sigma,\mu_{0})}(\eta)$ are approximately minimized at the same point $\eta_{\ast}=\operatorname{\mathsf{SNR}}_{\mu_{0}}^{-1}$ with high probability, it is plausible to conjecture that the aforementioned cross-validation methods $\widehat{\eta}^{\operatorname{\mathsf{GCV}}},\widehat{\eta}^{\operatorname{\mathsf{CV}}}$ could also yield optimal performance for estimation and in-sample risks, at least for ‘most’ signal $\mu_{0}$’s. We show in Theorems 4.2 and 4.3 that this is indeed the case: for ‘most’ signal $\mu_{0}$’s and all $\#\in\{\operatorname{\mathsf{pred}},\operatorname{\mathsf{est}},\operatorname{\mathsf{in}}\}$, with high probability,

Even more surprisingly, the optimality of $\widehat{\eta}^{\operatorname{\mathsf{GCV}}},\widehat{\eta}^{\operatorname{\mathsf{CV}}}$ extends to the much more difficult task of statistical inference. In fact, we show in Theorem 4.4 that in the so-called debiased Ridge scheme, these two adaptive tuning methods $\widehat{\eta}^{\operatorname{\mathsf{GCV}}},\widehat{\eta}^{\operatorname{\mathsf{CV}}}$ yield an asymptotically valid construction of confidence intervals for the coordinates of $\mu_{0}$ with the shortest possible length.

To the best of our knowledge, theoretical optimality properties for the cross-validation schemes beyond the realm of prediction accuracy has not been observed in the literature, either for Ridge regression or for other regularized regression estimators. Our findings here in the context of Ridge regression can therefore be viewed as a first step in understanding the potential broader merits of cross validation schemes for a larger array of statistical inference problems.

As previously mentioned, a significant body of recent research (cited above) has concentrated on various facets of the precise asymptotics of the prediction risk for the Ridgeless estimator $\widehat{\mu}_{0}$. These studies extensively utilize the explicit form of (1.3), and rely almost exclusively on techniques from random matrix theory (RMT) to relate the behavior of the bias and variance terms in the prediction risk and the spectral properties of the sample covariance.

Here, as (1.6) lacks a direct connection to the spectrum of the sample covariance, we adopt a different, two-step strategy for its proof:

1. (1)

   In the first step, we establish (1.6) under a Gaussian design $X$ via the so-called convex Gaussian min-max theorem (CGMT) approach [Gor85, Gor88].

2. (2)

   In the second step, we prove universality that lifts the Gaussianity $X$ via leveraging the recently developed comparison inequalities in [HS22].

Under a Gaussian design, the proof method for establishing distributional properties of regularized regression estimators via the CGMT has been executed for the closely related Lasso estimator with strict non-vanishing regularization, cf. [MM21, CMW22]¹¹1More literature on other statistical applications of the CGMT can be found in Section 6.. Here in our setting, the major technical hurdle is to handle the vanishing strong convexity in the optimization problem (1.2) as $\eta\downarrow 0$. We overcome this technical issue by establishing uniform localization of Gordon’s min-max optimization, valid down to $\eta=0$. The localization property, in a certain sense, allows us to conclude that the distributional properties of $\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}}(\gamma_{\eta,\ast};\tau_{\eta,\ast})$ are stable as $\eta\downarrow 0$. The readers are referred to Section 6.3 for a proof outline.

For the universality problem, the key step of our arguments is to reduce, in a proper sense, the difficult problem on the universality of $\widehat{\mu}_{0}$ to an easier problem on the universality of $\widehat{\mu}_{\eta}$ for small $\eta>0$, so that the comparison inequalities in [HS22] can be applied. This reduction is achieved by proving that (i) $\widehat{\mu}_{\eta}$ and related quantities are uniformly delocalized down to the interpolation regime $\eta=0$, and (ii) both the primal and Gordon optimization problems are ‘stable’ in suitable senses as $\eta\downarrow 0$. A more detailed proof outline is contained in Section 6.4.

The rest of the paper is organized as follows. In Section 2, we present our main results on the distributional characterizations (1.6) of the Ridge(less) estimator $\widehat{\mu}_{\eta}$, and the approximate $\ell_{q}$ risk formulae (1.7). In Section 3, we provide a number of approximate $\ell_{2}$ risk formulae via RMT, and establish the phase transition (1.8). In Section 4, we give a formal validation for the two cross validation schemes mentioned above, both in terms of (1.9) and statistical inference via the debiased Ridge estimator. A set of illustrative simulation results are presented in Section 5 to collaborate (some of) the theoretical results. Due to the high technicalities involved in the proof of (1.6), a proof outline will be given in Section 6. All the proof details are then presented in Sections 7-12.

For any positive integer $n$, let $[n]=[1:n]$ denote the set $\{1,\ldots,n\}$. For $a,b\in\mathbb{R}$, $a\vee b\equiv\max\{a,b\}$ and $a\wedge b\equiv\min\{a,b\}$. For $a\in\mathbb{R}$, let $a_{\pm}\equiv(\pm a)\vee 0$. For $x\in\mathbb{R}^{n}$, let $\lVert x\rVert_{p}$ denote its $p$-norm $(0\leq p\leq\infty)$, and $B_{n;p}(R)\equiv\{x\in\mathbb{R}^{n}:\lVert x\rVert_{p}\leq R\}$. We simply write $\lVert x\rVert\equiv\lVert x\rVert_{2}$ and $B_{n}(R)\equiv B_{n;2}(R)$. For a matrix $M\in\mathbb{R}^{m\times n}$, let $\lVert M\rVert_{\operatorname{op}},\lVert M\rVert_{F}$ denote the spectral and Frobenius norm of $M$, respectively. $I_{n}$ is reserved for an $n\times n$ identity matrix, written simply as $I$ (in the proofs) if no confusion arises. For a square matrix $M\in\mathbb{R}^{n\times n}$, we let $\mathrm{diag}(M)\equiv(M_{ii})_{i=1}^{n}\in\mathbb{R}^{n}$.

We use $C_{x}$ to denote a generic constant that depends only on $x$, whose numeric value may change from line to line unless otherwise specified. $a\lesssim_{x}b$ and $a\gtrsim_{x}b$ mean $a\leq C_{x}b$ and $a\geq C_{x}b$, abbreviated as $a=\mathcal{O}_{x}(b),a=\Omega_{x}(b)$ respectively; $a\asymp_{x}b$ means $a\lesssim_{x}b$ and $a\gtrsim_{x}b$, abbreviated as $a=\Theta_{x}(b)$. $\mathcal{O}$ and $\mathfrak{o}$ (resp. $\mathcal{O}_{\mathbf{P}}$ and $\mathfrak{o}_{\mathbf{P}}$) denote the usual big and small O notation (resp. in probability). For a random variable $X$, we use $\operatorname{\mathbb{P}}_{X},\operatorname{\mathbb{E}}_{X}$ (resp. $\operatorname{\mathbb{P}}^{X},\operatorname{\mathbb{E}}^{X}$) to indicate that the probability and expectation are taken with respect to $X$ (resp. conditional on $X$).

For a measurable map $f:\mathbb{R}^{n}\to\mathbb{R}$, let $\lVert f\rVert_{\mathrm{Lip}}\equiv\sup_{x\neq y}\lvert f(x)-f(y)\rvert/\lVert x-y\rVert$. $f$ is called $L$-Lipschitz iff $\lVert f\rVert_{\mathrm{Lip}}\leq L$. For a proper, closed convex function $f$ defined on $\mathbb{R}^{n}$, its Moreau envelope $\mathsf{e}_{f}(\cdot;\tau)$ and proximal operator $\operatorname{\mathsf{prox}}_{f}(\cdot;\tau)$ for any $\tau>0$ are defined by $\mathsf{e}_{f}(x;\tau)\equiv\min_{z\in\mathbb{R}^{n}}\big{\{}\frac{1}{2\tau}\lVert x-z\rVert^{2}+f(z)\big{\}}$ and $\operatorname{\mathsf{prox}}_{f}(x;\tau)\equiv\operatorname*{arg\,min\,}_{z\in\mathbb{R}^{n}}\big{\{}\frac{1}{2\tau}\lVert x-z\rVert^{2}+f(z)\big{\}}$.

Throughout this paper, for an invertible covariance matrix $\Sigma\in\mathbb{R}^{n\times n}$, we write $\mathcal{H}_{\Sigma}\equiv\operatorname{tr}(\Sigma^{-1})/n$ as the harmonic mean of the eigenvalues of $\Sigma$.

For $K>1$, let

This notation will be used throughout the paper for uniform-in-$\eta$ statements. In particular, in the overparametrized regime $\phi^{-1}\geq 1+1/K$, we have $\Xi_{K}=[0,K]$.

Next, for $\gamma,\tau\geq 0$, recall $\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}}(\gamma;\tau)$ in (1.2.1), and we define its associated estimation error $\operatorname{\mathsf{err}}_{(\Sigma,\mu_{0})}(\gamma;\tau)$ and the degrees-of-freedom $\operatorname{\mathsf{dof}}_{(\Sigma,\mu_{0})}(\gamma;\tau)$ as

We note that the $\operatorname{\mathsf{dof}}_{(\Sigma,\mu_{0})}(\gamma;\tau)$ defined above is naturally related to the usual notion of degrees-of-freedom (cf. [Ste81, Efr04]) for $\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}}(\gamma;\tau)$, in the sense that $\mathrm{df}\big{(}\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}}(\gamma;\tau)\big{)}\equiv\sum_{j=1}^{n}\frac{1}{\gamma^{2}/n}\operatorname{Cov}\big{(}(\Sigma^{1/2}\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}})_{j},y_{(\Sigma,\mu_{0}),j}^{\operatorname{\mathsf{seq}}}\big{)}=\frac{n}{\gamma^{2}}\operatorname{\mathbb{E}}\operatorname{\mathsf{dof}}_{(\Sigma,\mu_{0})}(\gamma;\tau)$.

Here ‘uniform sub-gaussianity’ means $\sup_{i\in[m],j\in[n]}\lVert Z_{ij}\rVert_{\psi_{2}}\leq C$ for some universal $C>0$, where $\psi_{2}$ is the Orlicz 2-norm (cf. [vdVW96, Section 2.2, pp. 95]).

We shall often write the Gaussian design as $Z=G$, where $G\in\mathbb{R}^{m\times n}$ consists of i.i.d. $\mathcal{N}(0,1)$ entries.

The requirement on the noise level $\sigma_{\xi}^{2}$ will be specified in concrete results below.

Fix $\eta\geq 0$. Consider the following fixed point equation in $(\gamma,\tau)$:

We first establish some qualitative properties for the solution of (2.1).

The above proposition combines parts of Propositions 8.1 and 11.2.

As an important qualitative consequence of (2), under the condition $\lVert\Sigma\rVert_{\operatorname{op}}\vee\mathcal{H}_{\Sigma}\leq K$, the effective regularization $\eta\mapsto\tau_{\eta,\ast}$ is a strictly increasing and concave function of $\eta$. Moreover, in the overparametrized regime $\phi^{-1}>1$, the quantity $\tau_{0,\ast}$—also known as ‘implicit regularization’ in the literature [BLLT20, BMR21, CM22, HMRT22, TB22]—is strictly bounded away from zero.

The claim in (3) offers a useful approximate representation of the effective noise $\gamma_{\eta,\ast}^{2}$ in terms of the original noise $\sigma_{\xi}^{2}$, the effective regularization $\tau_{\eta,\ast}$ and the signal energy $\lVert\mu_{0}\rVert$ without explicitly dependence of $\Sigma$. This representation will prove useful in the approximate $\ell_{q}$ risk formulae in Theorem 2.5, as well as in understanding some qualitative aspects of the risk curves in Section 3 ahead.

The second equation of (2.1) has a natural connection to random matrix theory (RMT). To detail this connection, let $\widehat{\Sigma}\equiv\Sigma^{1/2}G^{\top}G\Sigma^{1/2}/m\in\mathbb{R}^{n\times n}$ and $\check{\Sigma}\equiv G\Sigma G^{\top}/m\in\mathbb{R}^{m\times m}$ be the sample covariance matrix and its dimension flipped, companion matrix. For $z\in\mathbb{C}^{+}\equiv\{z\in\mathbb{C}:\Im z>0\}$, let $\mathfrak{m}_{n}(z)\equiv m^{-1}\operatorname{tr}\big{(}\check{\Sigma}-zI_{m}\big{)}^{-1}$ and $\mathfrak{m}(z)$ be the Stieltjes transforms of the empirical spectral distribution and the asymptotic eigenvalue density (cf. [KY17, Definition 2.3]) of $\check{\Sigma}$, respectively. It is well-known that $\mathfrak{m}(z)$ can be determined uniquely via the fixed point equation

See, e.g., [KY17, Lemma 2.2] for more technical details and historical references. We also note that while the above equation is initially defined for $z\in\mathbb{C}^{+}$, it can be straightforwardly extended to the real axis provided that $z$ lies outside the support of the asymptotic spectrum of $\check{\Sigma}$.

The following proposition provides a precise connection between the effective regularization $\tau_{\eta,\ast}$ defined via the second equation of (2.1), and the Stieltjes transform $\mathfrak{m}$. This connection will prove important in some of the results ahead.

While (2.3) appears somewhat purely algebraic, it actually admits a natural statistical interpretation. Suppose $\xi$ is also Gaussian. We may then compute

Now comparing the above display with (2.3), we arrive at the following intriguing equivalence between the averaged law in RMT, and the proximity of $\widehat{\mu}_{\eta}$ and $\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}}(\gamma_{\eta,\ast};\tau_{\eta,\ast})$ in terms of “degrees-of-freedom”:

Below we will show the above proximity of $\widehat{\mu}_{\eta}$ and $\widehat{\mu}_{(\Sigma,\mu_{0})}^{\operatorname{\mathsf{seq}}}(\gamma_{\eta,\ast};\tau_{\eta,\ast})$ can be taken as far as the distributions of the two themselves.

In addition to $\widehat{\mu}_{\eta}$, we will also consider the distribution of the (scaled) residual $\widehat{r}_{\eta}$, defined by

We define the ‘population’ version of $\widehat{r}_{\eta}$ as

Here $h\sim\mathcal{N}(0,I_{m})$ is independent of $\xi$.

We are now in a position to state our main results on the distributional results for the Ridge(less) estimator $\widehat{\mu}_{\eta}$ and the residual $\widehat{r}_{\eta}$.

First we work under the Gaussian design $Z=G$, and we write $\widehat{\mu}_{\eta}=\widehat{\mu}_{\eta;G},\widehat{r}_{\eta}=\widehat{r}_{\eta;G}$. Recall $\mathcal{H}_{\Sigma}=\operatorname{tr}(\Sigma^{-1})/n$.

The choice $\mu_{0}\in B_{n}(1)$ is made merely for simplicity of presentation; it can be replaced by $\mu_{0}\in B_{n}(R)$ with another constant $C$ that depends further on $R$. The assumption $\mathcal{H}_{\Sigma}\lesssim 1$ is quite common in the literature of Ridge(less) regression; see, e.g., [BMR21, Assumption 4.12] or a slight variant in [MRSY23, Assumption 1]. The major assumption in the above theorem is the Gaussianity on the design $X$. This may be lifted at the cost of a set of slightly stronger conditions.

Theorems 2.3 and 2.4 are proved in Section 9 and Section 10, respectively. Due to the high technicalities in the proof, a sketch is outlined in Section 6.

We mention two particular important features on the theorems above:

1. (1)

   The distributional characterizations for $\widehat{\mu}_{\eta}$ in both theorems above are uniformly valid down to the interpolation regime $\eta=0$ for $\phi^{-1}>1$. This uniform control will play a crucial role in our non-asymptotic analysis of cross-validation methods to be studied in Section 4 ahead.

2. (2)

   The distribution of the residual $\widehat{r}_{\eta}$ in (2) is formulated conditional on the noise $\xi$. A fundamental reason for adopting this formulation is that the distribution of $\widehat{r}_{\eta}$ is not universal with respect to the law of $\xi$. In other words, one cannot simply assume Gaussianity of $\xi$ in Theorem 2.3 in hope of proving universality of $\widehat{r}_{\eta}$ in Theorem 2.4.

In the context of distributional characterizations for regularized regression estimators in the proportional regime, results in similar vein to Theorem 2.3 have been obtained in the closely related Lasso setting for isotropic $\Sigma=I_{n}$ in [MM21], and for general $\Sigma$ in [CMW22], both under Gaussian designs and with strictly non-vanishing regularization. A substantially simpler, isotropic ($\Sigma=I_{n}$) version of Theorem 2.4 is obtained in [HS22] that holds pointwise in non-vanishing regularization level $\eta>0$. As will be clear from the proof sketch in Section 6, in addition to the complications due to the implicit nature of the solution to the fixed point equation (2.1) for general $\Sigma$, the major difficulty in proving Theorems 2.3 and 2.4 rests in handling the singularity of the optimization problem (1.2) as $\eta\downarrow 0$.

As a demonstration of the analytic power of the above Theorems 2.3 and 2.4, we compute below the weighted $\ell_{q}$ risk $\lVert\mathsf{A}(\widehat{\mu}_{\eta}-\mu_{0})\rVert_{q}$ for a well-behaved matrix $\mathsf{A}\in\mathbb{R}^{n\times n}$ and $q\in(0,\infty)$.