Dimension free ridge regression

Introduction

Statistical theory studies this and similar estimators in three different regimes:

Despite the wealth of fascinating technical results in this area, this state of affairs leaves open many important questions.

First, it would be important have a unified theoretical framework that does not require the statistician to decide which asymptotics to use. For instance, in order to apply sharp asymptotics in the classical or proportional regimes, it is often assumed that a given pair $(n,d)$ is in fact an element of a sequence $(n,d(n))$ with, respectively, either $d(n)\asymp 1$, or $d(n)\asymp n$.

In practice we are given a single pair, say $(n,d)=(1000,50)$: should we interpret this as $d\asymp 1$, $d\asymp n$, or yet another regime that is not covered by current theory (e.g., $d\asymp n^{2/3}$)?

In fact, the distinction between three types of asymptotics outlined above is rather the consequence of the technical tools used to derive them, rather than a fundamental statistical phenomenon.

Second, the restriction $d=O(n)$ (or $s=O(n)$ in sparse regression) which is implied both by the proportional and by the classical asymptotics is artificial. While this condition might seem necessary for consistency at fist sight (it might seem that at least $d$ observations are required to estimate $d$ parameters), as shown in [BLLT20, TB20] this is in fact not the case. Further, it is not even clear how to check in practice $d=O(n)$ for a given pair $n,d$.

Third, it would be important to remove the assumption of a well conditioned ${\bm{\Sigma}}$, and derive precise asymptotics for general covariances. We would argue that the ill-conditioned case is most important in practice, since high-dimensional data have often low-dimensional structures.

Fourth, the proportional asymptotics is somewhat un-natural from a statistical viewpoint. Most statisticians are used to think of the data distribution is fixed (in particular, $d$ is fixed), while we sample size $n$ increases. In a standard proportional setting, one instead assumes $n,d\to\infty$ together with $n/d\to\delta$: the data distribution changes with the sample size.

Recent progress on several of these issues was achieved in the context of ridge regression. Among others, [HMRT22] derived a characterization for bias and variance in the proportional regime that is non-asymptotic, i.e. holds up to an approximation error that is explicit and vanishes for large $n$, $d$. Using a different approach, [BLLT20, TB20] obtained bounds on bias and variance that hold for arbitrary (possibly infinite) dimension $d$, in terms of of the decay of eigenvalue of ${\bm{\Sigma}}$. These bounds allow to demonstrate ‘benign overfitting,’ i.e. choices of ${\bm{\Sigma}},{\bm{\beta}}$ (i.e. data distributions) such that minimum norm interpolator is consistent.

The results [HMRT22, BLLT20, TB20] have limitations. The characterization of the risk proved in [HMRT22] has sharp leading constants, but only holds for $C^{-1}\leq n/d\leq C$ with $C$ a constant, and holds up to an additive error. However, this error terms can be larger than the actual excess risk when the latter vanishes. The bounds of [BLLT20, TB20], on the other hand, hold up to unspecified multiplicative constants. The proof techniques in these two sets of results are furthermore very different.

In this paper we attempt to provide a unified picture that covers these gaps, by extending the sharp characterization of ridge regression of [HMRT22] beyond the proportional regime. This will allows to recover the benign overfitting results of [BLLT20, TB20] (in several cases) with sharp constants. In doing so, we will extend random matrix theory analysis to cases with $d\gg n$ or $d=\infty$, without restrictions on the condition number of ${\bm{\Sigma}}$. In the case $d=\infty$, the feature vectors ${\bm{x}}_{i}$ are random elements in a separable Hilbert space, whose distribution is fixed (does not change with $n$), and whose covariance ${\bm{\Sigma}}$ is a trace class self-adjoint operator.

The rest of the paper is organized as follows. The next section describes the setting for our analysis, the main assumptions and the resulting asymptotic characterization. It also provides some intuition and connects our results to earlier work. Section 3 contains the formal statement of our general results, while Section 4 specializes our theorem to regimes of interest and develops tools to check its assumptions. Section 5 evaluates our characterization for certain choices of ${\bm{\Sigma}},{\bm{\beta}}$, and compare the predictions with simulations. Finally, proof are presented in Sections 6 and 7, with most technical steps deferred to the appendices.

Setting and characterization

the response vector ${\bm{y}}=(y_{1},\cdots,y_{n})^{\sf T}$ and the noise vector ${\bm{\varepsilon}}=({\varepsilon}_{1},\cdots,{\varepsilon}_{n})^{\sf T}$, we can write in matrix form

For an estimator $\widehat{{\bm{\beta}}}=\widehat{{\bm{\beta}}}({\bm{X}},{\bm{y}})$ we define the excess risk as

where ${\bm{x}}_{\mathsf{new}}$ is an independent copy of ${\bm{x}}_{1},\cdots,{\bm{x}}_{n}$ and $\|{\bm{x}}\|_{\bm{\Sigma}}^{2}:={\bm{x}}^{\sf T}{\bm{\Sigma}}{\bm{x}}$. We will also refer to this as the ‘test error’ or the ‘generalization error’ (although the latter is actually given by the difference between $\mathscr{R}_{\bm{X}}$ and ts empirical version.) Let us emphasize that in this definition, $\mathscr{R}_{\bm{X}}(\widehat{{\bm{\beta}}};{\bm{\beta}})$ is a random quantity because it depends on the data ${\bm{X}}$: however, as we will prove, it concentrates around a non-random value.

The generalization error admits a variance-bias decomposition $\mathscr{R}_{\bm{X}}(\widehat{{\bm{\beta}}};{\bm{\beta}})=\mathscr{V}_{\bm{X}}(\widehat{{\bm{\beta}}};{\bm{\beta}})+\mathscr{B}_{\bm{X}}(\widehat{{\bm{\beta}}};{\bm{\beta}})$, with

For ridge regression, we can write explicit forms of variance and bias:

Assumptions on the covariates distribution.

We impose the following assumptions on the covariates ${\bm{x}}_{i}$ throughout the paper.

We further assume ${\bm{x}}_{i}={\bm{\Sigma}}^{1/2}{\bm{z}}_{i}$ where the following hold.

I. There exist $\mathsf{d}_{{\bm{\Sigma}}}:=\mathsf{d}_{{\bm{\Sigma}}}(n)\geq n$ such that, for all $1\leq k\leq\min\{n,d\}$

II. There exist $\mathsf{C}_{{\bm{x}}}>0$, such that one of the following condition holds:

The technical motivation for assumption II is to establish concentration of quadratic forms of ${\bm{z}}_{i}$, via Hanson-Wright inequality. We notice that the convex concentration property is implied by any of the following. $(i)$ By Talagrand inequality, convex concentration holds for random vectors ${\bm{z}}_{i}$ with independent bounded entries [BLM13, Theorem 7.12]. $(ii)$ By Herbst’s argument, concentration of Lipschitz functions (and hence in particular convex concentration) holds for random vectors ${\bm{z}}_{i}$ that satisfy a log-Sobolev inequality [BGL+14, Proposition 5.4.1]. $(iii)$ Finally, as a special case of the last point, vectors ${\bm{z}}_{i}$ with strongly log-concave probability density function satisfy this condition [BGL+14, Corollary 5.7.2].

The form of Hanson-Wright inequality that we will use is given below.

Effective variance and bias.

An important observation of [HMRT22] is that variance $\mathscr{V}_{\bm{X}}$ and bias $\mathscr{B}_{\bm{X}}$ concentrate around some non-random quantities, that can be interpreted in terms of an ‘effective’ regression problem. While [HMRT22] proves such characterization in the proportional regime $n\asymp d$, here we will extend its validity and prove stronger guarantees.

Define the effective regularization $\lambda_{\star}$ as the unique non-negative solution of

we that define the effective variance and bias as

Our main result —stated in the next section— will establish dimension-free guarantees of the form

where the term These improve over earlier work in two important directions. First, they are dimension free, and in particular do not assume $n\asymp d$. Second, they provide multiplicative approximations, and hence retain their utility when the risk is small.

Bounds, interpretation, benign overfitting.

Before stating our formal results relating $\mathscr{V}_{\bm{X}}$ to $\mathsf{V}_{n}$ and $\mathscr{B}_{\bm{X}}$ to $\mathsf{B}_{n}$, it is useful to develop some intuition about the expressions (6), (7) and their immediate consequences. Note that, by Eq. (5), we necessarily have

If we assume that inequality between the first and last term holds with a constant multiplicative factor, i.e. ${\rm{Tr}}\left({{\bm{\Sigma}}^{2}({\bm{\Sigma}}+\lambda_{\star}{\bm{I}})^{-2}}\right)\leq n(1-c_{\star}^{-1})$ for some constant $c_{\star}\in(0,\infty)$, then we get

Comparing these bounds with the bias and variance of general ridge regression in Eqs. (4a), (4a), we observe that the right hand sides are (modulo the factor $c_{\star}$) the bias and variance of a modified ridge regression in which:

The design matrix is non-random and given by ${\bm{\Sigma}}^{1/2}$ instead of ${\bm{X}}$.

The regularization parameter is $\lambda_{\star}$ instead of $\lambda$.

The noise level is $\tau/\sqrt{n}$ instead of $\tau$.

Even more explicit expressions can be obtained by writing the right-hand side of Eqs. (11), (12) in the basis that diagonalizes ${\bm{\Sigma}}$ as in the next proposition. A proof of this statement is in Appendix A.

Assume ${\rm{Tr}}\left({{\bm{\Sigma}}^{2}({\bm{\Sigma}}+\lambda_{\star}{\bm{I}})^{-2}}\right)\leq n(1-c_{\star}^{-1})$, for $c_{\star}\in(1,\infty)$. Let ${\bm{\Sigma}}:=\sum_{i\geq 1}\sigma_{i}{\bm{v}}_{i}{\bm{v}}_{i}^{{\sf T}}$ be the eigendecomposition of of ${\bm{\Sigma}}$, and denote by ${\bm{\beta}}_{\leq k}:=\sum_{i\leq k}\langle{\bm{\beta}},{\bm{v}}_{i}\rangle{\bm{v}}_{i}$ the orthogonal projection of ${\bm{\beta}}$ onto the span of ${\bm{v}}_{1},\dots,{\bm{v}}_{k}$, and by ${\bm{\beta}}_{>k}:={\bm{\beta}}-{\bm{\beta}}_{\leq k}$ its complement. Finally, let $k_{\star}:=\max\{k\,:\;\sigma_{k}\geq\lambda_{\star}\}$, and define the tail effective rank parameters by

Then, defining $b_{k}:=\sigma_{k}/\sigma_{k+1}$, we have $2n\geq k_{\star}+r_{1}(k_{\star})/b_{k_{\star}}$ and

(We notice that if the singular values $\sigma_{k}$ do not decay faster than exponentially, then $b_{k}$ is of order one.) While these are only bounds on the theoretical characterization $\mathsf{B}_{n}(\lambda)$, $\mathsf{V}_{n}(\lambda)$ for bias and variance, our main resuls (Theorem 1and Theorem 3) will allow to transfer them to the actual bias and variance $\mathscr{B}_{\bm{X}}(\lambda)$, $\mathscr{V}_{\bm{X}}(\lambda)$ (modulo additional error terms).

These bounds (more precisely, the bounds on $\mathscr{B}_{\bm{X}}(\lambda)$, $\mathscr{V}_{\bm{X}}(\lambda)$ that follow from these and Theorem 1) are closely related to the ones in [BLLT20, TB20], see in particular [TB20, Theorem 1]. It is worth pointing out two important differences. First, the bounds in Eqs. (14), (15) are somewhat more precise/explicit: there is no unspecified constant factorThe factor $c_{\star}$ is explicit and, if useful, can be replaced by the original expression., no dependence on the condition number of $\sigma_{1}/\sigma_{k_{\star}}$, and no multiplicative factor depending on the probability. Second, Eqs. (14), (15) are only proved for the specific value of $k_{\star}$ defined there.

The bounds of Eqs. (14), (15) allow to characterize settings in which the excess test error (as predicted by our theory) vanishes. Indeed, for $\mathsf{V}_{n}(\lambda)$ to vanish, it is sufficient that $k_{\star}/n\to 0$ and $\overline{r}(k_{\star})/n\to\infty$. A simple sufficient condition for $\mathsf{B}_{n}(\lambda)\to 0$ is that ${\bm{\beta}}\in{\rm span}({\bm{v}}_{1},\dots,{\bm{v}}_{k})$ with $\sigma_{k}/\sigma_{k_{\star}}\to\infty$.

We will discuss special examples in Section 4, and show how our general results allow to derive more precise estimates of the risk in those cases.

Equivalent sequence model.

The discussion above relies on the assumption ${\rm{Tr}}\left({{\bm{\Sigma}}^{2}({\bm{\Sigma}}+\lambda_{\star}{\bm{I}})^{-2}}\right)\leq n(1-c_{\star}^{-1})$, which implies the simple bounds (11), (12). However the interpretation in terms of a modified ridge regression problem holds for the exact formulas of Eqs. (6), (7). This interpretation was developed in the context of earlier work on the proportional asymptotics [DJM13], but it is useful to spell it out here for the present context.

In the modified model, we observe ${\bm{y}}^{s}$ that is related to ${\bm{\beta}}$ according to

Without loss of generality, we can work in the basis in which ${\bm{\Sigma}}$ is diagonal, and therefore rewrite the above as $y^{s}_{i}=\sigma_{i}^{1/2}\beta_{i}+(\omega/\sqrt{n})g_{i}$, which coincides with the definition of the classical sequence model [Tsy09].

We use ridge regression at regularization level $\lambda_{\star}$ as defined in Eq. (5):

Finally, choose the noise level $\omega$ to be the unique positive solution of

Then our theoretical prediction for the excess test error $\mathsf{R}_{n}(\lambda)$ coincides with the excess test error of the sequence model:

Summarizing, the predicted test error for the original model is equal to the test error in the sequence model, albeit at a different value of the ridge regularization parameter and of the noise level. Needless to say, studying the sequence model is significantly simpler than the original model (3).

Statement of main results

For two functions $f({\bm{x}})$ and $g({\bm{x}})$ (where ${\bm{x}}$ can be a scalar or a vector), we write $f({\bm{x}})=\mathcal{O}_{{\bm{\alpha}}}(g({\bm{x}}))$ if there exists a constant $\mathsf{C}_{{\bm{\alpha}}}$ depending only on the value of ${\bm{\alpha}}$ (also ${\bm{\alpha}}$ can be either a scalar or a vector) such that $|f({\bm{x}})|\leq\mathsf{C}_{{\bm{\alpha}}}|g({\bm{x}})|$ for all ${\bm{x}}$. In particular, if the constant is universal we write $f({\bm{x}})=\mathcal{O}(g({\bm{x}}))$. Similarly, we write $f({\bm{x}})=\Omega_{{\bm{\alpha}}}(g({\bm{x}}))$ if $|f({\bm{x}})|\geq\mathsf{C}_{{\bm{\alpha}}}|g({\bm{x}})|$ for all ${\bm{x}}$ and some constant $\mathsf{C}_{{\bm{\alpha}}}>0$. Finally, we write $f({\bm{x}})=\Theta_{{\bm{\alpha}}}(g({\bm{x}}))$ if we have both $f({\bm{x}})=\mathcal{O}_{{\bm{\alpha}}}(g({\bm{x}}))$ and $f({\bm{x}})=\Omega_{{\bm{\alpha}}}(g({\bm{x}}))$.

We wil state three theorems: the first one for ridge regression with positive regularization $\lambda>0$ (Theorem 1), and the other two for the ridgeless case $\lambda=0+$ (Theorem 2 for the overparametrized regime, and Theorem 3 for the underparametrized one). Our approximation guarantees will depend on the pair ${\bm{\Sigma}}$, ${\bm{\beta}}$ through the following three quantities (in the case $\lambda=0+$, these quantities will be modified as described below):

The ratio between effective dimension and regularization parameter:

Here $\eta$ is a constant that only depends on $\mathsf{C}_{{\bm{x}}}$, and hence we will leave it implicit.

The ratio between regularization and effective regularization

For a positive semi-definite operator ${\bm{Q}}$, define the modified population resolvent:

Letting ${\bm{\beta}}={\bm{\Sigma}}^{1/2}{\bm{\theta}}$, $\|{\bm{\theta}}\|<\infty$, we consider the ratio

We next present our master theorem for ridge regression: its proof is postponed to Section 6.

Under Assumption 1, for any positive integers $k$ and $D$, there exist constants $\eta=\eta(\mathsf{C}_{{\bm{x}}})\in(0,1/2)$ and $\mathsf{C}=\mathsf{C}(\mathsf{C}_{{\bm{x}}},D)>0$ such that the following hold. Define $\chi_{n}(\lambda),{\kappa},\rho(\lambda)$ as above (with $\eta=\eta(\mathsf{C}_{{\bm{x}}})$ in Eq. (20)).

then for all $n=\Omega_{k,D}(1)$, with probability $1-\mathcal{O}_{k}(n^{-D+1})$ we have:

Bias approximation. If we additionally have $\chi_{n}(\lambda)^{3}\log^{2}n\leq\mathsf{C}n{\kappa}^{4.5}\sqrt{\rho(\lambda)}$ and $\lambda kn^{-\frac{1}{k}}\leq n{\kappa}/2$, for all $n=\Omega_{k,D}(1)$, we have

The condition $\|{\bm{\beta}}\|_{{\bm{\Sigma}}^{-1}}<\infty$ in Assumption 1 amounts to requiring that the coefficients of ${\bm{\beta}}$ in the basis of eigenvectors ${\bm{v}}_{i}$ of ${\bm{\Sigma}}$ decay fast enough. Namely, it is equivalent to $\sum_{i}\langle{\bm{v}}_{i},{\bm{\beta}}\rangle^{2}/\sigma_{i}<\infty$. This condition appears to be a proof artifact and it would be interesting to relax it.

As mentioned above, the conditions on the isotropic random vectors ${\bm{z}}_{i}$ in Assumption 1 are mainly imposed to be able to apply Hanson-Wright inequality (Lemma 2.1). It is an interesting research question to analyze ridge regression for covariates which do not satisfy this inequality.

We next consider the ridgeless limit for in the overparametrized case: recall that $\widehat{{\bm{\beta}}}_{\lambda}$ coincides in this case with the minimum norm interpolator. In this case we need to modify the quantities defined above to measure the quality of our approximation. We begin by noting that Eq. (5) makes perfect sense in the case $\lambda=0$ and we have $\lim_{\lambda\downarrow 0}\lambda_{\star}(\lambda)=\lambda_{\star}(0)>0$. We then use the following definitions.

We replace $\chi_{n}(\lambda)$ of Eq. (20) by:

where ${\kappa}$ will be introduced in the theorem statement.

The quantity $\rho(\lambda)$ defined in Eq. (23) has a well defined limit as $\lambda\downarrow 0$, given by

Suppose Assumption 1 holds with $n<d$. Further assume $\sigma_{n}>0$, and let $s_{\sf min}$ be the minimum nonzero eigenvalue of the sample covariance ${\bm{X}}^{\sf T}{\bm{X}}/n$. For any positive integers $k$ and $D$, there exist constants $\eta=\eta(\mathsf{C}_{{\bm{x}}})\in(0,1/2)$ and $\mathsf{C}_{1}=\mathsf{C}_{1}(\mathsf{C}_{{\bm{x}}},D)>0$, $\mathsf{C}_{i}=\mathsf{C}_{i}(k,\mathsf{C}_{{\bm{x}}},D)>0$, $i\in\{2,3\}$, such that the following hold. Define $\chi_{n}^{\prime}({\kappa})$, $\rho(0)$, $\mathsf{C}_{{\bm{\Sigma}}}$ as above.

Let $\kappa>0$ be such that the following hold

Then, on the event $\{s_{\sf min}\geq 8\lambda_{\star}(0)\kappa\}$, the following hold with probability $1-\mathcal{O}_{k}(n^{-D+1})$:

Variance approximation. If in addition $\chi_{n}^{\prime}({\kappa})^{3}\log^{2}n\leq\mathsf{C}_{2}n^{1-\frac{1}{k}}{\kappa}^{9.5}$, then

Bias approximation. If in addition $\chi_{n}^{\prime}({\kappa})^{3}\log^{2}n\leq\mathsf{C}_{1}n{\kappa}^{4.5}\sqrt{\rho(0)}$, $\lambda_{\star}(0)kn^{-\frac{1}{k}}\leq 1/4$ and

The proof of this theorem is presented in Appendix D.

Our approach to proving Theorem 2 consists in reducing the ridgeless case $\lambda=0+$ to the case $\lambda>0$, and appealing to Theorem 1. For instance, when controlling the variance, we will use triangular inequality

We then use Theorem 1 to bound the first term by a quantity that diverges as $\lambda\downarrow 0$, and the main technical challenge is in bounding the other two terms by a quantity that vanishes faster than any polynomial as $\lambda\downarrow 0$.

In Theorem 2 we use the (random) minimum nonzero eigenvalue $s_{\sf min}$ of the sample covariance ${\bm{X}}^{\sf T}{\bm{X}}/n$. To apply the theorem, we need to choose $\kappa$ such that $\{s_{\sf min}\geq 8\lambda_{\star}(0)\kappa\}$ holds with high probability, and therefore we need a lower bound on $s_{\sf min}$ that holds with high probability. In this paper, we will provide such lower bonds in two cases: $(i)$ proportional regime and $(ii)$ bounded varying spectrum, cf. Section 4. In proportional regime, $s_{\sf min}=\Omega_{|d/n-1|^{-1}}(1)$ by Bai-Yin law; and for for bounded varying spectrum, $s_{\sf min}=\Omega(\sigma_{n})$ (cf. Lemma G.1).

Beyond the two cases in the paper, it would be interesting to apply Theorem 2 with results lower bounding $s_{\sf min}$ for other examples (e.g. for the kernel random matrices [HLCH19]).

In the underparameterized regime $d<n$, we have $\lim_{\lambda\downarrow 0}\lambda_{\star}(\lambda)=0$ and therefore the previous bounds do not apply. In this case, we trivially have $\mathscr{B}_{\bm{X}}(0)=\mathsf{B}_{n}(0)=0$. The proof for the variance approximation requires a different proof, which is presented in Appendix E.

Suppose Assumption 1 holds with $n>d$, and further assume

For any positive integers $k$ and $D$, there exist constants $\eta=\eta(\mathsf{C}_{{\bm{x}}})>0$ $\mathsf{C}_{1}=\mathsf{C}_{1}(\mathsf{C}_{{\bm{x}}},D)>0$, $\mathsf{C}_{2}=\mathsf{C}_{2}(k,\mathsf{C}_{{\bm{x}}},D)>0$, such that the following hold.

If ${\bm{X}}$ has rank $d$ and $s_{\sf min}$ is the minimum eigenvalue of the sample covariance ${\bm{X}}^{\sf T}{\bm{X}}/n$, then the following hold:

Variance approximation. Let $\varepsilon$ be such that

Then, on the event $\{s_{\sf min}\geq 2\varepsilon\}$, with probability $1-\mathcal{O}_{k}(n^{-D+1})$:

Bias approximation. $\mathscr{B}_{\bm{X}}(0)=\mathsf{B}_{n}(0)=0$ (this holds deterministically on the event ${\rm rank}({\bm{X}})=d$).

Applications

As a first application, we revisit the proportional regime that is defined by the following condition.

There exists a constant $M>1$ such that $M^{-1}\leq d/n\leq M$ and $\sigma_{d}\geq M^{-1}$.

This case is well studied and is not the main motivation of the present paper, but it is nevertheless important to compare our results to earlier work. We refer the reader to [Dic16, ASS20, DW18, WX20, RMR21] for background.

Among others, the results of [HMRT22] are more directly comparable to ours because they establish nonasymptotic bounds comparing variance and bias to the effective variance and bias of Eqs. (6) and (7), for both ridge and ridgeless regression. The proofs of [HMRT22] build on recent advances in random matrix theory, and in particular the anisotropic local law of [KY17].

Here we apply Theorems 1, 2 and 3 to the proportional regime. We note that, under assumption 2, the minimum eigenvalue of ${\bm{X}}^{\sf T}{\bm{X}}$ is, with high probability, of order $n$. In order for the ridge regularization to have a non-trivial effect, we need to choose $\lambda\asymp n$ as well, cf. (4a) and (4b). We will therefore assume $\lambda/n$ bounded above and below (there is no loss of generality in using the same constant as in Eq. (2)). We will address the case $\lambda=0+$ in a separate statement below.

Let Assumptions 1 and 2 hold, and further assume $\lambda/n\in[1/M,M]$. Then for any positive integers $k$ and $D$, if $n=\Omega_{k,M,\mathsf{C}_{{\bm{x}}},D}(1)$, with probability $1-\mathcal{O}_{k}(n^{-D+1})$ we have

The proof of this result is presented in Appendix F.

We note that the rates $\mathcal{O}(n^{-1})$ and $\mathcal{O}(n^{-1/2})$ are optimal for variance and bias approximation—corresponding to fluctuations of the average law and local law for the resolvent [AEK+14, KY17]. The most direct comparison is with [HMRT22, Theorem 5]: let us point out two ways in which the present result improves over the earlier [HMRT22].

In [HMRT22], the rate for variance approximation of ridge regression is $\mathcal{O}(n^{-1/2})$, while here we obtain the faster rate $\mathcal{O}(n^{-1})$.

Error terms in [HMRT22] are additive, while Proposition 4.1 provides multiplicative error terms: the quality of approximation does not deteriorate in the interesting case in which bias and variance become small.

Note that [HMRT22] informally claimed that $n^{-1/2}$ is the optimal rate in the above estimates. While this is correct for the bias, for the variance Proposition 4.1 yields a faster rate. As related phenomenon arises for linear eigenvalue statistics of random matrices (i.e. statistics of the form $n^{-1}\sum_{i=1}^{n}\varphi(\lambda_{i})$). While naively such statistics would have normal deviations of order $n^{-1/2}$, the actual deviations are of order $n^{-1}$ because of eigenvalues correlations [LP09].

Let Assumptions 1 and 2 hold for ${\bm{x}}_{i}={\bm{\Sigma}}^{1/2}{\bm{z}}_{i}$, where ${\bm{z}}_{i}$ has i.i.d. sub-Gaussian coordinates.

Overparameterized regime. If additionally $d/n\geq 1+M^{-1}$, then for all $n=\Omega_{M,\mathsf{C}_{{\bm{x}}},D}(1)$, with probability $1-\mathcal{O}(n^{-D+1})$ we have

Underparameterized regime. If additionally $d/n\leq 1-M^{-1}$, then for all $n=\Omega_{M,\mathsf{C}_{{\bm{x}}},D}(1)$, with probability $1-\mathcal{O}(n^{-D+1})$ we have

We do not expect the exponent $1/14$, $1/28$, $1/5$ in this statement to be tight. However, as in the positive $\lambda$ case, also in this case the error is multiplicative.

2 Bounded varying spectrum

We next consider the highly overparametrized case $d\gg n$. Overparametrized ridge (or minimum norm) regression attracted significant attention recently because of the realization that many deep learning models are overparametrized and overfit the training data. This connection is reviewed in [BMR21, Bel21].

Here we consider covariate vectors ${\bm{x}}_{i}$ taking values in a general Hilbert space with $d=\infty$, under Assumption 1 on the covariates distribution. This is most closely related to [BLLT20, TB20], and [KZSS21]. The last paper derives refined upper bounds using Gaussian width techniques, but is limited to the case of Gaussian covariates and, as for earlier results, is only accurate up to constant factors.

We will show that our general theory yields excess risk estimates that are accurate up to $1+o_{n}(1)$ multiplicative errors. We impose the following condition on the spectrum of ${\bm{\Sigma}}$.

Recall that, by definition, for any $j\leq i$, $\sigma_{j}/\sigma_{i}\geq 1$. The bounded varying condition requires that, if $i,j$ diverge proportionally, then the eigenvalue ratio $\sigma_{j}/\sigma_{i}$ stays bounded. Note that this assumption is equivalent to $\sup_{i\geq 1}\sigma_{\lfloor\delta i\rfloor}/\sigma_{i}<\infty$ for every $\delta\in(0,1]$, which is in turn equivalent to

As special case, Assumption 3 holds if the sorted eigenvalues $(\sigma_{1},\sigma_{2},\cdots)$ forms a so-called regularly varying sequence, namely for any $\delta\in(0,\infty)$,

where $\psi(\delta)$ is positive and finite for any $\delta$. In other words, in the regularly varying case, the ratio $\sigma_{j}/\sigma_{i}$ converges when $i,j$ diverge proportionally.

A special case of regularly varying spectrum is given by Zipf’s law whereby $\sigma_{i}=i^{-\alpha}$ for some $\alpha>1$ (in this case $\psi(\delta)=\delta^{-\alpha}$). Regularly varying functions were characterized by Karamata [Kar33] (for functions on the positive real line), and by Galambos and Seneta [GS73] (for the sequences, i.e. functions defined on the naturals). Namely all such sequences take the form

where $a_{i}$ are arbitrary and converge to a positive limit as $i\to\infty$ and $b_{i}\to 0$.

It is easy to see that Assumption 3 holds beyond the case of regularly varying sequences. Consider for instance $\sigma_{i}=3^{-s}$ for all $2^{s}\leq i<2^{s+1}$, $s=0,1,\cdots$.

Applying Theorems 1 and 2 to ${\bm{\Sigma}}$ with bounded varying spectrum, we obtain the following result, whose proofs are detailed in Appendix G.

Let Assumptions 1 and 3 hold. For any constants $M>0$, $\gamma\in(0,1/3)$, and positive integers $k$, $D$ the following holds. If $\mathsf{d}_{{\bm{\Sigma}}}\leq Mn^{1+\gamma}$ and $\lambda\in[n\lambda_{\star}(0)/M,n\lambda_{\star}(0)M]$, then for $n=\Omega_{k,M,\psi,\gamma,\mathsf{C}_{{\bm{x}}},D}(1)$, with probability $1-\mathcal{O}_{k}(n^{-D+1})$

If additionally $\mathsf{d}_{{\bm{\Sigma}}}=\mathcal{O}_{M,\psi,\mathsf{C}_{{\bm{x}}}}(n^{1+\gamma}\left({\rho(\lambda)}\right)^{1/6})$ $\lambda\in[n\lambda_{\star}(0)/M,n\lambda_{\star}(0)M]$ and $\lambda_{\star}(0)=\mathcal{O}(1)$ (cf. Theorem 1 for the function $\rho$), with the same probability we have

Applying Theorem 2, we have the following conclusion for ridgeless regression.

Let Assumptions 1 and 3 hold. Suppose ${\bm{\beta}}={\bm{\Sigma}}^{1/2}{\bm{\theta}}$ with $\left\|{{\bm{\theta}}}\right\|<\infty$. If we have $\lambda_{\star}(0)/\sigma_{n}=\mathcal{O}(\log^{\mathcal{O}(1)}n)$ and $\mathsf{d}_{{\bm{\Sigma}}}(n)=\mathcal{O}(n\log^{\mathcal{O}(1)}n))$, for any $n=\Omega_{\psi,\mathsf{C}_{{\bm{x}}},D}(1)$, it holds with probability $1-\mathcal{O}(n^{-D+1})$ that

The assumptions $\lambda_{\star}(0)/\sigma_{n}=\mathcal{O}(\log^{\mathcal{O}(1)}n)$ and $\mathsf{d}_{{\bm{\Sigma}}}=\mathcal{O}(n\log^{\mathcal{O}(1)}n)$ are primarily introduced to simplify the form of the statement. These two conditions can be relaxed to $lambdaefct(0)/\sigma_{n}=\mathcal{O}(n^{\bar{\gamma}})$ and $\mathsf{d}_{{\bm{\Sigma}}}=\mathcal{O}(n^{1+\bar{\gamma}})$ for a sufficiently small $\bar{\gamma}$, but we do not pursue this generalization here.

It is possible to apply the upper/lower bounds on the bias of Theorem 2 to prove bounds on the bias in the setting of Proposition 4.4. However the resulting error term is larger than $(\sigma_{n}^{2}\left\|{{\bm{\theta}}_{\leq n}}\right\|^{2}+\sigma_{n}\left\|{{\bm{\beta}}_{>n}}\right\|^{2})$, which is the size of the upper bound on $\mathsf{B}_{n}(0)$ in Proposition 2.2.

In order to illustrate the accuracy of our general framework, we apply Proposition 4.3 to derive sharp asymptotics for bias and variance in a number cases. In each of the case below, we scale the regularization parameter $\lambda$ as $\lambda=\lambda_{0}(n)\cdot\nu$ for a certain explicit function $\lambda_{0}(n)$. The scaling $\lambda_{0}(n)$ is chosen so that the bias and variance retain a non-trivial dependence on $\nu$ for large $n$. We expect that the excess risk achieved by optimal regularization is also covered by this scaling (up to negligible corrections), but do not prove it formally here.

Let Assumption 1 hold. Then, for a fixed constant $\nu>0$ and any positive integer $D$, the following events hold with probability $1-\mathcal{O}(n^{-D})$ (the $o_{n}(1)$ errors may depend on $D$):

Regularly varying spectrum with $\alpha>1$. Assume $(\sigma_{i})_{i\geq 1}$ is a regularly varying sequence with exponent $\alpha>1$. As a consequence, $\sigma_{i}=i^{-\alpha}a_{i}\exp\left\{{\sum_{j=1}^{i}b_{j}/j}\right\}$ with $a_{i}$ converging to a positive limit and $b_{i}\to 0$. Define ${\sf c}_{\star}={\sf c}_{\star}(\nu)>0$ as the unique positive solution of

Let $F_{\bm{\beta}}(x)=\sum_{k=1}^{\lfloor nx\rfloor}\langle{\bm{\beta}},{\bm{v}}_{k}\rangle^{2}$. If additionally ${\bm{\beta}}$ satisfies the following “polynomial-decay” property: for some $0<\theta\leq 1$ that

Regularly varying spectrum with $\alpha=1$. Next consider the case $\sigma_{i}=i^{-1}a_{i}(1+\log i)^{-\alpha^{\prime}}$ for some $\alpha^{\prime}>1$ with $a_{i}$ converging to a positive limit. Define ${\sf c}_{\star}={\sf c}_{\star}(\nu)>0$ as

Let $F_{\bm{\beta}}(x)=\sum_{k=1}^{\lfloor(n/\log n)x\rfloor}\langle{\bm{\beta}},{\bm{v}}_{k}\rangle^{2}$. If additionally ${\bm{\beta}}$ satisfies the following “rapid-decay” property: for some $0<\theta\leq 1$ that

A non-regularly varying spectrum. $\sigma_{i}=p^{-s}$ for all $q^{s}\leq i<q^{s+1}$, with $1<q<p$ and $s=0,1,\ldots$ Define $s_{\star}$ such that $q^{s_{\star}}\leq n<q^{s_{\star}+1}$, and for positive integer $r$ the following decreasing function in $t>0$,

Let $\rho_{\star}=n/(q^{s_{\star}+1}-q^{s_{\star}})\in[1/(q-1),q/(q-1))$. Then there exists a unique solution ${\sf c}_{\star}={\sf c}_{\star}(\nu)$ to the following equation

Let $F_{{\bm{\beta}}}(x)=\sum_{k=1}^{q^{\lceil x\rceil}-1}\langle{\bm{\beta}},{\bm{v}}_{k}\rangle^{2}$. If additionally ${\bm{\beta}}$ satisfies the following “rapid-decay” property: for some $0<\theta\leq 1$ that

The proof of this theorem is presented in Appendix H.

In the case of a regularly varying spectrum with $\alpha>1$, the bias vanishes with the sample size as $n^{-\alpha+o(1)}$ but the variance stays bounded away from zero as long as $\tau>0$, cf. Eq. (27). In other words in this case overfitting is not benign and Theorem 4 quantifies precisely this claim.

On the other hand, in the case $\alpha=1$, both bias and variance vanish for large $n$, an therefore we achieve benign overfitting. We must emphasize however that the variance decay is very slow, namely $\mathscr{V}_{\bm{X}}(\lambda)\asymp(\log n)^{-1}$, and hence the decay of the excess risk is at least as slow.

Numerical illustrations

In this section we evaluate numerically the theoretical prediction for variance and bias, cf. Eqs. (6), (7) and compare them with the results of numerical simulations with synthetic data. We carry out the simulations in the ridgeless limit $\lambda=0+$ (corresponding to min-norm interpolation). This case is interesting because it is not covered by some of our theorems. Our numerical experiments suggest that the theoretical predictions of Eqs. (6), (7) hold in a broader domain of validity than the one that we are able to control rigorously.

We use Gaussian covariates ${\bm{x}}_{i}$. By rotational invariance, we can limit ourselves to diagonal covariance ${\bm{\Sigma}}$. We will consider two eigenvalue structures:

This is defined by $\sigma_{i}=i^{-\alpha}$ for all $i\geq 1$. This fits within the first case of Theorem 4.

This model is defined by $\sigma_{i}=i^{-1}(1+\log i)^{-\alpha^{\prime}}$, with $\alpha^{\prime}>1$. This fits within the second case of Theorem 4.

In all numerical experiments, we generate data according to the model (3) with a true parameters vector ${\bm{\beta}}$ concentrated on the top $d_{0}=100$ eigenvectors of ${\bm{\Sigma}}$. More precisely, we will use ${\bm{\beta}}=(1,1,\dots,1,0,0,\ldots)$ where $\|{\bm{\beta}}\|_{0}=d_{0}=100$.

In Figure 1, we plot our theoretical predictions $\mathsf{V}_{n}$, $\mathsf{B}_{n}$, $\mathsf{R}_{n}$ for variance, bias and as a function of the sample size $n$, for the two models $(I)$ and $(II)$ defined above. We use $\lambda=0+$. In each case, we consider several values of the exponents $\alpha$, $\alpha^{\prime}$ that control the decay of eigenvalues of ${\bm{\Sigma}}$.

In Figure 2, we plot the same quantities at fixed sample size $n=500$ and vary the regularization parameter $\lambda$. A few facts emerge from these figures:

For both models, the bias of the minimum norm interpolator is a decreasing function of the sample size $n$, and appears to vanish as $n\to\infty$, see second row of Figure 1.

In contrast, the variance exhibits a strikingly different behavior in the two covariance models, see first row of Figure 1. For model $(I)$ (polynomial eigenvalue decay, with exponent $\alpha>1$), the variance increases with $n$, and eventually stabilizes to a limit value. For model $(II)$ (exponent $\alpha=1$), the variance decreases with $n$, and appears to vanish, albeit very slowly, as $n\to\infty$.

As a consequence of these points, the excess test error of minimum norm interpolation vanishes with sample size in model $(II)$ but does not vanish in model $(I)$. This behavior (and the one at previous points) is precisely quantified by Theorem 4 for $\lambda>0$.

Finally the dependence of bias and variance on $\lambda$ is the expected one. As $\lambda$ increases, bias increases but variance decreases. However, the balance between these two factors is non-trivial:

For the slowest eigenvalue decay (large $\alpha$ in model $(I)$ or large $\alpha^{\prime}$ in model $(II)$), the optimal $\lambda$ is strictly positive.

On the other hand, for the fastest eigenvalue decay, the optimal $\lambda$ vanishes. In these case interpolation is superior to ridge regression: we need to overfit to achieve the best test error.

The above discussion is based on evaluating the theoretical formulas for bias and variance, as given in Eqs. (6), (7). While our main result, Theorems 1, 2 guarantee that these formulas are accurate, it is important how accurate they are at small or moderate $n$, and whether random deviations modify the picture.

In Figure 3 we plot numerical simulations corroborating that $\mathscr{V}_{\bm{X}},\mathscr{B}_{\bm{X}}$ do concentrate around $\mathsf{V}_{n},\mathsf{B}_{n}$ in models $(I)$ and $(II)$. As mentioned above, the predictions $\mathsf{V}_{n},\mathsf{B}_{n}$ appear to be accurate beyond what is guaranteed by Theorem 2, and the error appears to be a $(1+o_{n}(1))$ multiplicative factor.

Proof of Theorem 1

Let $\mathcal{F}_{k}:=\sigma({\bm{x}}_{1},\cdots,{\bm{x}}_{k})$ be the $\sigma$-field generated by the first $k$ data points for $1\leq k\leq n$, and $\mathcal{F}_{0}$ the trivial $\sigma$-field. We then have $\mathscr{V}_{\bm{X}},\mathscr{B}_{\bm{X}}\in\mathcal{F}_{n}$ and $\mathsf{V}_{n},\mathsf{B}_{n}\in\mathcal{F}_{0}$. Extending the previous notation of $\mathscr{R}_{0}$ in Eq. (22) to $\mathscr{R}_{k}$, we let

For $\mu=0$, this equation reduces to Eq. (5), via the change of variables $\mu_{\star}=\lambda/\lambda_{\star}$. For $\mu>0$ existence and uniqueness follows by a similar argument to the case $\mu=0$. Indeed, setting $\xi:=(\mu_{\star}-\mu)^{-1}$, the equation is equivalent to $n\xi=1+{\rm{Tr}}({\bm{\Sigma}}({\bm{A}}+\xi{\bm{\Sigma}}))$, where ${\bm{A}}:=\lambda{\bm{I}}+\mu{\bm{\Sigma}}$. Existence and uniqueness follow since the left-hand side is monotone increasing and the right-hand side monotone decreasing in $\xi$.

In order to quantify the approximation errors $|\mathscr{V}_{\bm{X}}-\mathsf{V}_{n}|$ and $|\mathscr{B}_{\bm{X}}-\mathsf{B}_{n}|$, we will apply the following lemma (Lemma 6.1), which expresses the bias and variance $\mathscr{B}_{\bm{X}}$, $\mathscr{V}_{\bm{X}}$ in terms of derivatives of $\mathscr{F}_{n}$ and $\mathscr{F}_{0}$ w.r.t. $\lambda$ and $\mu$.

For any $\lambda>0,\mu\geq 0$, the quantity $\mu_{\star}>\mu$ is uniquely determined and we have

and $\mu_{\star}(\lambda,0)=\lambda/\lambda_{\star}$.

Our proof strategy proceeds in four parts: (I) We show that —due to the regularity properties of $\mathscr{F}_{0}$ and $\mathscr{F}_{n}$— a bound on $|\mathscr{F}_{0}-\mathscr{F}_{n}|$ implies a bound on the difference of their derivatives, and hence (via Lemma 6.1) on the error in approximating bias and variance; (II) We prove a bound on $|\mathscr{F}_{0}-\mathscr{F}_{n}|$ interpolating between $\mathscr{F}_{0}$ and $\mathscr{F}_{n}$ by adding one row at the time to ${\bm{X}}$; (III), (IV) We apply these general bounds respectively to controlling variance and bias.

Recall that we defined ${\bm{\theta}}:={\bm{\Sigma}}^{-1/2}{\bm{\beta}}$, and assumed $\|{\bm{\theta}}\|<\infty$. By homogeneity, we can and will assume $\left\|{{\bm{\theta}}}\right\|=1$ throughout the proof.

The following lemma reduces controlling the difference of derivatives of $\mathscr{F}_{0}$ and $\mathscr{F}_{n}$ to the less arduous task of bounding the difference in function values. Its proof is presented in Appendix B.2.

Before passing to bounding errors in function values, we provide upper bounds for higher order derivatives in Eqs. (37) and (38). Bounding the derivatives of $\mathscr{F}_{n}$ is easier as we can easily write an explicit formula for the $k$-th derivative for any $k$. (The proof of this lemma is presented in Appendix B.3).

Computing higher order derivatives of $\mathscr{F}_{0}$ is less straightforward because $\mathscr{F}_{0}$ depends on $\mu_{\star}$ which itself depends implicitly depending on $(\lambda,\mu)$. We postpone this proof to Appendix B.4.

and for all $\mu$ such that $0\leq\mu\leq\mu_{\star}(\lambda,\mu)/2$,

Part II: Bounding errors in function values.

We next proceed to bounding $|\mathscr{F}_{n}(\lambda,\mu;{\bm{Q}})-\mathscr{F}_{0}(\lambda,\mu_{\star}(\lambda,\mu);{\bm{Q}})|$ for a p.s.d. matrix ${\bm{Q}}$, which appears in Eqs. (37) and (38). Recall that $\mathscr{F}_{i}(\lambda,\mu;{\bm{Q}})=\lambda\mathscr{R}_{i}(\lambda,\mu;{\bm{Q}})$.

The next theorem bounds $|\mathscr{R}_{n}(\lambda,\mu;{\bm{Q}})-\mathscr{R}_{0}(\lambda,\mu_{\star}(\lambda,\mu);{\bm{Q}})|$ and is the most importan technical step in the proof of our main theorems. Its proof is outlined in Section 7, with several technical lemmas deferred to the appendices

Introduce the shorthand $\mathsf{R}_{0}({\bm{Q}}):=\mathscr{R}_{0}(\lambda,\mu_{\star}(\lambda,\mu);{\bm{Q}})$. Under Assumption 1, for any $\lambda>0,\mu\geq 0$, p.s.d. matrix ${\bm{Q}}$ with $\|{\bm{Q}}\|=1$ and positive integer $D$, there exists constants $\eta=\eta(\mathsf{C}_{{\bm{x}}})\in(0,1/2)$, $\mathsf{C}_{\alpha}=\mathsf{C}_{\alpha}(\mathsf{C}_{{\bm{x}}},D)>0$, $\mathsf{C}_{\beta}=\mathsf{C}_{\beta}(\mathsf{C}_{{\bm{x}}},D)>0$ and $\mathsf{C}_{\gamma}=\mathsf{C}_{\gamma}(\mathsf{C}_{{\bm{x}}},D)$ such that for

if $\alpha_{1}\leq\mathsf{R}_{0}({\bm{I}})/8$, $\beta_{1}\leq\mathsf{R}_{0}({\bm{Q}})/64$, $\gamma\beta_{2}(1+\mathsf{R}_{0}({\bm{I}}))\leq 1/64$ and $n^{-D}=\mathcal{O}(\alpha_{1}/(1+\mathsf{R}_{0}({\bm{I}})))$, for all $n=\Omega_{D}(1)$ with probability $1-\mathcal{O}(n^{-D+1})$ we have

Let us emphasize that this theorem holds under weaker assumptions than Theorem 1, but the error bounds it provides are quite implicit. We can obtain more explicit bounds by imposing the assumptions of Theorem 1. We first define the generalized version of $\rho(\lambda)$ in Eq. (23) for any p.s.d. matrix ${\bm{Q}}$ as

The proof of this corollary is given in Appendix C.6.

Under Assumption 1, for any positive integers $k$, $D$ and p.s.d. matrix ${\bm{Q}}$ with $\|{\bm{Q}}\|=1$, there exist constants $\eta=\eta(\mathsf{C}_{{\bm{x}}})\in(0,1/2)$ and $\mathsf{C}=\mathsf{C}(\mathsf{C}_{{\bm{x}}},D)>0$, such that the following hold. Define $\chi_{n}(\lambda),{\kappa},\rho(\lambda)$ as per Eqs. (20), (21), (39) (those quantities are defined for $\mu=0$). If it holds that $\mu_{\star}(\lambda,\mu)\leq(1-{\kappa}/2)^{-1}\mu_{\star}(\lambda,0)$, and

we then have for all $n=\Omega_{D}(1)$ with probability $1-\mathcal{O}(n^{-D+1})$ that

To further simplify the assumption $\mu_{\star}(\lambda,\mu)\leq(1-{\kappa}/2)^{-1}\mu_{\star}(\lambda,0)$ in Corollary 6.5, the next lemma will be helpful. We defer its proof to Appendix B.5.

For any fixed $\lambda>0$, the function $\mu_{\star}(\lambda,\mu)$ is increasing in $\mu$ for all $\mu\geq 0$. Assuming Eq. (21), if $0\leq\mu\leq n{\kappa}^{3}/2$, then

Part III: Approximation error for variance.

We are now ready to combine our results in Part I and Part II to obtain approximation errors $|\mathscr{V}_{\bm{X}}-\mathsf{V}_{n}|$ and $|\mathscr{B}_{\bm{X}}-\mathsf{B}_{n}|$. For the variance, we want to take $\delta$ in Eq. (37) such that $k\delta\leq\lambda$. In this case, $[\lambda,\lambda+k\delta]\subset[\lambda,2\lambda]$. Note that $\lambda_{\star}(\lambda)$ is an increasing function of $\lambda$. Further, by

we know $\lambda\mapsto\mu_{\star}(\lambda,0)=\lambda/\lambda_{\star}$ is an increasing function. Therefore $\lambda/\lambda_{\star}(\lambda)\leq 2\lambda/\lambda_{\star}(2\lambda)$, which implies $\lambda_{\star}(\lambda)\leq\lambda_{\star}(2\lambda)\leq 2\lambda_{\star}(\lambda)$. For $\lambda$ that satisfies Eq. (21), this guarantees that for any $\lambda^{\prime}\in[\lambda,2\lambda]$,

Hence, for any $\lambda^{\prime}\in[\lambda,2\lambda]$, Eq. (21) still holds but with constant ${\kappa}^{\prime}\geq{\kappa}/2$. Therefore, we can apply Corollary 6.5 for any $\lambda^{\prime}\in[\lambda,2\lambda]$ for ${\bm{Q}}={\bm{I}}$ and $\mu=0$, provided the following conditions hold

where the last equality used the fact that $\rho(\lambda^{\prime})=1$ when ${\bm{Q}}={\bm{I}}$. Finally, setting $\mathsf{C}^{\prime}:=2^{-4.5}\mathsf{C}$ and using the fact that $\chi_{n}(\lambda^{\prime})$ is decreasing in $\lambda$, it suffices to require

which holds by the theorem’s assumptions.

Hence, we can now apply Corollary 6.5 with ${\bm{Q}}={\bm{I}}$, and it follows that with probability $1-\mathcal{O}_{k}(n^{-D+1})$,

where in the last inequality we use that $\mathscr{R}_{0}(\lambda,\mu_{\star}(\lambda,0);{\bm{I}})=n/\mu_{\star}(\lambda,0)-1$ is a decreasing function in $\lambda$ as $\mu_{\star}(\lambda,0)$ is increasing in $\lambda$. Next by Lemmas 6.3 and 6.4 we obtain

where in (i) we use again that $\mathscr{R}_{n}(\lambda,\mu;{\bm{I}})$ and $\mathscr{R}_{0}(\lambda,\mu_{\star}(\lambda,0);{\bm{I}})$ are decreasing in $\lambda$ and in (ii) we apply Corollary 6.5. Substituting the above displays into Eq. (37), we have

We therefore have, by Lemma 6.1, $\tau^{2}\mathscr{F}_{0}(\lambda,\mu_{\star}(\lambda,0);{\bm{I}})\leq\lambda{\kappa}^{-1}\mathsf{V}_{n}(\lambda)$. Substituting in Eq. (41), we obtain

By setting $\delta=\lambda{\kappa}^{2}n^{-1/k}$, the condition $\delta k\leq\lambda$ is satisfied for all $n=\Omega_{k}(1)$, which completes the proof for variance approximation with

where we use $\chi_{n}(\lambda)\geq 1$ in the final bound.

Part IV: Approximation error for bias.

Note that all the terms on the right-hand side of Eq. (38) are evaluated at the same value of $\lambda$. Hence, Eq. (21) applies to each of these terms. We claim that the assumptions of Corollary 6.5 apply to all of these terms, provided the following conditions hold

Indeed, condition (42) implies $\mu_{\star}(\lambda,\mu)\leq(1-{\kappa}/2)^{-1}\mu_{\star}(\lambda,0)$ for all $\mu\in[0,k\delta]$ since $\mu\mapsto\mu_{\star}(\lambda,\mu)$ is monotone decreasing; finally, condition (43) is independent of $\mu$.

then we can apply Lemmas 6.3 and 6.4 and invoke Corollary 6.5 with ${\bm{Q}}={\bm{\theta}}{\bm{\theta}}^{\sf T}$. To be specific, by Corollary 6.5, we have with probability $1-\mathcal{O}_{k}(n^{-D+1})$,

as $\mathscr{F}_{0}(\lambda,\mu_{\star}(\lambda,\mu);{\bm{\theta}}{\bm{\theta}}^{\sf T})$ decreases with $\mu$. By Lemmas 6.3 and 6.4 we obtain

where in the bound (i) we use the fact that $\mu_{\star}(\lambda,\mu)$ is increasing in $\mu$ (cf. Lemma 6.6) and $\mathscr{F}_{k}(\lambda,\mu;{\bm{Q}})$ is decreasing in $\mu$ when $\mu\geq 0$; in (ii) we use that $\mu_{\star}(\lambda,0)=\lambda/\lambda_{\star}$. Combining the calculations above, we have from Eq. (38)

where in the last line we used the definition of $\mathsf{B}_{n}(\lambda)$ in Eq. (7). By Eq. (21), we have

which reduces the approximation bound for bias to

We again take $\delta=\lambda{\kappa}^{2}n^{-\frac{1}{k}}$ and the bound becomes

This bounds hold under the conditions (42) to (43), which are implied by the following:

For the first condition, we invoke Lemma 6.6 to obtain a sufficient requirement $\lambda kn^{-\frac{1}{k}}\leq n{\kappa}/2$. For the last condition, it suffices to have $\chi_{n}(\lambda)^{3}\log^{2}n\leq\mathsf{C}^{\prime}n{\kappa}^{4.5}\sqrt{\rho(\lambda)}$ for the same $\mathsf{C}^{\prime}$ defined in Part III.

Proof of Theorem 5

The proof is based on the following interpolating construction. We will construct a sequence of random variables $\mu_{i}\in\mathcal{F}_{i-1}$ for $i=0,1,\cdots,n+1$ (where, by convention, $\mathcal{F}_{-1}=\mathcal{F}_{0}$ is the trivial $\sigma$-algebra) such that, defining

we obtain that $\mathsf{R}_{i}({\bm{Q}})$ is approximately a martingale and, as a consequence, $\mathsf{R}_{0}({\bm{Q}})\approx\mathsf{R}_{n}({\bm{Q}})$. We will further have $\mu_{0}=\mu_{\star}(\lambda,\mu)$ and $\mu_{n}\approx\mu$, asd therefore we obtain the desired claim $\mathscr{R}_{0}(\lambda,\mu_{\star}(\lambda,\mu);{\bm{Q}})\approx\mathscr{R}_{n}(\lambda,\mu;{\bm{Q}})$.

Before formally defining the sequence $\{\mu_{0},\cdots,\mu_{n+1}\}$, we introduce some helpful notations. We first define the matrices ${\bm{A}}_{i},{\bm{B}}_{i}\in\mathcal{F}_{i}$ for $0\leq i\leq n$ as

Then we can write $\mathsf{R}_{i}({\bm{Q}})={\rm{Tr}}\left({{\bm{Q}}{\bm{A}}_{i}}\right)$. Similarly we define another sequence of functions by $\mathsf{S}_{i}({\bm{Q}}):={\rm{Tr}}\left({{\bm{Q}}{\bm{B}}_{i}}\right)$.

Now we are ready to define the sequence $\mu_{i}\in\mathcal{F}_{i-1}$. We set the initial value $\mu_{0}=\mu_{\star}(\lambda,\mu)\in\mathcal{F}_{-1}$ and thus $\mathsf{R}_{0}({\bm{Q}})=\mathscr{R}_{0}(\lambda,\mu_{\star}(\lambda,\mu);{\bm{Q}})$. The sequence $(\mu_{j})_{j\geq 1}$ is iteratively determined through the following equation

It is evident that if the solution $\mu_{i+1}$ exists and is unique (almost surely with respect to the random choice of $\mu_{i}$), since $\mu_{i}\in\mathcal{F}_{i-1}$ and ${\bm{X}}_{i}\in\mathcal{F}_{i}$, it follows that $\mu_{i+1}\in\mathcal{F}_{i}$. The next lemma shows that the iteration via (46) is indeed well-defined. Its proof is in Appendix C.1.