Asymptotics of Ridge (less) Regression under General Source Condition

Introduction

Understanding the generalisation properties of overparameterized model is a key question in machine learning, recently popularized by the study of neural networks with millions and even billions of parameters. These models perform well in practice despite perfectly fitting (interpolating) the data, a property that seems at odds with classical statistical theory . This observation has lead to the investigation of the generalisation performance of methods that achieve zero training error (interpolators) and, in the context of linear least squares, the unique least norm solution to which gradient descent converges . Overparameterized linear models, where the number of variables exceed the number of points, are arguably the simplest and most natural setting where interpolation can be studied. Moreover, in some specific regimes, neural networks can be approximated by suitable linear models .

The learning curve (test error versus model capacity) for interpolators has been shown to possibly exhibit a characteristic “Double Descent” shape, where the test error decreases after peaking at an “interpolating” threshold, that is, the model capacity required to interpolate the data. The regime beyond this threshold naturally captures the settings of neural networks , and thus, has motivated its investigation . Indeed, for least squares regression, sharp characterisations double descent have been obtained for the least norm interpolating solution in the case of isotropic or auto-regressive covariates and random features .

For least squares regression the structure of the features and data can naturally influence performance. Within kernel regression (or inverse problems), for instance, it is often assumed that the parameter of interest is regular with respect to a given basis so as to to ensure a well-posed problem . Meanwhile for neural networks, inductive biases can be encoded in the network architecture e.g. convolution layers for image classification . In each case, the problem is made easier by leveraging (through model design) that data encountered in practice exhibits lower dimensional structure owing to, for example, a set of simple physical laws governing the data generation. In contrast, the least squares models investigated beyond the interpolation threshold have focused on cases where the true regression parameter is isotropic , which is a single instance in the range possible of alignements between the parameter and population covariance. This has left open the natural questions of whether additional structure within the data generating distribution can be responsible for determining when interpolating is optimal.

In this work we investigate the performance of ridge regression, and its ridgeless limit, in a high dimensional asymptotic regime with a non-isotropic parameter. We show that one can naturally reduce to a parameter sampled from a prior that, in short, encodes how the signal strength is distributed across the principal components of the covariates. This structure has long been recognized as relevant in the statistics literature , and is analogous to standard smoothness condition used within kernel regression and inverse problems, see e.g. .

Specifically, a prior function encodes the parameter’s norm when it is projected onto eigenspaces of the covariates population covariance. Thus, it represents how aligned the ground truth is to the principle components in the data. When considering the expected test error of ridge regression, this assumption can then encode any deterministic parameter (Proposition 1). Following the classic name in inverse problems, we call these assumptions source conditions.

Given this assumption, we then study the test error of ridge regression in a high-dimensional asymptotic regime when the number of samples and ambient dimension go to infinity in proportion to one another. The limits of resulting quantities are then characterised by utilising tools from asymptotic Random Matrix Theory , with results specifically developed to characterise the influence of the prior function. This provides a natural and intuitive framework for studying the limiting test error of ridge regression, characterised by the signal to noise ratio, regularisation, overparmeterisation, and now, the structure of the regression parameter as encoded by the source condition.

We then illustrate our general framework and results in a simplified setting that highlights the role of model misspecification and its effect on prediction error and regularisation. Specifically, we consider a population covariance with two types of eigenvectors: strong features, associated with a common large eigenvalue (hence favored by the ridge estimator), as well as weak features, with a common smaller eigenvalue. This model is an idealization of a realistic structure for distributions, with some parts of the signal (associated for instance to high smoothness, or low-frequency components) easier to estimate than other, higher-frequency components. The use of source conditions allows to study situations where the true coefficients are either more or less aligned with the principal components, than implicitly postulated by the ridge estimator, a form of model misspecification which affects predictive performance. This encodes the difficulty of the problem, and allows to distinguish between “easy” and “hard” learning problems. We now summarise this work’s primary contributions.

Asymptotic prediction error under general source condition. An asymptotic characterisation of the test error under a general source condition on the regression parameter is provided. This required characterizing the limit of certain trace quantities, and provides a natural framework for investigating the performance of ridge regression. (Theorem 1)

Interpolating can be optimal even in noisy cases. In the overparameterised regime, we show that interpolation can lead to smaller risk than any positive choice of the regularisation parameter. This occurs in the favorable situation where the regression parameter is larger in high-variance directions of the data, and the signal-to-noise ratio is large enough (but finite). Previously, for least squares regression with isotropic prior, the optimal regularisation choice was zero only in the limit of infinite signal to noise ratio . (Section 3.1)

Our analysis of the strong and weak features model also provides asymptotic characterisations of a number of phenomena recently observed within the literature. That is, augmenting the data by adding noisy co-ordinates performs implicit regularisation and can recover the performance of optimally tuned regression restricted to the strong features . Also, we show an additional peak occurring in the learning curve beyond the interpolation threshold for the ridgeless bias and variance . These insights are presented in Sections 3.2 and 3.3, respectively.

The remainder of this work is organized as follows. Section 1.1 covers the related literature. Section 2 describes the setting, and provides the general theorem. Section 3 formally introduces the strong and weak features model, and presents the aforementioned insights. Section 4 gives the conclusion.

Due to the large number of works investigating interpolating methods as well as double descent, we next focus on works that consider the asymptotic regime.

Random matrix theory has found numerous applications in high-dimensional statistics . In particular, asymptotic random matrix theory has been leveraged to study the predictive performance of ridge regression under a well-specified linear model with an isotropic prior on the parameter, for identity population covariance and then general population covariance . More recently, considered the limiting test error of the least norm predictor under the spiked covariance model where both a subset of eigenvalues and the ratio of dimension to samples diverge to infinity. They show the bias is bounded by the norm of the ground truth projected on the eigenvectors associated to the subset of large eigenvalues. In contrast, our work follows standard assumption in kernel regression or inverse problems literature , by adding structural assumptions on the parameter through the variation of its coefficients along the covariance basis. Finally, we note the works that utilise tools from random matrix theory to characterise the prediction performance of linear estimators in the context of classification.

While interpolating predictors (which perfectly fit training data), are classically expected to be sensitive to noise and exhibit poor out-of-sample performance, empirical observations about the behaviour of artificial neural networks challenged this received wisdom. This surprising phenomenon, where interpolators can generalize, has first been shown for some local averaging estimators , kernel “ridgeless” regression , and linear regression, where characterised the variance of the ridgeless estimator up to universal constants. A “double descent” phenomenon for interpolating predictors, where test error can decrease past the interpolation threshold, has been suggested by .

This double descent curve has motivated a number of works established in the context of asymptotic least squares . The work considers either isotropic or auto-regressive features, while consider Random Features constructed from a non-linear functional applied to the product of isotropic covariates and a random matrix. In the data is assumed to be generated with an isotropic ground truth with some model mis-specification. The works considers recovery guarantees under sparsity assumptions on the parameter, with showing a peak in the test error when the number of samples equals the sparsity of the true predictor. The work considers recovery properties of interpolators in the non-asymptotic regime. In contrast to these works, we consider structural assumption on the ground truth in terms of the population covariance that directly follow from standard smoothness conditions in the kernel regression/ inverse problem literature.

The work gave empirical evidence showing additional peaks in the test error can occur beyond the interpolation threshold when the covariance and ground truth parameter are misaligned. These empirical observations are verified by the theory in this paper. Along these lines, we also note the concurrent work which shows a variety of different learning curves are possible for interpolating least squares regression when the sample size is fixed and dimension of the problem is varied.

We now review independent work, which appeared in parallel to or since the first version of this paper. The works also considers the asymptotic prediction performance of ridge regression with prior assumptions on the parameter. Similar to us, shows that interpolating is optimal when the parameter is sufficiently “aligned” to the population covariance and the signal to noise ratio is large. Our technical formulations are formally different but related: they express the alignment between the parameter and the population covariance in terms of the projections of $\beta$ on the eigenvectors of $\Sigma$, whereas we encode it through the source function $\Phi$; the correspondence between the two formulations is obtained through Proposition 1. They also include additional study of the sign of optimal ridge penalty. Meanwhile, has been recently updated to include refined non-asymptotic results that build upon both our work and , also accounting for the structure of the regression parameter along principal directions. They derived a general non-asymptotic bound, controlling the difference between the finite-sample risk and its high-dimensional limit.

Dense Regression with General Source Condition

In this section we formally introduce the setting as well as the main theorem. Section 2.1 introduces the linear regression setting. Section 2.2 shows the prior assumption we consider can encapsulate a general ground truth predictor. Section 2.3 introduces the functionals that arise from asymptotic random matrix theory. Section 2.4 presents the main theorem.

We start by introducing the linear regression setting and the general source condition.

For estimators linear in $Y$ (such as ridge regression), the expected risk only depends on the first two moments of the prior on $\beta^{\star}$, hence one can assume a Gaussian prior $\beta^{\star}\sim\mathcal{N}(0,r^{2}\Phi(\Sigma)/d)$. Under prior (3), $\Phi(\Sigma)^{-1/2}\beta^{\star}$ has isotropic covariance $I/d$, so that $\mathbf{E}\|\Phi(\Sigma)^{-1/2}\beta^{\star}\|^{2}=1$. This means that the coordinate $\beta_{j}:=\langle\beta^{\star},v_{j}\rangle$ of $\beta^{\star}$ in the $j$-th direction has standard deviation $\sqrt{\Phi(\tau_{j})/d}$. We note that, as $d\to\infty$, $\beta^{\star}$ has a “dense” high-dimensional structure, where the number of its components grows with $d$, while their magnitude decreases proportionally. This prior is an average-case, high-dimensional analogue of the standard source condition considered in inverse problems and nonparametric regression , which describes the behaviour of coefficients of $\beta^{\star}$ along the eigenvector basis of $\Sigma$. In the special case $\Phi(x)=x^{\alpha}$, $\alpha\geq 0$, one has $\mathbf{E}\|\Sigma^{-\alpha/2}\beta^{\star}\|^{2}=r^{2}$. For a Gaussian prior, $\Sigma^{-\alpha/2}\beta^{\star}\sim\mathcal{N}(0,r^{2}I/d)$, which is rotation invariant with squared norm distributed as $r^{2}\chi_{d}^{2}/d$ (converging to $r^{2}$ as $d\to\infty$), hence “close” to the uniform distribution on the sphere of radius $r$. In Section 2.2 we show, when considering the expected test error, that this source assumption can then encode any deterministic ground truth parameter.

The case of a constant function $\Phi(x)\equiv 1$ corresponds to an isotropic prior under the Euclidean norm used for regularisation, and has been studied by . In this case (see Remark 1 below), properly-tuned ridge regression (in terms of $r^{2}$) is optimal in terms of average risk. The influence of $\Phi$ can be understood in terms of the average signal strength in eigen-directions of $\Sigma$. Specifically, let $v_{j}$ be an eigenvector of $\Sigma$, with associated eigenvalue $\tau_{j}$. Then, given $\beta^{\star}$, the signal strength in direction $v_{j}$ (namely, the contribution of this direction to the signal) is $\mathbf{E}_{x}\langle\langle\beta^{\star},v_{j}\rangle v_{j},x\rangle^{2}=\tau_{j}\langle\beta^{\star},v_{j}\rangle^{2}$, and its expectation over $\beta^{\star}$ is $\tau_{j}\Phi(\tau_{j})$. When $\Phi$ is increasing, strength along direction $v_{j}$ decays faster as $\tau_{j}$ decreases, than postulated by the ridge regression penalty. In this sense, the problem is lower-dimensional, and hence “easier” than for constant $\Phi$; likewise, a decreasing $\Phi$ is associated to a slower decay of coefficients, and therefore a “harder”, higher-dimensional problem. While our results do not require this restriction, it is natural to consider functions $\Phi$ such that $\tau\mapsto\tau\Phi(\tau)$ is non-decreasing, so that principal components (with larger eigenvalue) carry more signal on average; otherwise, the norm used by the ridge estimator favours the wrong directions. In this respect, the hardest prior is obtained for $\Phi(\tau)=\tau^{-1}$, corresponding to the isotropic prior in the prediction norm induced by $\Sigma$: for this un-informative prior, all directions have same signal strength. Finally, note that in the standard nonparametric setting of reproducing kernel Hilbert spaces, source conditions are related to smoothness of the regression function .

The best linear (in $Y$) estimator in terms of average risk can be described explicitly. It corresponds to the Bayes-optimal estimator under prior $\mathcal{N}(0,r^{2}\Phi(\Sigma)/d)$ on $\beta^{\star}$, which writes:

This estimator requires knowledge of $\Sigma$ and $r^{2}\Phi$. In the special case of an isotropic prior with $\Phi\equiv 1$, the oracle estimator is the ridge estimator (2) with $\lambda=(\sigma^{2}d)/(r^{2}n)$.

2 Reduction to Source Condition

The equality in Proposition 1 holds for finite samples and deterministic $\beta^{\star}$ (and $\Sigma$), and provides a reduction to the setting of random $\beta^{\star}$ used in remaining sections.

On a technical side, the equality in Proposition 1 holds for the expected test error, while the remaining results within this work align with prior work where expectation is taken with respect to the parameter and noise only (conditionally on covariates $X$) i.e. $\mathbf{E}_{\epsilon,\beta^{\star}}[R(\widehat{\beta}_{\lambda})-R(\beta^{\star})]=\mathbf{E}_{\epsilon,\beta^{\star}}[\|\Sigma^{1/2}(\beta-\beta^{\star})\|_{2}^{2}]$. Note that convergence results on the conditional risk can be integrated under suitable domination assumptions, for instance with positive ridge parameter $\lambda$. In addition, framing our next convergence results in the context of deterministic $\beta^{\star}$ would lead to consider source functions $\Phi_{\beta^{\star},\Sigma}=\Phi_{d}$ depending on the dimension $d$, and converging to a fixed function $\Phi$ in a suitable sense as $d\to\infty$. For the sake of simplicity, we instead work in the setting of random parameter $\beta^{\star}$ with a fixed source function $\Phi$.

On another note, the generalised ridge estimator, which penalises with respect to a general covariance $\|P\beta\|_{2}^{2}$ for a positive definite matrix $P$, reduces after rescaling to standard ridge regression with an appropriate prior and covariate covariance. Namely, the problem instance with prior, penalisation and covariate covariances $(\Pi,\!P,\!\Sigma)$ is equivalent to using $(P^{1/2}\Pi P^{1/2},I,P^{-1/2}\Sigma P^{-1/2})$ with parameterisation $\widetilde{\beta}^{\star}=P^{1/2}\beta^{\star}$, $\widetilde{\beta}=P^{1/2}\beta$ and $\widetilde{X}=XP^{-1/2}$.

3 Random Matrix Theory

Let us now describe the considered asymptotic regime, as well as quantities and notions from random matrix theory that appear in the analysis.

We study the performance of the ridge estimator $\widehat{\beta}_{\lambda}$ under high-dimensional asymptotics , where the number of samples and dimension go to infinity $n,d\rightarrow\infty$ proportionally with $d/n\to\gamma\in(0,\infty)$. This setting enables precise characterisation of the risk, beyond the classical regime where $n\to\infty$ with fixed true distribution.

The ratio $\gamma=d/n$ plays a key role. A value of $\gamma>1$ corresponds to an overparameterised model, with more parameters than samples. Some care is required in interpreting this quantity: indeed, for a fixed sample size $n$, varying $\gamma$ changes $d$ and hence the underlying distribution. Hence, $\gamma$ should not be interpreted as a degree of overparmeterisation. Rather, it quantifies the sample size relatively to the dimension of the problem.

which is the limit of the trace quantity d^{-1}\operatorname{Tr}\big{(}\Phi(\Sigma)(\frac{X^{\top}X}{n}-zI)^{-1}\big{)} .

4 Main Theorem: Asymptotic Risk under General Source Condition

Let us now state the main theorem of this work, which provides the limit of the ridge regression risk.

The above theorem characterises the expected test error of the ridge estimator when the sample size and dimension go to infinity $n,d\rightarrow\infty$ with $d/n=\gamma\in(0,\infty)$, and $\beta^{\star}$ is distributed as (3). The asymptotic risk in Theorem 1 is characterised by the relative sample size $\gamma$, the limiting spectral distribution $H$, and the source function $\Phi$ (normalising $\sigma^{2}=r^{2}=1$). This provides a general form for studying the asymptotic test error for ridge regression in a dense high-dimensional setting. The source condition affects the limiting bias; to evaluate it we are required to study the limit of the trace quantity d^{-1}\operatorname{Tr}\big{(}\Sigma(\frac{X^{\top}X}{n}-zI)^{-1}\Phi(\Sigma)(\frac{X^{\top}X}{n}-zI)^{-1}\big{)}, which is achieved utilising techniques from both and (key steps in proof of Lemma 2 Appendix B). The variance term in Theorem 1 aligns with that seen previously in , as the structure of $\beta^{\star}$ only influences the bias.

We now give some examples of asymptotic expected risk in Theorem 1 for $3$ different structures of $\beta^{\star}$, namely $\Phi(x)=1$ (isotropic), $\Phi(x)=x$ (easier case) and $\Phi(x)=x^{-1}$ (harder case).

Consider the setting of Theorem 1. If $n,d\rightarrow\infty$ with $\gamma=d/n$, then almost surely

The three choices of source function $\Phi$ in Corollary 1 are cases where the asymptotic bias in Theorem 1 can be expressed in terms of the companion transform and its first derivative. The expression in the case $\Phi(x)\!=\!1$ was previously investigated in , while for $\Phi(x)\!=\!x$ the bias aligns with quantities previously studied in , and thus, can be simply plugged in. For $\Phi(x)\!=\!x^{-1}$, algebraic manipulations similar to the $\Phi(x)\!=\!x$ case allow $\Theta^{\Phi}(z)$ to be simplified. Finally, for $\Phi(x)\!=\!1$ it is clear how the bias and variance can be brought together and simplified yielding optimal regularisation choice $\lambda\!=\!\sigma^{2}\gamma/r^{2}$ , see also Remark 1. As noted in Section 2.1, $\Phi(x)\!=\!x^{-1}$ corresponds to a “harder" case, with no favoured direction, while $\Phi(x)\!=\!x$ corresponds to an “easier” case with faster coefficient decay.

Strong and Weak Features Model

We call elements of the span of rows of $U_{1}$ strong features, as they are associated to the dominant eigenvalue $\rho_{1}$. Similarly, $U_{2}$ is associated to the weak features. The size of $U_{1},U_{2}$ go to infinity $d_{1},d_{2}\!\rightarrow\!\infty$ with the sample size $n\!\rightarrow\!\infty$, with $d_{i}/d\!\rightarrow\!\psi_{i}\in(0,1)$ and thus $\psi_{1}\!+\!\psi_{2}\!=\!1$. The limiting population spectral measure is then atomic $dH(\tau)\!=\!\psi_{1}\delta_{\rho_{1}}\!+\!\psi_{2}\delta_{\rho_{2}}$.

Under the model just introduced, Theorem 1 provides the following asymptotic characterization for the expected test risk as $n,d\rightarrow\infty$

To gain insights into the performance of least squares when data is generated from the strong and weak features model, we now investigate the above limit in the overparameterised setting $\gamma>1$. The insights are summarised in the following sections. Section 3.1 shows that zero regularisation is optimal for easy problems with high signal to noise ratio. Section 3.2 shows how weak features can be used as a form of regularisation similar to ridge regression. Section 3.3 present findings related to the ridgeless bias and variance.

1 Interpolating can be optimal in the presence of noise

Consider the strong and weak features model with $\gamma=2$, $\psi_{1}=\psi_{2}=1/2$ and $\mathbf{E}[\|\beta^{\star}\|_{2}^{2}]=r^{2}$. If

Looking to Figure 1 plots for the performance of optimally tuned ridge regression (Left) and the optimal choice of regularisation parameter (Right) against (a monotonic transform) of the eigenvalue ratio $\rho_{1}/\rho_{2}$, for a coefficient ratios $\phi_{1}\geq\phi_{2}$ have been given. As shown in the right plot of Figure 1, for a fixed distribution of $X$ (characterised by $\psi_{1},\rho_{1},\rho_{2}$) and sample size (characterised by $\gamma$) as the ratio $\phi_{1}/\phi_{2}$ increases the optimal regularisation decreases. Following Corollary 2, if the ratio $\phi_{1}/\phi_{2}$ is large enough, the optimal ridge regularisation parameter $\lambda$ can be , corresponding to ridgeless interpolation. We note that the negative derivative at (Corollary 2) and the right plot of Figure 1, see also .

2 The Special Case of Noisy Weak Features

In this section we consider the special case where weak features are pure noise variables, namely $\phi_{2}=0$, while their dimension is large. Such noisy weak features can be artificially introduced to the dataset, to induce an overparameterised problem. We then refer to this technique as Noisy Feature Regularisation, and note it corresponds to the design matrix augmentation in . Looking to Figure 2, the ridgeless test error is then plotted against the eigenvalue ratio $\rho_{2}/\rho_{1}$ (Left) and the number of weak features with the tuned eigenvalue ratio (Right).

Observe (right plot) as we increase the number of weak features (as encoded by $1/\psi_{1}$), and tune the eigenvalue $\rho_{2}$, the performance converges to optimally tuned ridge regression with the strong features only. The left plot then shows the “regularisation path” as a function of the eigenvalue ratio $\rho_{2}/\rho_{1}$ for some numbers of weak features $1/\psi_{1}$.

3 Ridgeless Bias and Variance

In this section we investigate how the ridgeless bias and variance depend on the ratio of dimension to sample size $\gamma$. Looking to Figure 3 the ridgeless bias and variance is plotted against the ratio of dimension to sample size in the overparameterised regime $\gamma\geq 1$ .

Note an additional peak in the ridgeless bias and variance is observed beyond the interpolation threshold. This has only recently been empirically observed for the test error , as such, these plots now theoretically verify this phenomenon. The location of the peaks naturally depends on the number of strong and weak features as well as the ambient dimension, as denoted by the vertical lines. Specifically, the peak occurs in the ridgeless bias for the “hard” setting when the number of samples and number of strong features are equal $n=d_{1}$. Meanwhile, a peak occurs in the ridgeless variance when the number of samples and strong features equal $n=d_{1}$, and the eigenvalue ratio is large $\rho_{1}>\rho_{2}$. This demonstrates that learning curves beyond the interpolation threshold can have different characteristics due to the interplay between the covariate structure and underlying data. We conjecture this arises due to instabilities of the design matrix Moore-Penrose Pseudo-inverse, akin to the isotropic setting . Since variance matches prior work , the additional peak could be previously derived. Meanwhile the peak in the bias here uses of the source condition, and thus, as far as we aware is not encompassed in prior work.

Conclusion

In this work, we introduced a framework for studying ridge regression in a high-dimensional regime. We characterised the limiting risk of ridge regression in terms of the dimension to sample size ratio, the spectrum of the population covariance and the coefficients of the true regression parameter along the covariance basis. This extends prior work , that considered an isotropic ground truth parameter. Our extension enables the study of “prior misspecification”, where signal strength may decrease faster or slower than postulated by the ridge estimator, and its effect on ideal regularisation.

We instantiated this general framework to a simple structure, with strong and weak features. In this case, we show that “ridgeless” regression with zero regularisation can be optimal among all ridge regression estimators. This occurs when the signal-to-noise ratio is large and when strong features (with large eigenvalue of the covariance matrix) have sufficiently more signal than weak ones. The latter condition corresponds to an “easy” or “lower-dimensional” problem, where ridge tends to over-penalise along strong features. This phenomenon does not occur for isotropic priors, where optimal regularisation is always strictly positive in the presence of noise. Finally, we discussed noisy weak features, which act as a form of regularisation, and concluded by showing additional peaks in ridgeless bias and variance can occur for our model.

Moving forward, it would be natural to consider non-Gaussian covariates. Given universality results in Random Matrix Theory we expect that the results provided here extend to the case of random vectors with independent coordinates (and linear transformations thereof). Other structures for the ground truth and data generating process can be investigated through Theorem 1 by consider different functions $\Phi$ and the population eigenvalue distributions. The tradeoff between prediction and estimation error exhibited by in the isotropic case can be explored with a general source $\Phi$.

Acknowledgments

D.R. is supported by the EPSRC and MRC through the OxWaSP CDT programme (EP/L016710/1). Part of this work has been carried out at the Machine Learning Genoa (MaLGa) center, Universita di Genova (IT). L.R. acknowledges the financial support of the European Research Council (grant SLING 819789), the AFOSR projects FA9550-17-1-0390 and BAA-AFRL-AFOSR-2016-0007 (European Office of Aerospace Research and Development), and the EU H2020-MSCA-RISE project NoMADS - DLV-777826. We would also like to thank the anonymous reviewers for their feedback and suggestions.

References

Appendix A Proofs for Ridge Regression

In this section we provide the calculations associated to ridge regression. Section A.1 provides the proof of Proposition 1. Section A.2 provides the calculation for the oracle estimator presented in remark 1. Section A.3 provides some preliminary calculations related to random matrix theory. Section A.4 gives the proof of Theorem 1. Section A.5 provides the proof of Corollary 1. Section A.6 provides the calculations associated to the strong and weak features model.

Denote $\beta^{\prime}=U\beta$, as well as $X^{\prime}=XU^{-1}$ and $x^{\prime}=U^{-1}x$. Let $\widehat{\beta}_{\lambda}^{\prime}$ the Ridge estimator computed on data $(X^{\prime},Y)$, namely

Then, $y=\langle\beta,x\rangle+\sigma\epsilon=\langle\beta^{\prime},x^{\prime}\rangle+\sigma\epsilon$, hence the best linear predictor of $y$ based on $x^{\prime}$ is $\beta^{\prime}$. In addition, $x^{\prime}$ has distribution $\mathcal{N}(0,U^{-1}\Sigma U)=\mathcal{N}(0,\Sigma)$, where $U^{-1}\Sigma U=\Sigma$ comes from the fact that $U$ is an isometry on the eigenspaces $V_{j}$ of $\Sigma$. This implies that $(X^{\prime},\epsilon)$ has the same distribution as $(X,\epsilon)$, and thus $\mathbf{E}_{\epsilon,X}[\mathcal{E}_{\beta^{\prime}}(\widehat{\beta}^{\prime}_{\lambda})]=\mathbf{E}_{\epsilon,X}[\mathcal{E}_{\beta^{\prime}}(\widehat{\beta}_{\lambda})]$. On the other hand,

(note that $\|\Sigma^{1/2}U\cdot\|^{2}_{2}=\|U\Sigma^{1/2}\cdot\|^{2}_{2}=\|\Sigma^{1/2}\cdot\|^{2}_{2}$ as $U$ commutes with $\Sigma^{1/2}$ and is an isometry), so that $\mathbf{E}_{\epsilon,X}[\mathcal{E}_{\beta^{\prime}}(\widehat{\beta}^{\prime}_{\lambda})]=\mathbf{E}_{\epsilon,X}[\mathcal{E}_{\beta}(\widehat{\beta}_{\lambda})]$. This proves that $\mathbf{E}_{\epsilon,X}[\mathcal{E}_{\beta^{\prime}}(\widehat{\beta}_{\lambda})]=\mathbf{E}_{\epsilon,X}[\mathcal{E}_{\beta}(\widehat{\beta}_{\lambda})]$. ∎

We now turn to the proof of Proposition 1:

Let $V_{1},\dots,V_{k}$ denote the eigenspaces of $\Sigma$, with distinct eigenvalues $\tau_{1}^{\prime}>\dots>\tau_{k}^{\prime}$. Let $U_{1},\dots,U_{k}$ be independent random isometries, where $U_{j}$ is distributed according to the uniform (Haar) measure on the orthogonal group of $V_{j}$. Define $U$ to be the random isometry acting as $U_{j}$ on $V_{j}$, and let $\beta=U\beta^{\star}$ and $\Pi$ its distribution.

Note that $U$ is of the form of Lemma 1, hence $\mathbf{E}_{X,\epsilon}[\mathcal{E}_{U\beta^{\star}}(\widehat{\beta}_{\lambda})]=\mathbf{E}_{X,\epsilon}[\mathcal{E}_{\beta^{\star}}(\widehat{\beta}_{\lambda})]$ and thus

Now, let $\beta_{j}^{\prime}\in V_{j}$ be the orthogonal projection of $\beta^{\star}$ on $V_{j}$, so that $U\beta^{\star}=\sum_{j=1}^{k}U_{j}\beta_{j}^{\prime}$. We have $\mathbf{E}[U\beta^{\star}]=0$ since $\mathbf{E}[U_{j}]=0$ for all $j$. In addition, the distribution of $U_{j}\beta^{\star}$ is invariant by rotation (since $R_{j}U_{j}$ has the same distribution as $U_{j}$ for any fixed rotation $R_{j}$), hence $\mathbf{E}[(U_{j}\beta_{j}^{\prime})(U_{j}\beta_{j}^{\prime})^{\top}]=t_{j}I_{V_{j}}$ (with $I_{V_{j}}$ the identity on $V_{j}$), where letting $d_{j}=\dim(V_{j})$,

hence $t_{j}=\|\beta^{\prime}_{j}\|^{2}_{2}/d_{j}$. In addition, if $j\neq l$, by independence of $U_{j},U_{l}$,

A.2 Proof of Oracle Estimator (Remark 1)

Since the risk is quadratic, the average risk (integrated over the prior) of any estimator linear in $Y$ only depends on the first two moments of the prior, hence one can assume that the prior is Gaussian (namely, $\mathcal{N}(0,r^{2}\Phi(\Sigma)/d)$) without loss of generality. In this case, a standard computation shows that the posterior is $\mathcal{N}(\widetilde{\beta},[{X^{\top}X}+(\sigma^{2}d/r^{2})\Phi(\Sigma)^{-1}]^{-1})$, where $\widetilde{\beta}$ is the estimator defined in (4). Finally, since the risk is quadratic, the Bayes-optimal estimator is the posterior mean, which corresponds to $\widetilde{\beta}$.

A.3 Random Matrix Theory Preliminaries

We now introduce some useful properties of the Stieltjes transform as well as its companion transform. Firstly, we know the companion transform satisfies the Silverstein equation

Meanwhile from from the equality $\gamma(m(z)+1/z)=v(z)+1/z$ we note that we have the following equalities

which we will readily use to simplify/rewrite a number of the limiting functions.

A.4 Proof of Theorem 1

We begin with the decomposition into bias and variance terms following . The difference for the ridge parameter can be denoted

And thus taking expectation with respect to the noise in the observations $\epsilon$

Taking expectation with respect to $\mathbf{E}_{\beta^{\star}}$ we arrive at

It is now a matter of showing the asymptotic almost sure convergence of the following three functionals

The limit of the first trace quantity comes directly from meanwhile the limit of the second trace quantity is proven in . The third trace quantity depends upon the source condition $\Phi$ and computing its limit is one of the main technical contributions of this work. The limits for these objects is summarised within the following Lemma, the proof of which provides the key steps for computing the limit involving the source function.

Under the assumptions of Theorem 1 for any $\lambda>0$ we have almost surely as $n,d\rightarrow\infty$ with $d/n=\gamma$

The result is arrived at by plugging in the above limits and noting from the definition of the Companion Transform $v$ that $1-\gamma(1-\lambda m(-\lambda))=\lambda v(-\lambda)$, $1-\lambda m(-\lambda)=\gamma^{-1}(1-\lambda v(-\lambda))$ and, taking derivatives, $m(-\lambda)-\lambda m^{\prime}(-\lambda)=\gamma^{-1}(v(-\lambda)-\lambda v^{\prime}(-\lambda))$. The proof of Lemma 2, which is the key technical step in the proof of Theorem 1, is provided in Appendix B.

A.5 Proof of Corollary 1

In this section we provide the proof of Corollary 1. It will be broken into three parts associated to the three cases $\Phi(x)=x$, $\Phi(x)=1$ and $\Phi(x)=1/x$.

The purpose of this section is to demonstrate, in the case $\Phi(x)=x$, how the functional $\Theta^{\Phi}(-\lambda)+\lambda\frac{\partial\Theta^{\Phi}(-\lambda)}{\partial\lambda}$ can be written in terms of the Stieltjes Transform $m(z)$. For this particular choice of $\Phi$ the asymptotics were calculated in , see also Lemma 7.9 in . We therefore repeat this calculation for completeness. Now, in this case we have

Following the steps are the start of the proof for Lemma 2.2 in , consider $1+zm(z)$

Picking $z=-\lambda$ and differentiating with respect to $\lambda$ we get

where on the second equality we used (12). Multiplying through by $\lambda^{2}$ then yields the quantity presented.

A.5.2 Case: Φ​(x)=1Φ𝑥1\Phi(x)=1

The functional of interest in this case aligns with that calculated within , which we include below for completeness. In particular we have $\Theta^{\Phi}(-\lambda)=m(-\lambda)$ and as such we get

where on the second equality we used (12). Dividing by $v(-\lambda)^{2}$ as well as adding the asymptotic variance we get, from Theorem 1, the limit as $n,d\rightarrow\infty$

A.5.3 Case: Φ​(x)=1/xΦ𝑥1𝑥\Phi(x)=1/x

The functional in the case $\Phi(x)=1/x$ takes the form

Solving for $\Theta^{\Phi}(z)$ and plugging in the definition of the companion transform $v(z)$ we arrive at

Fixing $z=-\lambda$ the quantity of interest then has the form

which when differentiated with respect to $\lambda$ yields

Multiplying the above by $\lambda$ and adding $\Theta^{\Phi}(-\lambda)$ brings us to

Dividing the above by $v(-\lambda)^{2}$ and adding the limiting variance yields, from Theorem 1, the limit as $n,d\rightarrow\infty$

A.6 Strong and Weak Features Model

This section presents the calculations associated to the strong and weak features model. We begin giving the stationary point equation of the companion transform $v(t)$, after which we explicitly compute the limiting risk with the particular choice of $\Phi(x)$ in this case. Section A.6.1 there after gives explicit form for the companion transform in the ridgeless limit. Section A.6.2 gives the proof of Corollary 2 found within the main body of the manuscript.

We begin by recalling the limiting spectrum of the covariance $\Sigma$ for the two Bulks Model is $dH(\tau)=\psi_{1}\delta_{\rho_{1}}+\psi_{2}\delta_{\rho_{2}}$. Recall we have $\psi_{1}+\psi_{2}=1$ therefore we simply write $\psi_{2}=1-\psi_{1}$. Using the Silverstein equations (11) the companion transform must satisfy

as such given $v(t)$ we can compute the derivative. Rearranging (16) and denoting $v(t)=v$ the companion transform evaluated at $t$ satisfies

This cubic can then be solved computationally for different choices of $t$. In the case of the ridgeless limit $t\rightarrow 0$ in the overparameterised setting $\gamma>1$, the above simplifies to a quadratic which can be solved, as shown in Section A.6.1.

where on the last equality we used (12) to rewrite the above in terms of the companion transform. Plugging in the regularisation parameter $z=-\lambda$ we then get

To the end of computing $\frac{\partial\Theta^{\Phi}(-\lambda)}{\partial\lambda}$, we can differentiate the above to get

as required. The final form for the limiting risk is then

To consider the Ridgeless limit $t\rightarrow 0$ of the companion transform $v(t)$, some care must be taken about which regime $\gamma<1$ or $\gamma>1$ we are in.

Following the proof of Lemma 6.2 in we have in the underparameterised case $\gamma<1$ the limit $\lim_{t\rightarrow 0_{-}}tv(t)=1-\gamma$.

Following the proof of Lemma 6.2 in when $\gamma>1$ we have the limit $\lim_{t\rightarrow 0_{-}}v(t)=v(0)$. From dominated convergence theorem we can take the limit in the Silverstein equation (16) to arrive at the quadratic

Solving for $v$ with the quadratic formula immediately gives

Recall from we have that $v(z)\in\mathcal{S}$, as such we take the sign above which yields a non-negative quantity. Noting we we focus on the regime where $\gamma>1$, we see for the above to be non-negative we require the numerator to be negative, and thus, we take the negative sign.

A.6.2 Proof of Corollary 2

Meanwhile, recall by differentiating both sides of the silverstein equations (11) in $t$ we can get

Therefore, if we differentiate once more we get

and thus, multiplying through by $v(t)^{3}/v^{\prime}(t)^{2}$ and rearranging we arrive at

Furthermore, noting that $1-\gamma\int\frac{\tau^{2}v(z)^{2}}{(1+\tau v(z))^{2}}dH(\tau)=\frac{v^{\prime}(t)}{v(t)^{2}}$ means we get the following equality for the second derivative

Taking $t\rightarrow 0$ and plugging in the defintion of $v(0)$ yields the following, which will be required for the proof

Taking $\lambda\rightarrow 0$ and plugging in $\psi_{1}=\psi_{2}=1/2$, $\psi_{1}\phi_{1}+\psi_{2}\phi_{2}=1$, $\phi_{1}+\phi_{2}=2$ as well as $v(0)$ into $\rho_{i}v(0)$ for $i=1,2$ yields

Let us now plug in the second derivative $v^{\prime\prime}(0)$. In particular, note that we can write

where on the second equality we have used the equality for $\frac{v^{\prime}(0)}{v(0)^{2}}$ from above to note that

Appendix B Proof of Lemma 2

In this section we provide the proof for Lemma 2. We recall that the limits (13) and (14) have been computed previously. In particular, Lemma 2.2 of (the roles of $d,n$ are swapped in their work, and thus, one must swap $\gamma$ with $1/\gamma$) shows

This leaves us to show the limit (15), for which we build upon the techniques as well as .

Multiplying the above on the left by $\Phi(\Sigma)R(z)$, taking the trace and dividing by $d$ yields

where for $i=1,\dots,n$ we have plugged in (19) twice into for $R(z)$ to get

where the error term $\delta=\delta_{1}+\delta_{2}+\delta_{3}+\delta_{4}$ such that

As shown in section B.1 the error terms $|\delta_{1}|,|\delta_{2}|,|\delta_{3}|,|\delta_{4}|\rightarrow 0$ almost surely as $n,d\rightarrow\infty$. It is now a matter of computing the limits of the remaining terms. As discussed previously the limit of $\frac{1}{d}\operatorname{Tr}(\Sigma R(-\lambda))$ is known from . From the same work it is also known that

That leaves us to compute the limit of \frac{1}{d}\operatorname{Tr}\big{(}\Phi(\Sigma)R(-\lambda)^{2}\big{)}. If we are to write f_{d}(\lambda)=\frac{1}{d}\operatorname{Tr}\big{(}\Phi(\Sigma)R(-\lambda)\big{)} then note the derivative with respect to $\lambda$ is f^{\prime}_{d}(\lambda)=-\frac{1}{d}\operatorname{Tr}\big{(}\Phi(\Sigma)R(-\lambda)^{2}\big{)}. We wish to now study the limit of the $f^{\prime}_{d}(\lambda)$ through the limit of $f_{d}(\lambda)$. To do so we will follow the steps in , which will require some definitions and the following theorem.

Let $f_{1},f_{2},\dots$ be analytic on the domain $D$, satisfying $|f_{n}(z)|\leq M$ for every $n$ and $z$ in $D$. Suppose that there is an analytic function $f$ on $D$ such that $f_{n}(z)\rightarrow f(z)$ for all $z\in D$. Then it also holds that $f^{\prime}_{n}(z)\rightarrow f^{\prime}(z)$ for all $z\in D$

for all $\lambda\in\mathcal{S}:=\{u+iv:v\not=0,\text{ or }v=0,u>0\}$. Checking the conditions of Theorem 2 we have that $f_{d}(\lambda)$ is an analytic function of $\lambda$ on $\mathcal{S}$ and is bounded $|f_{d}(\lambda)|\leq\frac{\|\Phi(\Sigma)\|_{2}}{\lambda}$. To apply Theorem 2 it suffices to show that the limit $\Theta^{\Phi}(-\lambda)$ is analytical. To this end we invoke Morera’s theorem which states if

for any closed curve $\gamma$ in the region $\mathcal{S}$ then $\Theta^{\Phi}(-\lambda)$ is analytic. We see this is the case by applying Fubini’s Theorem as follows

and noting that the inner integral is zero from Cauchy Theorem as $\frac{1}{\tau(1-\gamma(1-\lambda m(-\lambda)))+\lambda}$ is an analytical function of $\lambda$ in $\mathcal{S}$ for any $\tau\in[h_{1},h_{2}]$. By Theorem 2 we have that

The final limit (15) is arrived at by considering the limit as $d,n\rightarrow\infty$ of (20). Specifically, with the fact that $\delta\rightarrow 0$, bringing together (21), (22) and (13). Noting that (13) is applied to the square of 1+\frac{1}{n}\operatorname{Tr}\big{(}\Sigma R(-\lambda)\big{)}=1+\gamma\frac{1}{d}\operatorname{Tr}\big{(}\Sigma R(-\lambda)\big{)}\rightarrow\frac{1}{1-\gamma(1-\lambda m(-\lambda))}.

To analyse these quantities we introduce the following concentration inequality from Lemma A.2 of with $\delta=1/3$.

Furthermore, we will use the fact that the maximal eigenvalues are upper bounded

We proceed to show that each of the error $\delta_{1},\delta_{2},\delta_{3},\delta_{4}$ converge to zero almost surely.

Begin with $\delta_{1}$. For $i=1,\dots,n$ by adding and subtracting $\operatorname{Tr}(\Sigma^{1/2}R(z)\Phi(\Sigma)R_{i}(z)\Sigma^{1/2})$ we can decompose

Using (19) and letting $A=\Phi(\Sigma)R_{i}(-\lambda)\Sigma$ we then get

An identical calculation with $A=\Phi(\Sigma)R(-\lambda)\Sigma$ yields the same bound. This then yields with the lower bound $(1+\operatorname{Tr}(\Sigma^{1/2}R(-\lambda)\Sigma^{1/2}))\geq 1$

and as such $\delta_{1}$ goes to zero as $n,d\rightarrow\infty$ so that $d/n\rightarrow\gamma$.

Now consider the term $\delta_{2}$. Note that for two positive numbers $a,b\geq 0$ we have

and as such $|(1+a)^{-2}-(1+b)^{-2}|\leq 2|b-a|$. Using this with $a=\frac{1}{n}\operatorname{Tr}(\Sigma^{1/2}R_{i}(-\lambda)\Sigma^{1/2})$ and $b=\frac{1}{n}\operatorname{Tr}(\Sigma^{1/2}R(-\lambda)\Sigma^{1/2})$ whom are both non-negative, allows us to upper bound

where for the final inequality we used the argument (23) with $A=\Sigma$. Now, since the eigenvalues in the following trace are non-negative we can upper bound

Combining these two facts yields the upper bound

which goes to zero as $n\rightarrow\infty$.

We now proceed to bound $\delta_{3}$ and $\delta_{4}$. With the bound on the trace (24) as well as using the bound $|(1+a)^{-2}-(1+b)^{-2}|\leq 2|b-a|$ we arrive at the bound for $\delta_{3}$

Meanwhile using that $1+\frac{1}{n}Z_{i}\Sigma^{1/2}R_{i}(-\lambda)\Sigma^{1/2}Z_{i}^{\top}\geq 1$ we arrive at the bound for $\delta_{4}$

We now show that \max_{1\leq i\leq n}\Big{|}Z_{i}\Sigma^{1/2}R_{i}(z)\Phi(\Sigma)R_{i}(z)\Sigma^{1/2}Z_{i}^{\top}-\operatorname{Tr}\big{(}\Sigma^{1/2}R_{i}(z)\Phi(\Sigma)R_{i}(z)\Sigma^{1/2}\big{)}\Big{|} converges to zero almost surely. Observe since we have the upper bound on the largest eigenvalue we have using Lemma 3 as well as union bound for $1\leq i\leq n$ we have for $0<t<\frac{\|\Sigma\|_{2}\|\Phi(\Sigma)\|_{2}}{\lambda^{2}}$

Let V_{n,d}:=\max_{1\leq i\leq n}\frac{1}{d}\Big{|}Z_{i}\Sigma^{1/2}R_{i}(z)\Phi(\Sigma)R_{i}(z)\Sigma^{1/2}Z_{i}^{\top}-\operatorname{Tr}\big{(}\Sigma^{1/2}R_{i}(z)\Phi(\Sigma)R_{i}(z)\Sigma^{1/2}\big{)}\Big{|} and, for any $t>0$, let $E_{n,d}(t)$ denote the event $\{V_{n,d}\geq t\}$ where $d=d_{n}$. Then, if $d=d_{n}$ satisfies $d_{n}/n\to\infty$, \mathbf{P}(E_{n,d})\leq 2n\exp\Big{\{}-\frac{dt^{2}\lambda^{4}}{6\|\Sigma\|_{2}^{2}\|\Phi(\Sigma)\|^{2}}\Big{\}}\leq 2n\exp\Big{\{}-\frac{\gamma nt^{2}\lambda^{4}}{12\|\Sigma\|_{2}^{2}\|\Phi(\Sigma)\|^{2}}\Big{\}} where the last inequality for $n$ large enough that $d/n\geq\gamma/2$. Hence,

so that, by the Borel-Cantelli lemma, almost surely, $V_{n,d}\geq t$ only holds for a finite number of values of $n$. This implies that, almost surely, $\limsup_{n\to\infty}V_{n,d}\leq t$. Note that this is true for every $t>0$; letting $t=1/k$ and taking a union bound over $k\geq 1$ shows that $\limsup_{n\to\infty}V_{n,d}=0$ almost surely, i.e. $V_{n,d}\to 0$ almost surely.