Benign overfitting in ridge regression

Introduction

The bias-variance tradeoff is well known in statistics and machine learning. The classical theory suggests that large models overfit the data and that one needs significant regularization to make them generalize. This intuition is, however, in contrast with the empirical study of modern machine learning techniques. It was repeatedly observed that even models with enough capacity to exactly interpolate the data can generalize with little regularization, or no regularization at all (Belkin et al., 2019a; Zhang et al., 2016). In some cases, the best value of the regularizer can be zero (Liang and Rakhlin, 2018) or even negative (Kobak et al., 2020) for such models.

The aim of this paper is to provide a theoretical understanding of these phenomena, and to do that we consider one of the simplest settings in which they can be observed—ridge regression in dimension $p$ with $n<p$ i.i.d. noisy observations. Despite being a classical statistical methodology, ridge regression and its ridgeless limit are still not completely understood in such a regime: when $n<p$ classical theory suggests that the regularization parameter should be large enough to provide additional capacity control (see, e.g., Hsu et al. (2014) and references therein). The basis of our work was set by Bartlett et al. (2020), who studied the variance term for ridgeless regression with $n<p$ under the additional assumption that the data vectors have independent components. The main discovery of their work is that the variance term can be small if and only if there exists $k^{*}\ll n$ such that if one removes the first $k^{*}$ largest eigenvalues of the covariance operator, the remaining tail of the sequence of eigenvalues has large effective rank compared to $n$. In our work we start afresh and use the same separation of eigendirections from the very beginning, which allows us to substitute the independence assumption by a weaker assumption on the condition number of the Gram matrix of the tails of the data vectors. Moreover, we show how the same separation of the eigenvalues gives tight bounds for the bias term too. Finally, by virtue of algebra, our argument extends very easily to the setting of ridge regression, which allows for comparison with the above mentioned classical results and investigation of the case when the regularization is even less than zero. We show that we extend (with different constants) the results of Hsu et al. (2014) to a larger range of regularization parameters, and give general conditions under which negative regularization is optimal and can provide arbitrarily high multiplicative gain in excess risk.

The structure of the paper is the following: in Section 1.2, we provide an overview of the field of overparameterized ridge regression. We postpone a more technical overview to Section 9, where we also explain how our paper relates to other works. We start the presentation of our results with introducing the setting of ridge regression in Section 2. After that, we use Section 3 to introduce the separation of eigendirections and define the relevant important objects: Subsection 3.1 shows two simple sketches aimed at building up intuition, Subsection 3.2 explains the results of Bartlett et al. (2020) in terms of that intuition and Subsection 3.3 explains how our work completes the story. The aim of this discussion is to elucidate the meaning behind the rigorous assumptions and results that we show in Section 4. Then Section 5 provides a more technical discussion of the main assumption. Section 6 provides an outline of the proof and explains where it uses the assumption that the data is sub-Gaussian. In Section 7 we note that as a side product of the proof an alternative form of the main bound arises, which makes it convenient to compare our bounds to the results of other papers. In Section 8, we derive the sufficient conditions for optimality of negative regularization. Finally, we conclude the paper with Section 10.

2 Related work

Motivated by the empirical success of overparametrized models, there has recently been a flurry of work aimed at understanding theoretically whether the corresponding effects can be seen in overparametrized linear regression; see, e.g., (Liang et al., 2019; Muthukumar et al., 2019; Belkin et al., 2019b; Bibas et al., 2019; Nakkiran, 2019; Xu and Hsu, 2019; Zhou et al., 2021; Negrea et al., 2020) and other references in this section.

The results that aim at characterizing the generalization performance of linear methods can be split roughly into three categories. The first category is results that give exact expressions of the excess risk in the asymptotic setting with ambient dimension and the number of data points going to infinity, while their ratio goes to a constant, and the spectral density of the covariance operator converges weakly to some limiting distribution (Dobriban and Wager, 2015; Hastie et al., 2019; Wu and Xu, 2020; Richards et al., 2020).

The second category is results that make strong assumptions on the distribution of data (e.g., that data vectors have i.i.d. components or come from a uniform distribution on a sphere) and derive bounds on excess risk of linear regression with some specific features, or kernel regression with a kernel that has some specific properties (Montanari and Zhong, 2020; Ghorbani et al., 2020b; Mei and Montanari, 2019; Ghorbani et al., 2020a; Liang et al., 2020). Some of these results are also asymptotic, and some are non-asymptotic.

The third category is results that prove non-asymptotic bounds depending on the arbitrary structure of the covariance of the data. This is the category to which this paper belongs. We already mentioned the work of Bartlett et al. (2020). The other works in this category are (Kobak et al., 2020), (Chinot and Lerasle, 2021), (Dereziński et al., 2019) and (Dereziński et al., 2020).

We provide more detailed comparison and discuss some technical aspects in Section 9.

There have been many related works since the arXiv version of this paper (Tsigler and Bartlett, 2020) was posted (Mei et al., 2021a, b; Ghosh et al., 2021; Misiakiewicz and Mei, 2021; Bartlett et al., 2021; Celentano et al., 2021; Muthukumar et al., 2021; Narang et al., 2021; McRae et al., 2021; Shamir, 2022; Koehler et al., 2021; Bunea et al., 2022) etc. Hastie et al. (2020) obtained a finite sample version of the asymptotic results of the old version of their paper (Hastie et al., 2019). In Section 7.3 we provide an explicit comparison with our results. More recently, Mei et al. (2021a) obtained generalization bounds for kernel ridge regression under similar assumptions to those we consider here (see their Assumption 1). Koehler et al. (2021) used the idea of separating the firs $k$ eigendirections of the covariance to study excess risk of minimum norm interpolators with arbitrary norms and Gaussian data. Bartlett et al. (2021) obtained results which belong to the intersection of the first and the second categories which we described in Section 1.2 (see their Theorem 4.1). Shamir (2022) constructed an example of a misspecified setting (i.e., the noise is not independent from the data) in which our results don’t hold even though the condition number of the matrix $A_{k}$ is a constant (see their Example 1).

Ridge regression setup

where $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{p}$ is the non-increasing sequence of eigenvalues of $\Sigma$.

We assume sub-Gaussianity: denote $Z:=X\Sigma^{-1/2}$ (whitened data matrix). Rows of $Z$ are isotropic centered i.i.d. random vectors. We assume that rows of $Z$ are sub-Gaussian with sub-Gaussian norm $\sigma_{x}$ as defined in Appendix A.1.

Sub-Gaussianity is a classical assumption, which provides a convenient framework for controlling deviations of various quantities of interest (see Vershynin (2018) for an introduction). We discuss whether it is actually needed in Section 6.4.

2 Response model

where $\varepsilon$ is the noise vector. We assume that components of $\varepsilon$ are i.i.d. centered random variables with variance $v_{\varepsilon}^{2}$.

3 Learning procedure

Ridge regression with regularization parameter $\lambda$ is a classical learning algorithm that estimates $\theta^{*}$ from $X,y$ according to the following formula:

See Appendix B for a discussion. The matrix $\lambda I_{n}+XX^{\top}$ will play an important role in our analysis, so we denote

In the ridgeless case ($\lambda=0$), $A$ is the Gram matrix of the data. Ridge regularization shifts all its eigenvalues by $\lambda$.

4 Excess risk and its bias-variance decomposition

The quantity of interest is excess risk that we define in the following way: recall that $x$ is a new data point from the same distribution as rows of $X$. The error that our predictor incurs on this data point is $x^{\top}(\hat{\theta}(y)-\theta^{*})$. We define excess risk as the average squared error over the population, i.e.,

where we define $\|x\|_{M}:=\sqrt{x^{\top}Mx}$ for any positive semi-definite (PSD) matrix $M$ and any vector $x$ of the corresponding dimension.

Note that $\hat{\theta}(y)$ is linear in $y$, which allows us to write

The term $\|\hat{\theta}(X\theta^{*})-\theta^{*}\|_{\Sigma}^{2}$ is the error in the noiseless regime; it is caused by rows of $X$ not spanning the whole space and by regularization. The term $\|\hat{\theta}(\varepsilon)\|_{\Sigma}^{2}$ is the error of learning the zero function from pure noise. One can see that these two terms nicely decouple from each other and can be studied separately. Moreover, note that $\|\hat{\theta}(\varepsilon)\|_{\Sigma}^{2}$ is a quadratic form in $\varepsilon$. Its expectation scales linearly with $v_{\varepsilon}^{2}$ (variance of the noise):

If the noise is sub-Gaussian with sub-Gaussian norm $\sigma_{\varepsilon}$, then by Lemma 22 from the appendix for some absolute constant $c$ and any $t>1$, with probability at least $1-ce^{-n/c}$,

These quantities don’t depend on the distribution of the noise. The goal of this paper is to provide sharp non-asymptotic bounds for them.

The story of separating the first k𝑘k eigendirections and our contribution

Before we present our results, we develop some intuition by considering two easy scenarios: "essentially low-dimensional" and "essentially high-dimensional". For each scenario we will do an informal computation of the excess risk and give a geometric interpretation.

where $\Pi_{X}$ is the projection on the span of columns of $X$. This allows us to write the following informal computation, which leads to the classical $k/n$ rate:

Here we used the informal transition $\|n^{-1}X^{\top}X-\Sigma\|\approx 0$ — the population covariance matrix is well-approximated by the sample covariance matrix uniformly in all directions. If $k\ll n$ this results holds with very few additional assumptions (see (Tikhomirov, 2017) and references therein).

Such a result leads to a classical bias-variance trade-off: the larger the model is, the better it can approximate the true dependence, but also the more noise it picks up. A classical cartoon is shown in Figure 1: Figures 1(a)–1(c) show the result of performing least squares regression with features $\{\cos(m\pi x)\}_{m=0}^{p}$. As the number of features grows, the ability of the model to approximate the signal grows too, but at the cost of increasing sensitivity to the noise. As the number of features approaches the number of data points (the "interpolation threshold"), this leads to overfitting.

According to our definitions of bias and variance from Equation (2) with $\lambda=0$,

Here we see the following: the matrix $X^{\top}(XX^{\top})^{-1}X$ is the projection on the span of the data. This is a random $n$-dimensional subspace in $p$-dimensional space. Thus, with high probability $\|X^{\top}(XX^{\top})^{-1}X\theta^{*}\|^{2}/\|\theta^{*}\|^{2}\approx n/p$, so the projection only preserves an $n/p$ fraction of the energy of the signal. When it comes to the variance term, we can use the same concentration result for the sample covariance as we did in the low-dimensional case, but for the transposed data matrix, meaning $XX^{\top}\approx pI_{n}$. Finishing the computation yields

We see that the signal is almost not learned at all in this regime (the bias term is close to the full energy of the signal), but the noise is also damped by the factor $p/n$.

The geometric interpretation is as follows: if $p\gg n$, the span of $n$ data points is almost orthogonal to $\theta^{*}$ with high probability. The data just does not measure $\theta^{*}$ in most directions, so almost the whole signal is lost. On the other hand, despite the noise fully propagating into in-sample predictions, a new data point $x$ is also almost orthogonal to all the old data points with high probability, so those noisy predictions don’t influence the prediction in $x$. Overall, despite interpolating the data, we effectively learn a zero estimate out of sample. The zero estimator can be a very good estimator, e.g., if the true signal is zero. This hints at the possibility of learning via high-dimensional interpolation: the model can use the directions in which the signal is not learned to smear the noise over them.

The learning cartoon for this regime is given in Figures 1(d)–1(e): as the number of cosine features becomes large compared to the number of data points, the learning procedure predicts zero out of sample, despite interpolating the values in sample. However, if we add multiplicative weights to the cosine features, down-weighting higher frequencies, it causes the minimum norm solution to learn the low frequency signal and interpolate the noise using the high frequency components.

2 The ground provided by the previous work

Bartlett et al. (2020) studied the variance term for ridgeless regression under the additional assumption that the data vectors have independent components. To give an overview of their results, introduce the following quantities: for any $k\in\{0,1,2,\dots,p-1\}$ define

In the ridgeless setting, meaning $\lambda=0$, $r_{0}$ is a well-known effective rank of the matrix $\Sigma$, $r_{k}$ is the effective rank of the same matrix, but after restricting it to the span of its last $p-k$ eigenvectors, and $\rho_{k}$ measures how large that effective rank is compared to the number of data points.

Given this notation, Bartlett et al. (2020) defined $k^{*}$ as the minimum $k$ for which $\rho_{k}$ is larger than a universal constant. Their result is then that if such a $k^{*}$ doesn’t exist or if $k^{*}/n$ is at least a constant, then $V$ is lower bounded by a constant. Otherwise, they show that with high probability $V$ is equal up to a constant factor to the following quantity:

Inspection of the proof shows that the "essentially low-dimensional" rate $k/n$ comes from the first $k$ components of the vector $\hat{\theta}(\varepsilon)$ and the term ${\left(n\sum_{i>k^{*}}\lambda_{i}^{2}\right)}/{\left(\sum_{i>k^{*}}\lambda_{i}\right)^{2}}$ comes from the rest of the components of $\hat{\theta}(\varepsilon)$. Note that if one plugs in $\lambda_{i}=\lambda_{j}$ for all $i,j>k^{*}$, then it becomes ${\left(n\sum_{i>k^{*}}\lambda_{i}^{2}\right)}/{\left(\sum_{i>k^{*}}\lambda_{i}\right)^{2}}=n/(p-k^{*})$ — exactly the variance term of the "essentially high-dimensional" regime of Section 3.1. The conclusion of Bartlett et al. (2020) is therefore that the only way that an interpolating solution can damp the noise by more than a constant factor is the following: the data is such that after removing $k$ components, it becomes "essentially high-dimensional", meaning that the effective rank of its covariance is large compared to the number of data points. After that the variance in the first $k$ components is the same as for the classical least squares, and the variance in the rest of the components corresponds to the "essentially high-dimensional" case, where you cannot learn but the noise is still damped. Note, however, that that story was not complete because only the variance term was bounded sharply in that work.

3 Our contribution

We complete the story of Bartlett et al. (2020) by providing sharp bounds on the bias term, extending the results to the setting of ridge regression with nonzero $\lambda$, and replacing the assumption of independence of the components by a much broader sufficient condition. From our point of view, $k^{*}$ is the main discovery of Bartlett et al. (2020). In our work we also start with separation of the first $k$ eigendirections and show that the same split leads to a bound for the bias term that is in full alignment with the intuitive explanation given above.

The central object in our analysis is the following matrix:

The matrix $X_{k:\infty}X_{k:\infty}^{\top}$ is the Gram matrix of the data after removing the first $k$ components. $A_{k}$ is obtained from that Gram matrix by shifting all eigenvalues by the ridge regularization parameter $\lambda$.

In Bartlett et al. (2020), the crucial step was to show that the singular values of $A_{k}$ are within a constant factor of each other for $k=k^{*}$ (see their Lemma 5). When the components of data vectors are independent, such control over the condition number is a consequence of high effective rank. In this paper, the roles of effective rank and condition number of $A_{k}$ are reversed. We prove sharp bounds assuming that there is some oracle that guarantees that with high probability all eigenvalues of $A_{k}$ are within a constant factor of each other. Independence of components is not needed. Moreover, such control implies that $\rho_{k}$ is at least a constant, which, in turn, implies sharpness of the bounds. In other words, we provide a more general condition under which the tail of the data is "essentially high dimensional" — instead of assuming independent components and high effective rank, only oracle control of condition number of $A_{k}$ is needed. In Section 5 we provide an extensive discussion of this assumption: we show that a version of a small-ball condition for the tails of the data is required and that a stronger version of the same condition is sufficient if the data is sub-Gaussian.

The bound that we obtain for the bias term is given informally by the following expression:

One can see how it aligns with the intuition of “essentially low-dimensional” and “essentially high-dimensional” parts: one cannot estimate the signal in the high dimensional part, so almost all of its energy $\|\theta^{*}_{k:\infty}\|_{\Sigma_{k:\infty}}^{2}$ goes into the error. When it comes to the low-dimensional part, the high-dimensional part acts as a ridge regularizer for it, so the bias in the first $k$ components is the same as that of ridge regression with regularization coefficient $\lambda+\sum_{i>k}\lambda_{i}$ (i.e., the full regularization is equal to the explicitly imposed part $\lambda$ plus "implicit regularization", which is equal to the energy of the tail.)

Our extension of the results to the ridge regression scenario allows us to answer the following question: can it happen that the "essentially high dimensional part" is too high dimensional, meaning that it provides too much regularization and negative $\lambda$ is needed to compensate for that? In Section 8, we show that this indeed can happen and that the following is sufficient for it to be true: the noise and the energy of the signal in the tail (components ${k:\infty}$) are small compared to the signal in the spikedHere we use the word ”spiked” as in the ”spiked covariance models”, which usually assume that the eigenvalues of $\Sigma_{k:\infty}$ are all equal and of smaller order than eigenvalues of $\Sigma_{0:k}$. One way to interpret our results is that only spiked-covariance-like models can exhibit benign overfitting, and we derive general conditions for a model to be spiked-covariance-like. part (components ${0:k}$), but the effective rank of the tail abruptly becomes much larger than $n$.

4 Additional notation

Throughout the paper the following objects will be needed: for any $i$ denote $z_{i}$ to be the $i$-th column of $Z$. Then define

an analogue of the matrix $A$, but we throw away the $i$-th component of the data vectors. Denote also

the ratio of the effective rank of the tail to the number of data points without taking regularization $\lambda$ into account.

For the readers convenience, we compile all the notation in Appendix A.

Main results

As we explained in the previous section, the central objects in our proof are $A_{k}$ and $\rho_{k}$. In principle, any control of the spectrum of $A_{k}$ leads to some upper bound on $B$ and $V$ (see our Theorem 5), the question is when that bound is tight. The intuitive answer is the following: the bound is tight when the condition number of $A_{k}$ is a constant and $k$ is chosen correctly, meaning that either $\rho_{k}$ is a constant or $k$ is the smallest number such that $\rho_{k}$ is larger than a constant (i.e., $k=k^{*}$).Note that there may be several values of $k$ that satisfy these conditions. Applying our upper bound for any of those $k$ will yield the same result up to a constant factor. Our arguments, however, only support this intuition when the following technical assumption holds for some constant $\gamma<1$:

Assume that $\lambda>-\gamma\sum_{i>k}\lambda_{i}$.

The focus of our work was to obtain the tight upper bound on the excess risk under minimal assumptions. Such minimal assumption turns out to be

Assume that with probability at least $1-\delta$ the matrix $A_{k}$ is positive-definite (PD) with condition number at most $L$.

We provide a thorough discussion of this assumption in Section 5, for example we derive sufficient and almost matching necessary conditions for it to hold when the distribution is sub-Gaussian. The reason why we don’t just assume those sufficient conditions is that we believe that sub-Gaussianity is not essential for our results to hold, as we discuss in Section 6.4. Moreover, the matrix $A_{k}$ is the central object in our argument, and making an assumption on its condition number explicitly makes presentation easier.

A careful reader will notice that we have just stated that another condition is needed for the bound to be tight: $k$ should be chosen in the right way. This, however, can be achieved by shifting $k$ to $k^{*}$ if necessary: indeed, assumptions NoncritReg($k,\gamma$)$(k,\gamma)$ and CondNum($k,\delta,L$)$(k,\delta,L)$ imply a constant lower bound on $\rho_{k}$ (see Corollary 6). That means that either $\rho_{k}$ is a constant, or it is more than a constant, i.e., $k>k^{*}$. In the latter case one can shift from $k$ to $k^{*}$ meaning that Assumption CondNum($k,\delta,L$)$(k^{*},\delta^{\prime},L^{\prime})$ also holds with modified constants $\delta^{\prime},L^{\prime}$ (see Lemma 11 for the exact statement). Now applying the upper bound (Corollary 6) with $k=k^{*}$ gives tight result, as given by the following

Fix any constants $b>0,$ $\gamma\in[0,1),$ $L>0.$ Denote

There exists a constant $c$ which only depends on $\sigma_{x}$, $b$, $\gamma$, $L$ such that the following holds: suppose NoncritReg($k,\gamma$)$(\bar{k},\gamma)$ and CondNum($k,\delta,L$)$(\bar{k},\delta,L)$ are satisfied for some $\bar{k}<n/c$ and $\delta<1-ce^{-n/c}$. Take $k=\min(\bar{k},k^{*})$. Then with probability at least $1-ce^{-n/c}-\delta$

Moreover $\rho_{k}\geq c^{-1}$, NoncritReg($k,\gamma$)$(k,\gamma)$ holds, and there exist $L^{\prime},c^{\prime}$ that only depend on $\sigma_{x},b,\gamma,L$ s.t. CondNum($k,\delta,L$)$(k,\delta+c^{\prime}e^{-n/c^{\prime}},L^{\prime})$ holds.That is, the assumptions still hold if we substitute $\bar{k}$ by $k$, but with different $L,\delta$. Further we will see that satisfaction of these assumptions implies tightness of the bounds for the chosen $k$.

Proof In this proof let’s call any quantities that only depend on $\sigma_{x}$, $\gamma$, $b$ and $L$ "constants". First of all, if $\bar{k}\leq k^{*}$ then $k=\bar{k}$. Since we are given that NoncritReg($k,\gamma$)$(\bar{k},\gamma)$ and CondNum($k,\delta,L$)$(\bar{k},\delta,L)$ are satisfied, we immediately get that NoncritReg($k,\gamma$)$(k,\gamma)$ and CondNum($k,\delta,L$)$(k,\delta+c^{\prime}e^{-n/c^{\prime}},L^{\prime})$ are satisfied with $L^{\prime}=L$ and any $c^{\prime}>0$. However, if $\bar{k}>k^{*}$ then $k=k^{*}$ and by Lemma 11 NoncritReg($k,\gamma$)$(k,\gamma)$ and CondNum($k,\delta,L$)$(k,\delta+c^{\prime}e^{-n/c^{\prime}},L^{\prime})$ are still satisfied for some constants $c^{\prime},L^{\prime}$. Note that the larger the constants, the looser the assumptions, so we can take our final choice of $c^{\prime},L^{\prime}$ to be the maximum over two cases.

Now that we know that NoncritReg($k,\gamma$)$(k,\gamma)$ and CondNum($k,\delta,L$)$(k,\delta+c^{\prime}e^{-n/c^{\prime}},L^{\prime})$ are satisfied, by Corollary 6, there is a constant $c_{1}$ such that $\rho_{k}>1/c_{1}$ and with probability at least $1-c_{1}e^{-n/c_{1}}-c^{\prime}e^{-n/c^{\prime}}-\delta$

Taking $c\geq c_{1}+c^{\prime}$ gives the first part.

Algebraically, under Assumption CondNum($k,\delta,L$)$(k,\delta,L)$ all eigenvalues of $A_{k}^{-1}$ are within a constant factor of each other, so one can pull its operator norm from the expressions and obtain an upper bound without losing tightness. This strategy, however, doesn’t produce lower bounds, so we derive them in a different way: we decompose bias and variance into sums with respect to individual coordinates of the predictor, and bound each term in each sum from below. Because of that, we impose different assumptions, namely

Assume that all elements of matrix $X$ are independent (i.e., data vectors have independent coordinates).

Assume that the sequence of coordinates of $\Sigma^{-1/2}x$ is exchangeable (any deterministic permutation of the coordinates of whitened data vectors doesn’t change their distribution).

Assume that $\theta^{*}$ is sampled from a prior distribution in the following way: one starts with vector $\bar{\theta}$ and flips signs of all its coordinates with probability $0.5$ independently.

for the bias term. Because of this mismatch in assumptions, our lower bounds don’t show that our upper bound is always tight. What they show is that one needs some specific knowledge about the distribution to obtain better bounds. We provide a more detailed discussion of the relations between those assumptions in Section 6.2. The lower bounds themselves are given by the following

Fix any constants $b>a>0$, $\gamma\in[0,1),$ $L>0.$ Denote

There exists a constant $c$ which only depends on $\sigma_{x}$, $a$, $b$, $\gamma$, $L$ such that all the following hold:

For any $k\in\{0,1,\dots,k^{*}\}$ under assumptions IndepCoord, NoncritReg($k,\gamma$)$(k,\gamma)$, if $\rho_{k}>a$ then with probability at least $1-2\delta-ce^{-c/n}$

For any $k\in\{1,2,\dots,k^{*}\}$ under assumptions NoncritReg($k,\gamma$)$(k,\gamma)$, CondNum($k,\delta,L$)$(k,\delta,L)$, PriorSigns($\bar{\theta}$)$(\bar{\theta})$ and ExchCoord, if $\rho_{k}>a$ then with probability at least $1-2\delta-ce^{-c/n}$

Proof Lemma 7 gives a lower bound for $V$, and Lemmas 8 and 9 give the lower bound for B. Those lower bounds have the desired probability, but different algebraic form. To bring them to the same form as the upper bounds one needs the right $k$ to be chosen. We assumed that $\rho_{k}>a$. Moreover, since $k\leq k^{*}$ by definition of $k^{*}$ we either have $\rho_{k}\leq b$ or $k=k^{*}$. In both of those cases Theorem 10 guarantees that these lower bounds are the same as what we need up to multiplicative constants that only depend on $\sigma_{x}$, $\gamma$, $a$, $b$ and $L$.

One can notice from this proof that having separate arguments for the lower bounds results in a different algebraic form of the same bound. This different form turns out to be convenient to draw explicit connections between our results and results from earlier works. We do that in Section 7.

The central assumption that we need to compute the excess risk is Assumption CondNum($k,\delta,L$)$(k,\delta,L)$, which provides control over condition number of $A_{k}$. In this section we discuss when this assumption is known to be satisfied and what are the necessary conditions for it to happen.

Recall that $A_{k}=X_{k:\infty}X_{k:\infty}^{\top}+\lambda I_{n}$, so its spectrum is the shift by $\lambda$ of the spectrum of $X_{k:\infty}X_{k:\infty}^{\top}$, the random matrix that is equal to the Gram matrix of the projected data. There are therefore three ways of establishing a constant upper bound on the condition number of $A_{k}$:

Establish an upper bound $\bar{\mu}$ on $\mu_{1}(X_{k:\infty}X^{\top}_{k:\infty})$ and take $\lambda>\bar{\mu}/c$ for some constant $c>0$. In this case, the singular values of $A_{k}$ are all equal to $\lambda$ (and greater than $\bar{\mu}$) up to a constant multiplier.

Establish upper and lower bounds $\bar{\mu}$ and $\underline{\mu}$ on $\mu_{1}(X_{k:\infty}X^{\top}_{k:\infty})$ and $\mu_{n}(X_{k:\infty}X^{\top}_{k:\infty})$ respectively, such that $\bar{\mu}/\underline{\mu}$ is a constant. Then take $\lambda>-\underline{\mu}/c$ for some constant $c>1$. In this case, the singular values of $A_{k}$ are all equal to $\bar{\mu}$ (or $\underline{\mu}$) up to a constant multiplier.

Establish upper and lower bounds $\bar{\mu}$ and $\underline{\mu}$ on $\mu_{1}(X_{k:\infty}X^{\top}_{k:\infty})$ and $\mu_{n}(X_{k:\infty}X^{\top}_{k:\infty})$ respectively, and take $\lambda=-\underline{\mu}+\Diamond$, where $\Diamond\geq c(\bar{\mu}-\underline{\mu})$ for a constant $c>0$. In this case, the singular values of $A_{k}$ are all equal to $\Diamond$ up to a constant multiplier. This case can be substantially different from the previous case when the singular values of $X_{k:\infty}X^{\top}_{k:\infty}$ are very well concentrated, i.e., the gap $\bar{\mu}-\underline{\mu}$ is of smaller order than $\underline{\mu}$ itself. In this case $\Diamond$ can be a smaller order term.

Our bounds are sharp when assumption NoncritReg($k,\gamma$)($\gamma$) is satisfied for some $\gamma<1$, i.e., in the first and the second case above. The third case is quite rare because it requires very good concentration of the spectrum of $X_{k:\infty}X^{\top}_{k:\infty}$. Moreover, in this case $\lambda$ is very close to the critical negative value under which it is impossible to even guarantee that $A_{k}$ is PD as it becomes negative definite in expectation. We use this regime to investigate how negative regularization can improve excess risk by more than a constant factor in Section 8. However, we don’t expect our bounds to always be sharp in this regime.

Therefore, we focus our attention on the first two cases. In Section 5.2 we discuss informally what conditions on the distribution are necessary to bound $\mu_{1}(X_{k:\infty}X^{\top}_{k:\infty})$ and $\mu_{n}(X_{k:\infty}X^{\top}_{k:\infty})$, and show how notions of high effective rank and norm concentration condition arise. In Section 5.3 we combine those bounds for sub-Gaussian data with the choice of $\lambda$ to provide necessary and almost matching sufficient conditions for the condition number of $A_{k}$ to be constant under sub-Gaussianity. In Section 5.4 we show that sub-Gaussianity is not actually required for the condition number of $A_{k}$ to be controlled with high probability: Theorem 4 states that norm concentration condition and a modified version of high effective rank condition are sufficient even if the data only has bounded $4+\varepsilon$ moments.

2 Informal necessary conditions

There are several easy observations that help understand what is needed for the condition number of $A_{k}$ to be bounded.

The first observation is that $X_{k:\infty}X^{\top}_{k:\infty}\succeq\lambda_{k+1}z_{k+1}z_{k+1}^{\top}$, where $z_{k+1}$ is the first column of $Z_{k:\infty}$ —a vector with $n$ i.i.d. coordinates with unit variance. By the law of large numbers, $\|z_{k+1}\|^{2}\approx n$, meaning that $\|\lambda_{k+1}z_{k+1}z_{k+1}^{\top}\|\approx\lambda_{k+1}n$. Therefore, $\bar{\mu}\gtrsim\lambda_{k+1}n$.

The second observation is that the diagonal elements of $X_{k:\infty}X^{\top}_{k:\infty}$ are squared norms of the tails of data vectors. Recall that we denoted the data points to be $\{x^{i}\}_{i=1}^{n}$. We can write

The third observation is that the diagonal elements of a PD matrix themselves provide bounds on the singular values:

Therefore, to control condition number of $X_{k:\infty}X^{\top}_{k:\infty}$ by a constant $L$ with probability $1-\delta$, it is necessary to guarantee that

i.e., $n$ independent random draws of the random variable $\|x_{k:\infty}\|^{2}$ should all lie within a constant factor of some value, meaning that the norm of the tail of a data vector should be within a constant factor of a fixed value with probability $(1-\delta)^{1/n}$.

3 Controlling condition number under sub-Gaussianity

Sub-Gaussianity of the data implies an upper bound on $\mu_{1}(A_{k})$, but doesn’t help with $\mu_{n}(A_{k})$. To see this one can consider a well-known construction: take a sub-Gaussian distribution and construct another distribution in the following way: to sample from this new distribution take a vector from the old distribution and multiply it by $\sqrt{2}$ with probability $1/2$ and by zero otherwise. The new distribution is still sub-Gaussian with the same covariance, but the Gram matrix of $n$ i.i.d. samples from it is degenerate with probability at least $1-2^{-n}$. Therefore, an additional assumption is needed to lower bound $\mu_{n}(A_{k})$. As we already mentioned in Section 5.2, we need norm concentration. Since sub-Gaussianity allows to bound the norm from above, it reduces to a version of the small-ball condition: $\|x_{k:\infty}\|$ should be lower-bounded with high probability. The formal result is given by the following

For any $\gamma\in[0,1)$ and $\sigma_{x}>0$ there exists $c>0$ that only depends on $\sigma_{x}$ and $\gamma$ such that under Assumption NoncritReg($k,\gamma$)$(k,\gamma)$ the following holds: for any $L\geq 1$

If $\rho_{k}\geq L^{2}$ and with probability at least $(1-\delta)^{1/n}$

then with probability at least $1-\delta-ce^{-n/c}$

Suppose that it is known that with probability at least $ce^{-n/c}$ $\mu_{n}(A_{k})\geq L^{-1}\mu_{1}(A_{k})$. Then $\rho_{k}\geq\frac{1}{cL}$ and with probability at least $\left(1-ce^{-n/c}\right)^{1/n}$

The proof is given in Appendix D. One can see that both the necessary and the sufficient conditions are that $\rho_{k}$ is lower bounded by a constant and a version of small-ball condition that says that the regularized squared norm of the data exceeds a constant fraction of its expectation with probability $(1-\delta)^{1/n}$. There is, however, a gap in those constants.

4 Heavy-tailed case

The following is a direct corollary of Theorem 2.1 from Guédon et al. (2017)

Suppose that the distribution of the tail satisfies the following two assumptions:

Norm concentration: For some $\delta\in(0,1/n)$, $L>1$ and $M>0$

Heavy-tailed effective rank: for some $h>4$ denote $r_{h,k}>0$ to be the maximum number such that for any $a\in\mathcal{S}^{p-k-1}$ and $t>0$

There exists a constant $c$ that only depends on $h$ such that with probability at least $1-cn^{1-h/4}-n\delta$

Proof First, note that by union bound with probability at least $1-n\delta$ all the diagonal elements of the matrix $X_{k:\infty}X_{k:\infty}^{\top}$ belong to the segment $[L^{-2}M^{2},L^{2}M^{2}]$. Next, take the bound on $B_{k}$ from the Case 1 of Theorem 2.1 from Guédon et al. (2017) with the following choice of their parameters: $k=N$, $\tau=1$, $\lambda=p$, $\sigma=1+p/4$, $t=\sqrt{n}$. Use that bound for vectors $\sqrt{r_{h,k}}x^{i}_{k:\infty}/M$. Note that that $B_{k}$ is exactly the operator norm of the off-diagonal part of $r_{h,k}X_{k:\infty}X_{k:\infty}^{\top}/M^{2}$.

The quantity $r_{h,k}$ that we introduced in Theorem 4 can be interpreted as a notion of effective rank for heavy tailed distributions. Indeed, one can write

— the ratio of the typical norm of the random vector to the scale of the worst case deviations of its one-dimensional projection. This is completely analogous to our usual definition of the effective rank: $r_{k}=\lambda_{k+1}^{-1}\sum_{i>k}\lambda_{i}$. Indeed, in sub-Gaussian case $\sqrt{\sum_{i>k}\lambda_{i}}$ is the typical value of the norm of the vector $x_{k:\infty}$, and $\sqrt{\lambda_{k+1}}$ is up to constant the largest sub-Gaussian norm of its one-dimensional projection. We see that the conditions under which the eigenvalues of $X_{k:\infty}X_{k:\infty}^{\top}$ are within a constant factor of each other with high probability remain the same even in the heavy-tailed case: the norm of $\|x_{k:\infty}\|$ concentrates within a constant factor of a fixed quantity, and the heavy-tailed effective rank $r_{h,k}$ should be large compared to the number $n$ of data points.

Structure of the proof and role of sub-Gaussianity

The core of our argument is Theorem 5 given below. There are two important things to note about it: first, it only requires sub-Gaussianity and matrix $A_{k}$ being positive semidefinite (which always holds with probability $1$ for non-negative $\lambda$). Second, its proof decomposes very clearly into two parts: an algebraic part, which only requires $A_{k}$ being PD and holds with probability $1$ conditionally on this event, and a probabilistic part, where standard concentration results are directly plugged into the algebraic bounds. Because of this decomposition, it is straightforward to track how the sub-Gaussianity is used and how it can be relaxed. We provide the sketch of the proof to show these details.

There exists a (large) constant $c$, which only depends on $\sigma_{x}$, s.t. for any $k<n/c$ with probability at least $1-ce^{-n/c}$, if the matrix $A_{k}$ is PD, then

Proof sketch The full proof of Theorem 5 can be found in Section I.1 of the appendix. The following is a sketch of its derivation.

Recall the following notation: for any $y$

As explained in Section 3.2, Bartlett et al. (2020) introduced the notion of $k^{*}$ for which the behaviour of the variance term in the first $k^{*}$ coordinates is qualitatively different than in the rest of the coordinates. Their argument, however, relies crucially on independence of the components of the data. The main idea that allowed us to get rid of that assumption and to obtain the tight bound for the bias term was to separate the first $k$ coordinates from the very beginning and to use some sort of uniform convergence argument in that low-dimensional subspace.

The crucial tool that allowed us to realise this idea turned out to be the following algebraic identity that we prove in Section F of the appendix:

This identity allows convenient access to the error in the first $k$ coordinates (the spiked part).

The argument decomposes clearly into two parts: algebraic and probabilistic. The algebraic part is to decompose the excess risk (up to a constant multiplier) into four terms and show that the following inequalities hold on the event that the matrix $A_{k}$ is PD: (1) Bias error in the spiked part:

The probabilistic part of the argument is to control the quantities that arise in the algebraic bound with high probability. Namely, we plug in

Concentration of $k$-dimensional sample covariance with $n$ samples: w.h.p.

Concentration of norm of vectors with i.i.d. components: w.h.p.

After plugging in the probabilistic bounds, the final result is obtained by a straightforward computation.

Note that the only probabilistic statements that are used in this proof are concentration of sample covariance in dimension $k$ and concentration of the sum of $n$ i.i.d. random variables. The same concentration results hold with weaker assumptions, but with larger probability. For example, under rather weak moment assumptions only a linear in dimension number of samples is needed for the sample covariance matrix to concentrate within a constant factor of the population covariance, see Tikhomirov (2017) and references therein. It is also interesting to point out that the "uniform convergence" result that we mentioned in the beginning of this sketch is nothing but the convergence of the empirical covariance matrix $\Sigma_{0:k}^{-1/2}X_{0:k}^{\top}X_{0:k}\Sigma_{0:k}^{-1/2}/n$ to its expectation $I_{k}$, which is exactly the uniform convergence result that gives the bound in the "essentially low-dimensional" regime from Section 3.1.

Despite the fact that the bounds of Theorem 5 apply under very general assumptions, we don’t expect them to be tight if the condition number of $A_{k}$ is not bounded by a constant. When some oracle control of the condition number of $A_{k}$ is provided, the bound becomes the following.

Fix any constants $\gamma\in[0,1)$ and $L>0$. There exists a constant $c$ that only depends on $\sigma_{x}$, $\gamma$, $L$ s.t. for any $k<n/c$ and $\delta<1-ce^{-n/c}$ under assumptions NoncritReg($k,\gamma$)$(k,\gamma)$ and CondNum($k,\delta,L$)$(k,\delta,L)$, it holds that $\rho_{k}>c^{-1}$, and with probability at least $1-\delta-ce^{-n/c}$,

It is also worth mentioning that the story about "essentially high-dimensional" and "essentially low-dimensional" parts is not just an interpretation of the final result: the whole proof strategy is in accordance with it, as we explicitly separate the two parts and bound errors in them separately.

2 Lower bounds

Our lower bounds have a different form from the upper bounds. We show separately that they match if the condition on effective rank is satisfied. One benefit of this approach is that the lower bounds provide a different form of the same result, which allows for different analysis. We employ it in Section 7.

The lower bound for the variance term is given by the following lemma, whose proof is given in Appendix E.1:

Fix any constant $\gamma\in[0,1)$. There exists a constant $c$ that only depends on $\sigma_{x}$ and $\gamma$ s.t. for any $k<n/c$ under assumptions NoncritReg($k,\gamma$)$(k,\gamma)$ and IndepCoord w.p. at least $1-ce^{-n/c}$

One can see that the assumptions under which the lower bound is proved are different from the assumptions required for the upper bound: we require independent components here. On the one hand, it means that there could be a gap between upper and lower bounds in some particular cases where one can control the condition number of $A_{k}$ without independence of components. On the other hand, it means that even such strong additional assumption as independence of components does not allow the upper bounds to be improved, which suggests that those specific cases for which the bound is not tight are rare and require even stronger additional assumptions.

The most general lower bound for the bias term that we prove requires the following assumption

Assume that for any $j\in\{1,2,\dots,p\}$ with probabilityNote that the condition on probability is separate for every $j$, i.e., we don’t assume that events hold simultaneously for all $j$. at least $1-\delta$

and that $\lambda>-\sum_{i>k}\lambda_{i}$.

Then the bound is given by the following lemma, whose proof is given in Appendix E.2

Fix any constant $L>0$. There exists $c$ that only depends on $\sigma_{x}$ and $L$ s.t. for any $k\in\{1,2,\dots,p\}$ under assumptions PriorSigns($\bar{\theta}$)$(\bar{\theta})$ and StableLowerEig($k,\delta,L$)$(k,\delta,L)$ w.p. at least $1-2\delta-ce^{-n/c}$

Assumption StableLowerEig($k,\delta,L$)$(k,\delta,L)$ is formally not comparable to Assumption CondNum($k,\delta,L$)$(k,\delta,L)$, but informally if $k\geq 1$ then StableLowerEig($k,\delta,L$)$(k,\delta,L)$ is weaker: indeed, the matrix $A_{-i}$ is obtained from the matrix $A$ by subtracting $\lambda_{i}z_{i}^{\top}z_{i}$, while the matrix $A_{k}$ is obtained from $A$ by subtracting $\sum_{i=1}^{k}\lambda_{i}z_{i}^{\top}z_{i}$, i.e., the sum of $k$ "largest" of the terms $\lambda_{i}z_{i}^{\top}z_{i}$. Therefore, the matrix $A_{-i}$ is "larger" than $A_{k}$, and controlling its lowest singular value should be easier. The following lemma, whose proof is given in Appendix E.2, formalizes this argument under Assumption ExchCoord:

For any $\gamma<1$ there exists a constant $c$ that only depends on $\gamma$ and $\sigma_{x}$ such that if assumptions CondNum($k,\delta,L$)$(k,\delta,L)$, NoncritReg($k,\gamma$)$(k,\gamma)$ and ExchCoord are satisfied for some $L\geq 1$ and $k\in\{1,2,\dots,p\}$, then StableLowerEig($k,\delta,L$)$(k,\delta+2e^{-n/c},cL)$ is also satisfied.

When it comes to averaging over the prior given by the assumption PriorSigns($\bar{\theta}$)$(\bar{\theta})$, it just means that it is impossible to obtain a better lower bound without some specific knowledge of how signs of components of $\theta^{*}$ interact with the probability distribution of the data.

3 Connecting upper and lower bounds

One slight inconvenience with our approach of imposing oracle control over the spectrum of $A_{k}$ via Assumption CondNum($k,\delta,L$)$(k,\delta,L)$ is the following: what if the oracle provides control for the wrong value of $k$? There can in principle be many values of $k$ for which such oracle control is possible, with not all of them giving the right point where the behaviour changes from "essentially low-dimensional" to "essentially high-dimensional". As an example, consider the isotropic setting with $p\gg n$: one can exclude any number $k$ of components such that $p-k\gg n$ and still be able to control the condition number.

First of all, in accordance with the result of Bartlett et al. (2020), the following theorem shows that the "right $k$" is the $k$ that is not larger than $k^{*}$.

Fix constants $a>0$ and $b>1/n$. There exists a constant $c>0$ that only depends on $a,b$, s.t. the following holds: if either $\rho_{k}\in(a,b)$ or $k=\min\{\kappa:\rho_{\kappa}>b\}$, then

Proof The proof is a rather straightforward comparison of pairs of sums term by term. It is given in Appendix I.2.

Secondly, if the data is sub-Gaussian, then oracle control for any $k<n$ results in tight bounds, but with worse constants. This happens because of the following lemma.

Fix any constants $\gamma\in[0,1)$, $b>0$, $L>0$. Denote

There exist constants $c,L^{\prime}$ that only depend on $\sigma_{x}$, $\gamma$, $b$, $L$ s.t. the following holds: suppose assumptions NoncritReg($k,\gamma$)$(k,\gamma)$ and CondNum($k,\delta,L$)$(k,\delta,L)$ hold for some $k\in[k^{*},n]$. Then assumptions NoncritReg($k,\gamma$)($k^{*},\gamma$) and CondNum($k,\delta,L$)$(k^{*},\delta+ce^{-n/c},L^{\prime})$ hold too.

Proof sketch Since $k\geq k^{*}$, $\mu_{n}(A_{k})$ provides a lower bound for $\mu_{n}(A_{k^{*}})$. When it comes to $\mu_{1}(A_{k^{*}})$, it can be bounded with high-probability because the data is sub-Gaussian. The full proof is given in Appendix D.

4 The role of sub-Gaussianity

As can be seen from the proof of Theorem 1, the strategy to obtain a tight bound is the following: ask the oracle to control the condition number of $A_{k}$, if that $k$ is too large, shift it to $k^{*}$, and then apply the bound from Corollary 6. In Section 5.4 we showed that if the norm $\|x_{k:\infty}\|$ concentrates, and the effective rank $r_{h,k}$ is high enough, then the control over the condition number of $A_{k}$ is possible even if we have very weak moment assumptions instead of sub-Gaussianity. Moreover, as we have discussed in the proof sketches, if we didn’t shift from $k$ to $k^{*}$, we would only need the usual concentration results such as the law of large numbers or concentration of $k$-dimensional empirical covariance matrix with $n$ samples, which also hold under weak moment assumptions. Therefore, sub-Gaussianity is not essential to obtain the bound in the form given in Corollary 6, one just needs to substitute the sub-Gaussian concentration results with their heavy-tailed analogues. However it may not necessarily give a tight result unless the oracle is guaranteed to choose the appropriate $k$ (e.g., $k=k^{*}$). To shift from $k$ to $k^{*}$ we also need an upper bound on $\|A_{k^{*}}\|$, which we derive from sub-Gaussianity. According to Section 5.4, an analogous bound is still possible under weak moment assumptions, but additional work is required: to use Theorem 4 for $k=k^{*}$ one would need to obtain a high-probability upper bound on $\|x_{k^{*}:\infty}\|$ under moment assumptions and to relate $r_{k}$ which we use in definition of $k^{*}$ to $r_{h,k}$, which is introduced in Theorem 4.

Alternative forms of the bounds and effect of increasing regularization

Theorem 10 reveals an alternative form of the bounds: when $\rho_{k}$ is lower- and upper-bounded by constants or when $k=k^{*}$, the bounds on the bias and variance respectively become equal to the following up to a constant multiplier:

These expressions closely resemble the classical expressions for the in-sample bias and variance of ridge regression. Indeed, a straightforward computation gives

where $\{\hat{\lambda}_{i}\}_{i=1}^{p}$ are eigenvalues of the empirical covariance $n^{-1}X^{\top}X$ and $\{v_{i}\}_{i=1}^{p}$ are the corresponding eigenvectors. Recall that $\rho_{k}\lambda_{k+1}=\left(\lambda+\sum_{i>k}\lambda_{i}\right)/n$. One can see that Equations (6)–(7) can be obtained from the classical equations for the in-sample risk by substituting the empirical eigenvalues with population eigenvalues and increasing the regularization level $\lambda$ by $\sum_{i>k}\lambda_{i}$ — the energy of the tail of the data.

2 Dependence on λ𝜆\lambda

The alternative form of the bounds presented in Section 7.1 provides a convenient way to investigate the dependence on $\lambda$, which is cumbersome in the initial form because increasing $\lambda$ may decrease $k^{*}$. This effect, however, is negligible when Equations (6)–(7) are considered. Indeed, in Appendix I.3 we show the following

Suppose $k<n/c$ for some $c>1$ and $k^{*}<k$. Then

Because of this lemma, any $k\in[k^{*},n/c]$ gives the same result (up to a constant factor) in Equations (6)–(7). One can, therefore, start with some $\lambda$ and the corresponding $k=k^{*}$ and then consider larger values of $\lambda$ without decreasing $k$ in Equations (6)–(7). The result will give sharp (up to a constant factor) bounds, which depend on $\lambda$ as follows:

which are obtained by simply plugging in the definition of $\rho_{k}$ into (6)–(7).

A particularly interesting case arises when $\lambda$ is large enough that it dominates $\sum_{i>k}\lambda_{i}$ and all eigenvalues of $A_{k}$ are equal to $\lambda$ up to a constant multiplier. The corresponding result is given by the following corollary.

There is a large positive constant $c$ that only depends on $\sigma_{x}$ such that if

Proof The full proof is given in Appendix I.3; the following is its outline:

Use Lemma 3 to control the eigenvalues of $A_{\lfloor n/c\rfloor}$.

Use Theorem 1 to obtain the bounds for $k=k^{*}$.

Use Theorem 10 to convert the bounds into the form given in Equations (6)–(7).

Use Lemma 12 to substitute $k^{*}$ back with $\lfloor n/c\rfloor$.

Since $\lambda>2\sum_{i>k}\lambda_{i}$, $\lambda/n$ is equal to $\rho_{k}\lambda_{k+1}$ up to a multiplicative constant.

Note that the statement of Corollary 13 does not require the notion of $k^{*}$.

3 Comparison with other results

As we saw in the previous section, the alternative form given by Equations (6)–(7) has milder dependence on the choice of $k^{*}$ than our main bounds (4)–(5) and allows to compare to classical results for in-sample error of ridge regression. In this section we use it to compare with more recent developments: the non-asymptotic bounds in Hsu et al. (2014) and Hastie et al. (2020).

First of all, we follow Hsu et al. (2014) and introduce the following notion of effective dimension of the problem:

where $\bar{\lambda}$ is a parameter which can informally be understood as effective level of regularization. Hsu et al. (2014) provide non-asymptotic bounds for $B$ and $V$ in the regime when

where $c$ is some constant that depends on the concentration properties of the data. This is the same as the result of Corollary 13, but with different constants. However, our Corollary 13 covers a wider range of $\lambda$ if $n$ is large enough. This follows from the following lemma, which is proven in Appendix I.3:

Suppose that $n\geq c^{2}+c$ for some $c>0$ and take

Indeed, $d(\lambda/n)$ is a decreasing function of $\lambda$, and due to Lemma 14 the range of $\lambda$ for which Corollary 13 is applicable when $d(\lambda/n)=O(n)$, while Equation (8) restricts to the range $d(\lambda/n)=O(n/\log n)$.

After we posted the first preprint of this paper, the following non-asymptotic bound for the interpolating regime (i.e., $\lambda=0$) appeared in (Hastie et al., 2020): informally

Suppose that $k<n/c$ and $\rho_{k}>c$ for some constant $c>1$ . Then

which implies $a\rho_{k}+1\geq\rho_{k}$, so $a\geq 1-1/\rho_{k}>1-1/c$.

The similarity of Equations (11)–(12) with our results should not be taken for granted, and it is actually quite surprising. As we explain in Section 9, the regime considered in Hastie et al. (2020) is significantly different, so it is rather unclear why the results would have the same form.

Negative regularization

The aim of this section is to find a family of regimes in which the optimal level of ridge regularization is negative. Since we are comparing different values of $\lambda$ in this section, the following notation will be useful: recall that for any $k$

the value of $\rho_{k}$ for $\lambda=0$. Intuitively, the components of the tail provide regularization for the first $k$ components, and the larger $\rho_{k}$ is, the more is that regularization. Thus, one could expect that if there is an abrupt jump in the sequence $\{\rho_{k}(0)\}_{k=0}^{p}$, then that additional regularization is too large and negative $\lambda$ may be optimal.

As we investigate further, a jump in $\rho_{k}(0)$ is indeed one of the sufficient conditions for optimality of negative regularization, but not the only one: the strength of the noise and how the signal is distributed among the principal components of the data also play an important role.

Now let’s look at the role of the signal in the tail. It contributes to error in two ways: first — the components in the tail are not getting estimated themselves, second — the signal that comes from those components acts as additional noise for estimation of the first $k$ components. When $\lambda$ is non-negative, the error of the first type dominates the error of the second type, but negative $\lambda$ can amplify the noise and result in error of the second type dominating. Therefore, the signal in the tail also needs to be sufficiently small in order for negative regularization to be optimal.

The final observation is the following: since we only compute the bounds up to a constant multiplier, the bound in Theorem 1 cannot distinguish between negative and zero regularization. To see this, consider the form of the bound given in Section 7: up to a constant factor the bound is a weighted combination in each component with weight $\lambda+\sum_{i>k}\lambda_{i}$, and as $\lambda$ increases there is no need to change $k$. Now it is easy to see that for all $\lambda$ in range from $-\gamma\sum_{i>k}\lambda_{i}$ to zero, that weight is the same up to a constant factor. Thus, negative regularization can only decrease the excess risk by more than a constant factor in the critical regime, i.e., $\lambda=-\sum_{i>k}\lambda_{i}+\Diamond$ where $\Diamond$ is of smaller order than $\sum_{i>k}\lambda_{i}$. To consider such $\lambda$ and have $A_{k}$ PD we need tight concentration of eigenvalues of $X_{k:\infty}X_{k:\infty}^{\top}$ around $\sum_{i>k}\lambda_{i}$. To ensure such tight control we restrict ourselves to the case of independent components, i.e., when Assumption IndepCoord is satisfied. In this case, the eigenvalues of $X_{k:\infty}X_{k:\infty}^{\top}$ can be bounded according to the following statement that was shown as an intermediate step in the proof of Lemma S.9 in (Bartlett et al., 2020).

Under assumption IndepCoord there exists a constant $c$ that only depends on $\sigma_{x}$ s.t. with probability at least $1-ce^{-n/c}$,

The fluctuations $n\lambda_{k+1}+\sqrt{n\sum_{i>k}\lambda_{i}^{2}}$ will be of smaller order than $\sum_{i>k}\lambda_{i}$ if $\rho_{k}(0)$ is larger than a constant, which is shown by the following bounds:

Using this lemma allows us to obtain following two lemmas. See Appendix J for the proofs.

There exist constants $b,c$ that only depend on $\sigma_{x}$ such that the following holds: suppose that assumptions IndepCoord and PriorSigns($\bar{\theta}$)$(\bar{\theta})$ hold. Take $k=\min\{\kappa:\rho_{\kappa}(0)>b\}$ and suppose that $k>0$. Then with probability at least $1-ce^{-n/c}$ for any $\lambda\geq 0$

There exists a constant $c$ that only depends on $\sigma_{x}$ such that the following holds: suppose that assumptions PriorSigns($\bar{\theta}$)$(\bar{\theta})$ and IndepCoord hold and that $\rho_{k}(0)>c$ for some $k<n/c$. Assume also that

Then there exists such $\lambda<0$ that with probability at least $1-ce^{-n/c}$

Lemma 17 provides a lower bound on the expected (over noise and $\theta^{*}$) excess risk which holds w.h.p. uniformly over all non-negative $\lambda$. Lemma 18 provides an upper bound that can be achieved by some negative $\lambda$. Combining these two lemmas gives a sufficient condition for the optimal $\lambda$ to be negative, which is given by the following theorem.

Proof It is easy to see that by taking $c$ large enough, the conditions of Lemmas 18 and 17 are satisfied, and the upper bound in Lemma 18 becomes lower than the lower bound in Lemma 17.

We see that the conditions indeed align with the intuition outlined in the beginning of this section: we need small variance, small signal in the tail, and a sharp jump in effective rank. However, we do not have matching lower bounds in the critical regime when Assumption NoncritReg($k,\gamma$)$(k,\gamma)$ is not satisfied for a constant $\gamma>1$. Thus, we don’t know whether these conditions are also necessary.

Comparison to other works

As we mentioned in Section 1.2, recently there has been a number of papers studying population risk of interpolating solutions of linear regression, and we gave a rough split of those results into three categories there. Here we elaborate on the comparison between the approaches and results.

Results from the first category (Dobriban and Wager, 2015; Hastie et al., 2019; Wu and Xu, 2020; Richards et al., 2020) compute exact asymptotic expressions for the excess risk assuming that $p/n$ goes to some constant as $p,n$ go to infinity, and that the spectral distribution of $\Sigma$ converges to some limiting distribution. From the point of view of our approach, such distributions are indistinguishable from isotropic: indeed, the very existence of limiting spectral measure implies that almost all eigenvalues are within a constant factor of each other. Many of those works even assume explicitly that the spectrum of $\Sigma$ is upper- and lower-bounded by two constants (Richards et al., 2020, page 7), (Wu and Xu, 2020, Assumption 1), (Hastie et al., 2019, Theorem 3). Our results don’t need any asymptotic set up, and apply to $p=\infty$ with some fixed summable sequence $\lambda_{i}$, which has no meaningful notion of limiting distribution, and no separation from zero is needed. For example, our setup covers kernel regression with a fixed kernel and increasing number of data points. On the other hand, when all $\lambda_{i}$ are within a constant factor of each other, our lower bounds become $B\geq\|\theta^{*}\|_{\Sigma}/c$ and $V\geq 1/c$, so the constant part of the whole signal doesn’t get learned and the variance term is at least a constant, i.e., the asymptotic expressions obtained in the works from this category are all just different constants and our approach cannot distinguish them. Therefore, we answer significantly different questions: while the asymptotic work distinguishes between constant error rates, we investigate when the error can be less than a constant. The final difference with our work is rather technical but quite strong: all the works in this category assume that the coordinates of the data become independent if multiplied by the inverse square root of the covariance. This assumption stems from asymptotic random matrix theory techniques, on which these papers are based. To the best of our knowledge, it is not known how to extend these techniques beyond random matrices with independent elements. Our approach, however, does not require the coordinates to be independent.

When it comes to the second category, featurized or kernel regression (Montanari and Zhong, 2020; Ghorbani et al., 2020b; Mei and Montanari, 2019; Ghorbani et al., 2020a; Liang et al., 2020), the difference from our approach is that we do not assume any particular mechanism for data generation or how the features are constructed, but we directly make assumptions about feature vectors. Our results can in principle be applied in this setting if one computes the spectrum of the population covariance for particular features or kernels and the corresponding sub-Gaussian norms. The major difficulty that precludes such a direct comparison is that that computation is not straightforward. The works from this category operate in a more particular setting and circumvent the computation of the spectrum of $\Sigma$. On the other hand, it is not hard to trace strong similarities with our approach on the level of the proof. First of all, all the papers in this category that we are aware of assume that the data comes from a very regular distribution: either $d$-dimensional isotropic data with i.i.d. coordinates (Liang et al., 2020, Assumption 1), or data from the uniform distribution on the sphere (Mei and Montanari, 2019; Ghorbani et al., 2020a, abstracts), (Montanari and Zhong, 2020, Section 3.2), or data from the product of two uniform distributions on spheres (Ghorbani et al., 2020b, Section 2.1). Second, in all those papers the kernel is either spherically symmetric (Ghorbani et al., 2020b, Section 2.2), (Liang et al., 2020, Equation 4) or close to being spherically symmetric due to isotropic initialization of the neural network or isotropic choice of random features (Ghorbani et al., 2020a, Assumption 1), (Mei and Montanari, 2019, Thorem 2), (Montanari and Zhong, 2020, Section 3.2). After that, they consider the regime where $n$ is large compared to $d^{\alpha}$ for some $\alpha$ (Montanari and Zhong, 2020, Assumption 1), (Ghorbani et al., 2020b, Theorem 1), (Mei and Montanari, 2019; Ghorbani et al., 2020a; Liang et al., 2020, abstracts)In (Mei and Montanari, 2019) $\alpha=1$.. Finally, all those papers derive that kernel regression works effectively as ridge regression with polynomial features up to degree $\alpha$ (Montanari and Zhong, 2020; Ghorbani et al., 2020a, abstracts), (Ghorbani et al., 2020b, Theorem 1), (Liang et al., 2020, Proposition 1 and Section 2.3). The only exception is Mei and Montanari (2019), who derive asymptotic expressions for excess risk when the true function is affine (i.e., a polynomial of degree $1$) plus Gaussian misspecification. The connection with our results is that in such a regime (uniform distribution on the sphere, spherically symmetric kernel) polynomials are exactly the eigenfunctions of the kernel operator, which plays the role of the covariance operator, and there are $k\approx d^{\alpha}$ of polynomials of degree at most $\alpha$. Thus, their approach is similar to ours: separate the first $k$ eigendirections (or their approximations) and show that other directions act as regularization.

Conclusions

We studied the excess risk of ridge regression and showed how geometry of the data can influence both which part of the signal is learned and how the noise is damped. For a range of values of the regularization parameter we showed that learning can be seen as the composition of two parts: classical ridge regression in the first $k$ components (the "essentially low-dimensional part") and learning the zero estimator in the rest of the components (the "essentially high-dimensional part"). We introduced a general assumption under which the data is “essentially high-dimensional”, and provided geometric sufficient conditions for its satisfaction. Moreover, we investigated the regime in which the “essentially high-dimensional part” is too high-dimensional, and derived general sufficient conditions for negative regularization to be optimal: small noise, small energy of the "essentially low-dimensional part", but an abrupt jump in the effective rank.

On the technical side, our proof decouples cleanly into an algebraic part, which holds with probability 1 for non-negative regularization,For the case of negative regularization we need to condition on the event that all the necessary symmetric matrices are PD. and the probabilistic part, where we plug in well-known concentration results from high-dimensional probability. This makes it easy to trace how different terms in the bound correspond to the parts of the estimator, and supports the geometric interpretation given above.

We provided a thorough overview of the related papers, and explained how our results are significantly different from them despite some optical similarities. Those similarities, however, are intriguing, and hint at the task of developing a unified treatment of different regimes of overparameterized linear regression as a promising direction of future work.

Acknowledgements

We gratefully acknowledge the support of the NSF through grants DMS-2023505 and DMS-2031883 and of the Simons Foundation through award #814639.

Appendix A Definitions and Notation

A random variable $z$ is sub-Gaussian if it has a finite sub-Gaussian norm

The sub-Gaussian norm of a random vector $Z$ is

A.2 Standard mathematical objects

$M[i,j]$ denotes the element of the matrix $M$ which stands at the intersection of the $i$-th row and $j$-th column.

$I_{m}$ is the $m\times m$ identity matrix.

A.3 Data and the learning procedure

components of $\varepsilon$ are independent and have variance $v_{\varepsilon}$,

$x$ denotes a new random draw from the data distribution, i.e., $x$ is independent from $X,\varepsilon$ and $x$ has the same distribution as $x^{1}$.

$Z=X\Sigma^{-1/2}$, $\{z_{i}\}_{i=1}^{p}$ are columns of $Z$.

the rows of $Z$ are sub-Gaussian with sub-Gaussian norm at most $\sigma_{x}$,

ridge regression outputs $\hat{\theta}(y):=X^{\top}(\lambda I_{n}+XX^{\top})^{-1}y$.

A.4 Splitting the coordinates

For some $k<n$ we spit the coordinates into two groups: the first $k$ components and the rest of the components. Thus we introduce the following notation. Consider integers $a,b$ from to $\infty$ (we always either take $a=0$ and $b=k$ or $a=k$ and $b=\infty$).

$A_{k}=\lambda I_{n}+X_{k:\infty}X_{k:\infty}^{\top}$.

$r_{k}=\frac{1}{\lambda_{k+1}}\left(\lambda+\sum_{i>k}\lambda_{i}\right).$