Kernel regression in high dimensions: Refined analysis beyond double descent

Introduction

Interpolation learning [MM19, HMRT19, BLLT20] has recently attracted growing attention in the machine learning community. This is mainly because current state-of-the-art neural networks appear to be models of this type: they are able to interpolate the training data while still generalize well on test data, even in the presence of label noise [ZBH+16]. It has been empirically observed that other models including random features, decision trees, and as simple as linear regression also exhibit similar phenomenon [BLLT20, BHMM19, LHCS20]. This is somewhat striking as it goes against the conventional wisdom of bias-variance trade-off[CS02]: predictors that generalize well must trade off the model complexity against training data fitting. The double descent theory [BHMM19] resolves this paradox by revisiting the bias-variance trade-off and showing that the model generalization error exhibits a phase transition at the interpolation point: moving away from this point on either side tends to reduce the generalization error.

The double descent phenomenon has recently inspired intense theoretical research [MM19, GLK+20, WX20, LCM20] and has been further extended to multiple descent [CMBK20, LRZ19, AP20] on various models. One line of work formalized the argument that, even when no explicit regularization is imposed, implicit regularization is encoded in the model via the choice of optimization algorithms and techniques, e.g., stochastic gradient descent (SGD) [HRS16], dropout [SHK+14], early stopping [AKT19], and ensemble methods [LJB20]. Different from these “external” schemes, the kernel interpolation estimator [LR20, BRT19] directly benefits from its intrinsic kernel structure that serves as an implicit regularization to help both interpolate and approximate. In fact, (strictly) positive-definite kernels can interpolate an arbitrary number of data points [Wen04], and thus kernel spaces contain (nearly) optimal interpolants [GMMM19, Li20]. Although the kernel space is rich enough to contain models that generalize well, the generalization property of kernel method, for example how it depends on the choice of kernel, its interplay with the data and the level of regularization, still remains unclear. In particular, the question whether the double descent phenomenon exists in the kernel regression models is still unanswered [LR20, BCP20]. As such, refined analyses are needed to have a thorough understanding of kernel estimators, notably in the high dimensional regime of interest. This is indeed the objective of the article.

where an explicit Tikhonov regularization term induced by a reproducing kernel Hilbert space (RKHS) $\mathcal{H}$ is added to the least-squares objective. In statistical learning theory [CZ07], the regularization parameter $\lambda>0$ is generally taken to depend on the sample size $n$ in such a way that $\lim_{n\rightarrow\infty}\lambda(n)=0$. Here we assume that $\lambda:=\bar{c}n^{-\vartheta}$ with some $\vartheta\geq 0$ and $0\leq\bar{c}\leq 1$ to cover the interpolation case.

In this paper, we propose a novel bias-variance decomposition of the KRR expected excess risk, and derive non-asymptotic bounds for both bias and variance. This precise assessment leads to fruitful discussions as a function of different data eigenvalue decays and regularization schemes. Our main findings include:

The explicit regularization $\lambda:=\bar{c}n^{-\vartheta}$ largely affects the peak point of the variance: a large $\lambda$ decreases the model complexity, and thus corresponds to a small value of interpolation point $n_{*}\equiv n_{*}(\lambda)$. Table 1 shows that, under a small (or zero) regularization so that $r_{*}\leq n_{*}$ with $r_{*}:=\operatorname{rank}(\bm{X}\bm{X}^{\!\top}/d)$: the error bound for variance ${\tt{V}}$ monotonically increases with $n$ until $n:=r_{*}$, as in the red curve of Figure 1(a). Under a moderate regularization with $n_{*}\leq r_{*}$: ${\tt{V}}$ first increases with $n$ until $n:=n_{*}$ and then decreases. In this case, the peak point will move to the left due to $n_{*}<d$, see the blue curve in Figure 1(a). Under a large regularization with $n_{*}\leq c$ for some constant $c$, ${\tt{V}}$ monotonically decreases with $n$, as in the green curve of Figure 1(a).

Our error bounds for the bias and the variance exhibit different characteristics. More specifically, the bias bound is (almost) independent of the data/feature dimension $d$ and monotonically decreases with $n$ at a certain $\mathcal{O}(\lambda)$ (learning) rate as in the classical learning theory [CZ07, WZ11, SS07]. Besides, the variance bound depends on $n$ and $d$, and exhibits monotonic decreasing or unimodal with $n$ under different regularizations. Hence, the expected excess risk, as the sum of bias and variance, can be double descent (Figure 1(b)), bell-shaped (Figure 1(c)), or monotonic decreasing (Figure 1(d)), depending on the level of implicit and explicit regularizations. This is in agreement with empirical findings in neural networks [YYY+20].

Our non-asymptotic results show that, for large but fixed $d$, both the variance and bias tends to zero as $n\rightarrow\infty$ under $\lambda:=\bar{c}n^{-\vartheta}$, implying that the excess risk approaches zero. Based on this, in the double descent case particularly, the minimum of the expected error in the over-parameterized $n>d$ regime is lower than that in the $n<d$ regime. This claim cannot be obtained from [LR20].

The rest of the paper is organized as follows. We briefly introduce problem settings in Section 2. In Section 3, we present our main results on the generalization property of KRR in high dimensions and briefly sketch the main ideas of the proof. Discussions on the derived error bounds are given in Section 4. In Section 5, we report numerical experiments to support our theoretical results and the conclusion is drawn in Section 6.

Problem Settings and Preliminaries

We work in the high dimensional regime for some large $d,n$ with $c\leq d/n\leq C$ for some constants $c,C>0$. For notational simplicity, we denote by $a(n)\lesssim b(n)$: there exists a constant $\widetilde{C}$ independent of $n$ such that $a(n)\leq\widetilde{C}b(n)$, and analogously for $\asymp$ and $\gtrsim$.

2 Background on RKHS

Now we characterize the integral operators defined by a kernel. Given a kernel $k$, its integral operator $L_{K}:\mathcal{L}_{\rho_{X}}^{2}\rightarrow\mathcal{L}_{\rho_{X}}^{2}$ admits

Since $L_{K}$ is compact, positive definite and self-adjoint, by the spectral theorem (see, Theorem A.5.13 in [SA08]), there exists countable pairs of eigenvalues and eigenfunctions $\{\mu_{i},\psi_{i}\}_{i=1}^{\infty}$ of $L_{K}$ such that $L_{K}\psi_{i}=\mu_{i}\psi_{i}$, where $\{\psi\}_{i=1}^{\infty}$ are orthogonal basis of $\mathcal{L}_{\rho_{X}}^{2}(X)$ and $\mu_{1}\geq\mu_{2}\cdots>0$ with $\lim\limits_{i\rightarrow\infty}\mu_{i}=0$. Accordingly, by Mercer’s theorem, we have $k(\bm{x},\bm{x}^{\prime})=\sum_{i=1}^{\infty}\mu_{i}\psi_{i}(\bm{x})\psi_{i}(\bm{x}^{\prime})$, and there exists a constant $\kappa\geq 1$ such that $\sup_{\bm{x}\in X}\sum_{i=1}^{\infty}\mu_{i}\psi_{i}^{2}(\bm{x})\leq\kappa^{2}$. It holds by $\kappa:=\max\{1,\sup_{\bm{x}\in X}\sqrt{k(\bm{x},\bm{x})}\}$. Based on the data matrix $\bm{X}$ and the integral operator $L_{K}$, the empirical integral operator is given by $L_{K,\bm{X}}=\frac{1}{n}\sum_{i=1}^{n}k(\cdot,\bm{x}_{i})\otimes k(\cdot,\bm{x}_{i})$, which converges to the data-free limit $L_{K}$ at an $\mathcal{O}(1/\sqrt{n})$ rate [DMDVR09].

Main Results

In this section, we state our main result under some basic/technical assumptions, compare it with existing results, and sketch the main ideas of our proof.

To illustrate our analysis, we need the following three standard assumptions.

(Existence of $f_{\rho}$) We assume $f_{\rho}\in\mathcal{H}$.

This is a standard assumption in learning theory and assumes that the target function $f_{\rho}$ defined in Eq. (2.1) is indeed realizable, see also [RR19, RR17, CZ07, SS07].

This is a broad model for the noise in the output $y$, containing uniformly bounded or sub-Gaussian noise; and is in fact weaker than the standard Bernstein condition, e.g., in [BK10].

This is a standard setting in high-dimensional statistics and random matrix theory [EK10, DW18, LR20, HMRT19, EKZ+20] that assumes that the data are drawn from some not-too-heavy-tailed distribution, with possibly (involved) structure between the entries.

To aid our proof, we need some extra results. In [EK10], it has been shown that the kernel matrix $\bm{K}$ in high dimensions can be well approximated by $\widetilde{\bm{K}^{\operatorname{lin}}}$ in spectral norm, i.e., $\|\bm{K}-\widetilde{\bm{K}^{\operatorname{lin}}}\|_{2}\rightarrow 0$ as $n,d\to\infty$

with non-negative parameters $\alpha$, $\beta$, $\gamma$, and the additional matrix $\bm{E}$ given in Table 2, see some typical examples in Appendix A. Here $\gamma$ is the implicit regularization parameter in kernel estimator that depends on the nonlinear function $h$ in the kernel $k$ and the data structure $\bm{\Sigma}_{d}$. According to Eq. (3.1), denote the shortcut $\widetilde{\bm{X}}:=\beta{\bm{X}\bm{X}^{\!\top}}/{d}+\alpha\bm{1}\bm{1}^{\!\top}$, we show in high dimensions that, $\bm{K}$ admits the same eigenvalue decay as $\widetilde{\bm{X}}$ and $\bm{X}\bm{X}^{\!\top}/d$ (see details in Appendix B). Subsequently, we introduce the following quantity function

which is associated with various quantity functions in [AKT19, DW18, LR20, JŞS+20b, NVKM20] and, as we shall see, plays an important role in determining the variance behavior. We will discuss at length $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ based on different data eigenvalue decays in Section 4.

Formally, our main results of KRR in a high-dimensional regime are stated as follows.

(Basic result) Under Assumptions 1-3, let $0<\delta<1/2$, $\theta=\frac{1}{2}-\frac{2}{8+m}$, $d$ large enough, taking the regularization parameter $\lambda:=\bar{c}n^{-\vartheta}$ with $0\leq\vartheta\leq 1/2$, for any given $\varepsilon>0$, it holds with probability at least $1-2\delta-d^{-2}$ with respect to the draw of ${\bm{X}}$ that

with ${\tt{V}}_{1}:=\frac{\sigma^{2}\beta}{d}\mathcal{N}^{n\lambda+\gamma}_{\widetilde{\bm{X}}}$ and the residual term ${\tt{V}}_{2}$

Remark: The first term in Eq. (3.3) is the bound of the bias, which is independent of $d$ and monotonically decreases with $n$. The sum ${\tt{V}}_{1}+{\tt{V}}_{2}$ is the bound of the variance that depends on both $n$ and $d$. Note that ${\tt{V}}_{2}$ monotonically decreases with $n$, and approaches to zero for a large $n$. Therefore, the error bound for ${\tt{V}}_{1}\asymp\frac{1}{d}\mathcal{N}^{n\lambda+\gamma}_{\widetilde{\bm{X}}}$ is the key part of estimates for the variance and will be discussed in in Section 4, where $n\lambda$ corresponds to the explicit regularization and $\gamma$ the implicit regularization. We will demonstrate that ${\tt{V}}_{1}$ can be monotonically decreasing or unimodal under different regularization schemes. Such monotonic bias and unimodal variance can lead to various behaviors of the excess risk, including monotonically decreasing, double descent, and bell-shaped risk curve, as illustrated in Figure 1 of introduction.

2 Refined result

Based on the basic result, if we consider two additional assumptions, i.e., extending Assumption 1 by considering the regularity of $f_{\rho}$ and studying spectral decay of $k$ via complexity of $\mathcal{H}$, we can obtain a refined result.

(Source condition [CZ07]) For some $0<r\leq 1$, there exists $g_{\rho}\in\mathcal{L}_{\rho_{X}}^{2}$ satisfying $\|g_{\rho}\|_{\mathcal{L}_{\rho_{X}}^{2}}\leq R$ such that $f_{\rho}=L_{K}^{r}g_{\rho}$.

(Capacity condition [CZ07]) For any $\lambda>0$, there exist $Q>0$ and $\eta\in$ such that

The notation $\mathcal{N}(\lambda)$ denotes the “effective dimension” and can be regarded as a “measure of size” of the RKHS. This is a natural and widely used assumption in the literature [CZ07, ZDW13, RR17]. Assumption 5 always holds for $\eta=1$ and $Q=\kappa$ where $\kappa:=\max\{1,\sup_{\bm{x}\in X}\sqrt{k(\bm{x},\bm{x})}\}$ as $L_{K}$ is a trace class operator. Its kernel matrix form is $d_{\bm{K}}^{\lambda}:=\operatorname{tr}\left((\bm{K}+\lambda\bm{I}_{n})^{-1}\bm{K}\right)=\sum_{i=1}^{n}\frac{\lambda_{i}(\bm{K})}{\lambda_{i}(\bm{K})+\lambda}$ [AKM+17, LTOS19]. While Assumption 5 can be further refined to obtain a bound that depends on $d$ [PRDVR20], here we focus on the eigenvalue decay of $\bm{K}$, see Section 4 for details.

Based on the above discussion, we obtain a refined result of Theorem 1 as below.

(Refined result) Under Assumptions 2-5, let $0<\delta<1/2$, $\theta=\frac{1}{2}-\frac{2}{8+m}$, and $d$ large enough, taking $\lambda:=\bar{c}n^{-\vartheta}$ with $0\leq\vartheta\leq\frac{1}{1+\eta}$, then for any given $\varepsilon>0$, it holds with probability at least $1-2\delta-d^{-2}$

where ${\tt{V}}_{1}$ and ${\tt{V}}_{2}$ are the same as in Theorem 1.

Remark: Compared to classical learning theory results [FS17] achieving $\mathcal{O}(n^{-\frac{2r+1}{2r+1+\eta}})$ learning rates, the parameter $\eta$ in our results only effects the selection range of $\lambda$, which is nearly independent of the learning rates to some extent. That means, the spectral decay of a kernel function $k$ in high dimensions is almost irrelevant to its kernel type. In fact, the eigenvalue decay of the kernel matrix in our model largely depends on the data, which is in essence different from classical learning theory results. Therefore, our result reflects a certain “universality” on the kernel function in high dimensional problems, which shows consistency to [EK10].

3 Related work

We provide non-asymptotic results that systematically analyze both implicit and explicit regularization schemes within a unified framework.

Implicit regularization in kernel/linear interpolation: Implicit regularization can be induced by minimum norm solutions in linear interpolation [DLM19, KLS20], or the curvature of the kernel function in kernel interpolation [LR20]. Compared to the risk curve in [LR20] that converges to a non-zero constant, the risk curve in our results tends to zero when $n\gg d$. Hence our result demonstrates that, in the double descent case, the minimum of the expected risk in the second descent is lower than the first descent; while the same claim cannot be obtained from [LR20]. Besides, under the basic $f_{\rho}\in\mathcal{H}$ case, our bias bound is based on the eigen-decay (trends) of the kernel matrix $\bm{K}$ and thus can be (almost) independent of $d$, achieving an optimal learning rate $\mathcal{O}(\lambda)$ in a minimax case. This is different from [LR20] that corresponds to the sum of tailed eigenvalues of $\bm{K}$. Specifically, if we directly set $\lambda$ to zero, our result for the bias still holds, which can be bounded by $\|L_{K,\bm{X}}-L_{K}\|_{\mathcal{L}_{\rho_{X}}^{2}}\lesssim\mathcal{O}(1/\sqrt{n})$.

Explicit regularization in kernel/linear regression: We provide non-asymptotic results that refine a series of asymptotic analyses, e.g., the Stieltjes transform approach in [HMRT19, EKZ+20, YYY+20, JŞS+20a] and the statistical mechanic approach in [CBP20]. In fact, by considering the limiting eigenvalue distribution of $\bm{X}\bm{X}^{\!\top}/d$ via its Stieltjes transform $\frac{1}{n}\mathcal{N}^{b}_{\bm{X}\bm{X}^{\!\top}/d}\!\approx\!m(-b)-bm^{\prime}(-b)$, for $m(b)$ the solution to the popular Marc̆enko–Pastur equation [MP67], our error bound recovers [HMRT19, Theorem 5] with $b:=\lambda$ and isotropic features $\bm{\Sigma}_{d}=\bm{I}_{d}$. Finite sample analyses are often based on a finer control of the Stieltjes transform [JŞS+20b] or the effective rank [BLLT20, CL20]. However, the aforementioned results are generally limited to Gaussian [JŞS+20b, NVKM20] and sub-Gaussian data [BLLT20, CL20, CC20], or Gaussian covariates [RMR20]. Here we consider a much broader family of distributions. Besides, under some specific situations, the regularization parameter $\lambda$ in (generalized) linear regression can be negative [KLS20] or optimal tuned [NVKM20, WX20] so as to generalize well. Recent research [GMMM19, LRZ19, BCP20] on kernel regression in $n:=\mathcal{O}(d^{c})$ shows different trends.

4 Proof framework

The proof of our results is fairly technical and lengthy, and we briefly sketch some main ideas of Theorem 2 here. Note that, Theorem 1 is a special case of Theorem 2 by taking $r=1/2$ and $\eta=1$. The modified error decomposition, the error bounds of variance for radial kernels, and estimates for bias are the main elements of novelty in the proof.

we have $f_{\bm{X},\lambda}=(L_{K,\bm{X}}+\lambda I)^{-1}L_{K,\bm{X}}f_{\rho}$. Accordingly, the variance-bias decomposition is stated in the following lemma, with proof deferred to Appendix C.

It is clear that, the variance term does not depend on the target function $f_{\rho}$, and the bias is independent of the residual error $\bm{\epsilon}$. Proof for the bias {\tt{B}}\lesssim n^{-2\vartheta r}\log^{4}\big{(}\frac{2}{\delta}\big{)} can be found in Appendix D. Proof for the variance ${\tt{V}}\lesssim{\tt{V}}_{1}+{\tt{V}}_{2}$ refers to Appendix E.

Discussion on Error Bounds

In this section, we discuss our Theorem 2 for different eigenvalue profiles of $\widetilde{\bm{X}}$ in the two regimes of $n<d$ and $n>d$. Since $\bm{K}$ shares the same eigenvalue decay as $\bm{X}\bm{X}^{\!\top}/d$ and $\widetilde{\bm{X}}$ (see Proposition B.1 in Appendix B), we do not distinguish the eigen-decay of these two data matrices in the subsequent discussions. We first focus on the variance ${\tt{V}}$ that can be unimodal or monotonically decreasing with $n$ under different regularization schemes. Subsequently, we investigate the total risk curve as the sum of bias and variance. Note that $\widetilde{\bm{X}}$ has different numbers of non-zero eigenvalues under the two regimes, we denote $r_{*}:=\operatorname{rank}(\widetilde{\bm{X}})\leq\min\{n,d\}$, which, as we shall see, plays a significant role in characterizing the different cases of our bounds.

We consider here three eigenvalue decays of $\widetilde{\bm{X}}$: harmonic, polynomial, and exponential decay [Bac13, LTOS19].

Under the three eigenvalue decays in Table 3, denote $r_{*}=\operatorname{rank}(\widetilde{\bm{X}})$, then the quantity function $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ with $b:=n\lambda+\gamma$ can be bounded by 1) harmonic decay: $\mathcal{N}^{b}_{\widetilde{\bm{X}}}\leq\frac{n}{b^{2}}\ln\frac{n+(r_{*}+1)b}{n+b}=\mathcal{O}(\frac{n}{b^{2}})$. 2) polynomial decay: $\mathcal{N}^{b}_{\widetilde{\bm{X}}}\leq\frac{\widetilde{C}}{2ab}\left(\frac{n}{b}\right)^{\frac{1}{2a}}$, where $\widetilde{C}$ is some constant. 3) exponential decay: $\mathcal{N}^{b}_{\widetilde{\bm{X}}}\!\leq\!\frac{1}{a}\left(\frac{1}{b+ne^{-a(r_{*}+1)}}\!-\!\frac{1}{b+ne^{-a}}\right)$.

Proof The proof can be found in Appendix F. ∎

According to Proposition 4.1, we summarize our results in Table 1 and discuss them as follows:

Harmonic decay: ${\tt{V}}_{1}\leqslant\mathcal{O}(\frac{n}{b^{2}d})$.

For $\lambda=0$, i.e., the ridgeless case, we have $b=\gamma=\mathcal{O}(1)$, and ${\tt{V}}_{1}\leq\mathcal{O}(\frac{n}{d})$, which indicates ${\tt{V}}_{1}$ increases with $n$ in the $n<d$ regime. For $\lambda\neq 0$, taking $\lambda:=\bar{c}n^{-\vartheta}$, we have ${\tt{V}}_{1}\leqslant\mathcal{O}(\frac{n}{d(\bar{c}n^{1-\vartheta}+\gamma)^{2}})$. To investigate the monotonicity of $g(n):=\frac{n}{d(\bar{c}n^{1-\vartheta}+\gamma)^{2}}$, define $n_{*}:=\left(\frac{\gamma}{2-2\vartheta-\bar{c}}\right)^{\frac{1}{1-\vartheta}}$, we find that, a large $\lambda$ leads to a small $n_{*}$. According to the relationship between $r_{*}$, $n_{*}$, and $d$, we can conclude that (see Table 1 and the red curve in Figure 1(a)):

When $\vartheta\geq\frac{1}{2(2-\bar{c})}$, ${\tt{V}}_{1}$ will increase with $n$ until $n:=r_{*}$ and then remain unchanged when $r_{*}<n<d$. When $\vartheta<\frac{1}{2(2-\bar{c})}$, there are various trends as follows: 1) if $d<n_{*}$, this is the same as the $\vartheta\geq\frac{1}{2(2-\bar{c})}$ case; 2) if $r_{*}<n_{*}<d$, ${\tt{V}}_{1}$ will increase with $n$ until $n:=r_{*}$, and then remain unchanged when $r_{*}<n<d$; 3) if $n_{*}<r_{*}<d$, ${\tt{V}}_{1}$ will increase with $n$ until $n:=n_{*}$ and then decrease with $n$ until $n:=r_{*}$, and stay unchanged on $r_{*}<n<d$; 4) If $n_{*}<c$ such that $n>c$ always holds for some constant $c$, we have ${\tt{V}}_{1}$ increases with $n$ until $n:=r_{*}$, and then stays unchanged on $r_{*}<n<d$. Remark that, for $\gamma<2-2\vartheta-\bar{c}$, we have $n_{*}<1$ and thus $n>n_{*}$, so that ${\tt{V}}_{1}$ always decreases with $n$ until $n:=r_{*}$.

Polynomial decay: ${\tt{V}}_{1}\leqslant\mathcal{O}(\frac{1}{bd}(\frac{n}{b})^{\frac{1}{2a}})$.

Similar to above, define $n_{*}=\left(\frac{\gamma}{2a\bar{c}[1-(1+\frac{1}{2a})\vartheta]}\right)^{\frac{1}{1-\vartheta}}$, we obtain results similar to the case of harmonic decay, but with different thresholds: $\vartheta\geq(1+\frac{1}{2a})^{-1}$ and $\vartheta<(1+\frac{1}{2a})^{-1}$, see Table 1 for details.

Exponential decay: {\tt{V}}_{1}\leq\frac{\widetilde{C}\beta}{ad}\big{(}\frac{1}{b+ne^{-a(r_{*}+1)}}\!-\!\frac{1}{b+ne^{-a}}\big{)}.

Here we consider the monotonicity of the function G(n):=\big{(}\frac{1}{b+ne^{-a(r_{*}+1)}}\!-\!\frac{1}{b+ne^{-a}}\big{)} with $b:=n\lambda+\gamma$ to study the trend of ${\tt{V}}_{1}$ regarding to $n$. Let $n_{*}$ be the solution of the equation $G^{\prime}(n)=0$, then we have the similar conclusion with that of harmonic decay and polynomial decay by the relationship between $n_{*}$, $r_{*}$, and $d$, see Table 1 for details. More specifically, under some certain conditions, ${\tt{V}}_{1}$ is able to monotonically decrease with $n$, refer to Appendix F.1 for details.

2 Variance trends for n>d𝑛𝑑n>d and total risk

Different from the above $n<d$ case, the current $n>d$ regime admits that $\widetilde{\bm{X}}$ has at most $d$ non-zero eigenvalues. In this under-parameterized regime, we are particularly interested in the behavior as $n\to\infty$. In Appendix F.2, we prove that ${\tt{V}}_{1}$ approaches to zero as $n\to\infty$ under the above three eigenvalue decays.

Based on the above discussions in the $n>d$ and $n<d$ regimes, we conclude that, the variance can be unimodal (small regularization) or decreasing (large regularization) as $n$ grows, which, together with the fact that the bias is monotonically decreasing with $n$, leads to the following three configurations for the total risk: (i) if the bias dominates at small $n$ and then decays fast (i.e., with a small regularization), we observe a double descent curve as in Figure 1(b); (ii) if the bias dominates but decays slowly (with a large regularization), the risk curve will be monotonic decreasing as in Figure 1(d); (iii) if the variance dominates, a bell-shaped risk curve as in Figure 1(c) will be observed.

Numerical Results

In this section, experiments are conducted to validate our theoretical resultsThe source code of our implementation can be found in http://www.lfhsgre.org.. Polynomial kernel of degree $3$ and Gaussian kernel are evaluated on 1) a synthetic dataset that satisfies our technical assumptions and 2) a subset of the YearPredictionMSD dataset [Cha08] with 1,000 data samples and $d=90$, to study our derived error bounds for the bias and variance. More experimental results can be found in Appendix G.

Eigenvalue decay equivalence: Here we study the eigenvalue decay of the original polynomial/Gaussian kernel matrices and their linearization $\bm{X}\bm{X}^{\!\top}/d$ on the subset of YearPredictionMSD dataset. Note that, polynomial kernels $k(\bm{x},\bm{x}^{\prime}):=\left(1+\langle\bm{x},\bm{x}^{\prime}\rangle/d\right)^{p}$ admit $\beta:=p$ independent of $\bm{\Sigma}_{d}$ (see in Table 4), so we use the linearization $\beta\bm{X}\bm{X}^{\!\top}/d$ for this kernel. Results in Figure 2 demonstrate that, the original nonlinear kernels admit the same eigenvalue decay as $\bm{X}\bm{X}^{\!\top}/d$. More experimental results on various dataset can be found in Appendix G.1.

Risk curves on synthetic dataset: To quantitatively assess our derived error bounds for the bias and variance, we generate a synthetic dataset under a known $f_{\rho}$, with harmonic decay for the data as an illustrating example. More experimental results on different eigenvalue decays refer to Appendix G.2. To be specific, we assume $y_{i}=f_{\rho}(\bm{x}_{i})+\varepsilon$ with target function $f_{\rho}(\bm{x})=\sin(\|\bm{x}\|^{2}_{2})$ and Gaussian noise $\varepsilon$ having zero-mean and unit-variance. The feature dimension $d$ is set to 500. The samples are generated from $\bm{x}_{i}=\bm{\Sigma}_{d}^{1/2}\bm{t}_{i}$ (and thus $\bm{X}^{\!\top}\bm{X}=\bm{T}^{\!\top}\bm{\Sigma}_{d}\bm{T}$ with $\bm{T}=[\bm{t}_{1},\bm{t}_{2},\cdots,\bm{t}_{n}]^{\!\top}$) by the following steps: (i) take $\bm{\Sigma}_{d}$ as a diagonal matrix with its diagonal entries following with harmonic decay, i.e., $(\bm{\Sigma}_{d})_{ii}\propto n/i$. (ii) take $\bm{T}$ as a random orthogonal matrixWe generate a random Gaussian matrix and use the QR decomposition to obtain an orthogonal matrix [YSC+16]. such that $\bm{T}^{\!\top}\bm{\Sigma}_{d}\bm{T}$ also has a harmonic eigen-decay with $\bm{T}$ having almost i.i.d entries.

Accordingly, the above generation process satisfies Assumption 3, and also $\bm{X}\bm{X}^{\!\top}/d$ admits the same eigenvalue decay as $\bm{\Sigma}_{d}$, which can be used to validate our discussion in Section 4. In this setting, the expected excess risk, the bias, and the variance can be directly computed to validate our derived error bounds. The experimental results are validated across 10 trials. Specifically, to disentangle the implicit regularization effect of KRR on the final result, we apply the linearization of the polynomial/Gaussian kernel by setting $\gamma=0$ in Eq. (3.1). In this case, the explicit $\lambda:=\bar{c}n^{-\vartheta}$ is the only regularization in KRR. In our model, $\bar{c}$ is empirically set to $0.01$ to avoid a large $\lambda$ when $n$ is small.

Figures 3 and 4 show results under the harmonic decay setting for the linearization of the polynomial/Gaussian kernel, respectively. We observe that: 1) our error bound ${\tt{V}}_{1}\asymp\frac{1}{d}\mathcal{N}^{n\lambda}_{\widetilde{\bm{X}}}$ exhibits the same trend as the true variance; 2) in this case, the variance dominates and we thus obtain a bell-shaped risk curve that first increases and then decreases; 3) as $\vartheta$ decreases, $\lambda$ increases and the peak point of the variance occurs at smaller and smaller $n$; 4) the bias monotonically decreases with $n$, which corresponds to our error bound for the bias at a certain $\mathcal{O}(n^{-2\vartheta r})$ rate in Theorem 2 by taking $r=1$ as the used $f_{\rho}$ is smooth enough to achieve a good approximation error; 5) in our high-dimensional regimes, different kernels lead to the same convergence rates of the bias, which verifies our results but is different from those in classical learning theory.

Risk curves on the real-world datasets: Figure 5(a) shows the relative mean squared error (RMSE) of kernel ridgeless regression and its linearization in Eq. (3.1) on a subset (1,000 examples) of the YearPredictionMSD dataset averaged on 10 trials. Figure 5(b) shows the classification accuracy of such two methods on the MNIST dataset [LBBH98]. To evaluate the effectiveness of our error bounds, we plot the re-scaled ${\tt{V}}_{1}\asymp\frac{1}{d}\mathcal{N}^{\gamma}_{\widetilde{\bm{X}}}$ with $\lambda=0$. It can be found that, kernel interpolation estimator generalizes well due to the implicit regularization, i.e., $\gamma\neq 0$, which also exhibits a bell-shaped risk curve as our theoretical results suggest. However, in Figure 5(b), the risk curve monotonically decreases with $n$ on the MNIST dataset [LBBH98], and at the same time kernel interpolation estimator and its linearization appear to generalize well. This observation may due to the implicit regularization parameter $\gamma$ in Eq. (3.1) (of $10^{-3}$ order on this dataset) that plays a fundamental role of “self-regularization”. Accordingly, the proposed analysis provides access to the high-dimensional classification problem that may establish more involved behavior than double descent, despite a clear mismatch between real-world data and the technical Assumption 3, thereby conveying a strong practical motivation for the present analysis.

Conclusion

We derived non-asymptotic expressions for the expected excess risk of kernel ridge regression estimators in the under- and over-determined regimes. The used linearization technique of nonlinear smooth kernel allows us to discuss the impact of implicit and explicit regularization in a systematic manner. Our refined analysis demonstrates that the monotonic bias and unimodal variance are able to exhibit various trends of risk curves. Since it is enough to require that the kernel function is differentiable in a neighborhood, our results further extend to the case of Laplace kernels [RZ19].

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation program / ERC Advanced Grant E-DUALITY (787960). This paper reflects only the authors’ views and the Union is not liable for any use that may be made of the contained information. This work was supported in part by Research Council KU Leuven: Optimization frameworks for deep kernel machines C14/18/068; Flemish Government: FWO projects: GOA4917N (Deep Restricted Kernel Machines: Methods and Foundations), PhD/Postdoc grant. This research received funding from the Flemish Government (AI Research Program). This work was supported in part by Ford KU Leuven Research Alliance Project KUL0076 (Stability analysis and performance improvement of deep reinforcement learning algorithms), EU H2020 ICT-48 Network TAILOR (Foundations of Trustworthy AI - Integrating Reasoning, Learning and Optimization), Leuven.AI Institute.

References

Appendix A Examples of kernels and their linearizations

In this section, we present linearization of some typical kernels by Eq. (3.1). Here we assume that $\alpha,\beta,\gamma\geq 0$ to ensure the positive definiteness of the approximated kernel matrix $\widetilde{\bm{K}}^{\text{lin}}$. Table 4 reports the results of three inner-product kernels including polynomial kernel, linear kernel, exponential kernel; as well as a radial kernel: the common-used Gaussian kernel. We can find that $\alpha,\gamma\geq 0$. Specifically, $\beta>0$ avoids a trivial solution.

Appendix B Eigenvalue decay equivalence

In this section, we demonstrate that, in high dimensions, a kernel matrix induced by inner-product kernels or radial kernels admits the same eigenvalue decay as $\widetilde{\bm{X}}=\beta{\bm{X}\bm{X}^{\!\top}}/{d}+\alpha\bm{1}\bm{1}^{\!\top}$ and ${\bm{X}\bm{X}^{\!\top}}/{d}$.

For notational simplicity, denote the inner-product kernel matrix $\bm{K}_{\operatorname{inner}}$ and its linearization $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{inner}}}$; the radial kernel matrix $\bm{K}_{\operatorname{radial}}$ and its linearization $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{radial}}}$.

The inner-product kernel matrix $\bm{K}_{\operatorname{inner}}$ admits the same eigenvalue decay as $\widetilde{\bm{X}}$ and ${\bm{X}\bm{X}^{\!\top}}/{d}$.

Proof According to Theorem 2.1 in [EK10], the inner-product kernel matrix $\bm{K}_{\operatorname{inner}}$ can be well approximated by $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{inner}}}$ with

in a spectral norm sense, where $\alpha$, $\beta$, $\gamma$ are given in Table 2. As a result, with high probability, the inner-product kernel matrix $\bm{K}_{\operatorname{inner}}$ and its linearization $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{inner}}}$ has the same eigenvalue. That means, $\bm{K}_{\operatorname{inner}}$ admits the same eigenvalue decay as $\widetilde{\bm{X}}:=\beta{\bm{X}\bm{X}^{\!\top}}/{d}+\alpha\bm{1}\bm{1}^{\!\top}$ via a constant shift $\gamma$.

Next, we shall demonstrate that $\bm{K}_{\operatorname{inner}}$ admits the same eigenvalue decay as $\bm{X}\bm{X}^{\!\top}/d$. Since $\bm{1}\bm{1}^{\!\top}$ is a rank-one matrix with $\lambda_{1}(\bm{1}\bm{1}^{\!\top})=n$, with Weyl’s inequality and $\lambda_{n}\leq\lambda_{n-1}\leq\ldots\leq\lambda_{1}$, we have

so that the eigenvalue of $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{inner}}}$ interlaced with those of $\beta{\bm{X}\bm{X}^{\!\top}}/{d}+\gamma\bm{I}$. We can thus conclude that the eigenvalue decay of $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{inner}}}$ is the same as that of ${\bm{X}\bm{X}^{\!\top}}/{d}$ with a constant shift and scaling, which do not effect the trend of eigenvalue decay. Accordingly, the inner-product-type kernel matrix ${\bm{K}}_{\operatorname{inner}}$ and its linearization $\widetilde{{\bm{K}}^{\operatorname{lin}}_{\operatorname{inner}}}$, $\widetilde{\bm{X}}$ admit the same eigenvalue decay as ${\bm{X}\bm{X}^{\!\top}}/{d}$, which concludes the proof. ∎ Proposition B.1 also provides a justification to study the eigenvalue decay of a radial kernel matrix. According to Theorem 2.2 in [EK10], the radial kernel matrix $\bm{K}_{\operatorname{radial}}$ can be well approximated by $\widetilde{{\bm{K}}_{\operatorname{radial}}^{\operatorname{lin}}}$ with

Accordingly, ${\bm{K}}_{\operatorname{radial}}$ admits the same eigenvalue decay as ${\bm{X}\bm{X}^{\!\top}}/{d}$.

Appendix C Proof of Lemma 3.1

Proof By virtue of the closed form of the KRR estimator in Eq. (2.2) and $\bm{\epsilon}:=\bm{y}-f_{\rho}(\bm{X})$, we have

Based on the definition of ${\tt{B}}$, we decompose ${\tt{B}}$ as

Appendix D Proof for the bias

The error bound for the bias is given by the following theorem.

(Bias) Under Assumption 4 (source condition with $0<r\leq 1$), Assumption 5 (capacity condition with $0\leq\eta\leq 1$), let $0<\delta<1/2$, taking the regularization parameter $\lambda:=\bar{c}n^{-\vartheta}$ with $0\leq\vartheta\leq\frac{1}{1+\eta}$, there holds with probability at least $1-2\delta$, we have

In our error decomposition, $\|f_{\lambda}-f_{\rho}\|^{2}_{\mathcal{L}^{2}_{\rho_{X}}}$ is independent of data ${\bm{X}}$ that corresponds to the approximation error in learning theory [CZ07]; while the first term $\|f_{\bm{X},\lambda}-f_{\lambda}\|^{2}_{\mathcal{L}^{2}_{\rho_{X}}}$ depends on ${\bm{X}}$, termed as bias-sample error. To prove Theorem 3, we need to bound the approximation error and the bias-sample error as follows.

In learning theory, the approximation error $\|f_{\lambda}-f_{\rho}\|_{\mathcal{L}^{2}_{\rho_{X}}}$ can be estimated by the source condition in Assumption 4.

(Lemma 3 in [SZ07]) Under the source condition in Assumption 4 with $0<r\leq 1$, the approximation error can be given by

D.2 Bound bias-sample error

To bound the bias-sample error $\|f_{\bm{X},\lambda}-f_{\lambda}\|_{\mathcal{L}^{2}_{\rho_{X}}}$, we need the following lemma.

(Lemma 17 in [LGZ17]) For any $0<\delta<1$, it holds with probability at least $1-\delta$ that

where $\kappa:=\max\{1,\sup_{\bm{x}\in X}\sqrt{k(\bm{x},\bm{x})}\}$.

Then the bias-sample error can be decomposed into several parts.

Proof [Proof of Lemma D.3] According to the definition of $f_{\bm{X},\lambda}$ and $f_{\lambda}$, we have

Due to $(A+\lambda I)^{-1}A=I-\lambda(A+\lambda I)^{-1}$ for any bounded positive operator $A$, we have

Further, by virtue of the first order decomposition of operator difference: $A^{-1}-B^{-1}=A^{-1}(B-A)B^{-1}$ for any invertible bounded operator and using the source condition in Assumption 4, the above equation can be further expressed as

Besides, using $\|AB^{t}\|\leq\|A\|^{1-t}\|AB\|^{t}$ with $t\in$ for any bounded linear operator $A$ and positive semi-definite operator $B$ in Proposition 9 in [RR17], we have

where we choose $A:=(L_{K}+\lambda I)^{-1/2}(L_{K}-L_{K,\bm{X}})$, $B:=(L_{K}+\lambda I)^{-1}$, and $t:=1-r\in[0,1)$. Accordingly, we can conclude our proof due to $\|(L_{K,\bm{X}}+\lambda I)^{-1/2}\|\leq 1/\sqrt{\lambda}$ and $\|(L_{K}+\lambda I)^{-r}L_{K}^{r}\|\leq 1$. ∎ Remark: The proof framework of Lemma D.3 is similar to Lemma 4 in [RR17] but we consider a more general case $0<r\leq 1$ than $1/2\leq r\leq 1$ in [RR17]. Although $0<r<1/2$ appears to be unattainable as claimed in [RR17], we follow with [LGZ17, GSW17] on a quite general case with $r>0$.

To prove Theorem 3, we also need the following two lemmas.

(Proposition 6 in [RR17]) Let $\delta\in(0,1/2]$, it holds with probability at least $1-2\delta$ that

For any $0<\delta<1$, with probability at least $1-\delta$, we have

Proof [Proof of Lemma D.5] By virtue of a second order decomposition of operator difference in Lemma 16 [LGZ17], we have

Accordingly, denote $A:=L_{K,\bm{X}}+\lambda I$ and $B:=L_{K}+\lambda I$, we can derive that

where $\mathcal{A}:=\frac{2\kappa}{\sqrt{n\lambda}}\left\{\frac{\kappa}{\sqrt{n\lambda}}+\sqrt{\mathcal{N}(\lambda)}\right\}\log(2/\delta)$ by Lemma D.2. The first inequality holds by $\|A^{s}B^{s}\|\leq\|AB\|^{s}$ with $0\leq s\leq 1$ for positive operators $A$ and $B$ on Hilbert spaces [BK10]. The second inequality can be derived by Eq. (D.1), $\|(L_{K,\bm{X}}+\lambda I)^{-1}\|\leq 1/\lambda$ and $\|(L_{K}+\lambda I)^{-1/2}\|\leq 1/\sqrt{\lambda}$. ∎ Remark: Lemma 7.2 in [RCR13] gives $\|(L_{K,\bm{X}}+\lambda I)^{-1/2}(L_{K}+\lambda I)^{1/2}\|\leq\sqrt{2}$ by assuming $\lambda>\frac{9}{n}\log\frac{n}{\delta}$; whereas our result does not require extra conditions on $\lambda$.

Based on the above lemmas, we are ready to prove Theorem 3.

Proof [Proof of Theorem 3] We first estimate $\|(L_{K,\bm{X}}+\lambda I)^{-1/2}(L_{K}+\lambda I)^{1/2}\|$ in Lemma D.5 by taking $\lambda:=\bar{c}n^{-\vartheta}$ and the capacity condition in Assumption 5: $\mathcal{N}(\lambda)\leq Q^{2}\lambda^{-\eta}$ with $\eta\in$. Accordingly, we have

where we use $\log^{r}(2/\delta)\leq\log(2/\delta)$ due to $\log(2/\delta)>1$ in the last inequality. Since $\|(L_{K,\bm{X}}+\lambda I)^{-1/2}(L_{K}+\lambda I)^{1/2}\|$ converges to zero when $n$ is large enough, we require $\vartheta<\frac{1}{1+\eta}$ to ensure a positive convergence rate, which implies $\vartheta\leq 1$. Then we bound $\|(L_{K}+\lambda I)^{-1/2}(L_{K}-L_{K,\bm{X}})\|^{r}$ by Lemma D.2. By virtue of $(a+b)^{r}\leq a^{r}+b^{r}$ for any $r\in(0,1]$ and $a,b\geq 0$, we have

where the second one admits by the capacity condition in Assumption 5. Similarly, to bound $\|(L_{K}+\lambda I)^{-1/2}(L_{K}-L_{K,\bm{X}})(L_{K}+\lambda I)^{-1}\|^{1-r}$ by Lemma D.4, we can derive that

Combining the above three inequalities, we have

where $\widetilde{C_{R,Q,\kappa,\bar{c}}}:=4R\kappa^{3}(Q+\kappa)^{2}/\bar{c}^{3}$ is independent of $n$ and $d$.

where the third inequality holds by $2\vartheta r\leq 1-\vartheta(\eta r+1-2r)$ due to $\vartheta\leq\frac{1}{1+\eta}$, and $\widetilde{C},\widetilde{C_{1}}$ are some constants independent of $n$ and $d$. Accordingly, we can conclude the proof. ∎

Appendix E Proof for the variance

Formally, we have the following theorem to bound the variance.

(Variance) Under Assumptions 2, 3, then for $0<\delta<1$ with probability $1-\delta-d^{-2}$, $\theta=\frac{1}{2}-\frac{2}{8+m}$, and $d$ large enough, for any given $\varepsilon>0$, we have

where ${\tt{V}}_{1}:=\frac{\sigma^{2}\beta}{d}\mathcal{N}^{n\lambda+\gamma}_{\widetilde{\bm{X}}}$ and ${\tt{V}}_{2}$ is the residual term with

For inner-product kernels, our proof framework follows [LR20], and is briefly discussed in Section E.1. Nevertheless, error bound on radial kernels has not been investigated in [LR20] and is more subtle to handle (than that of inner-product kernels) due to the additionally introduced $\bm{A}$ and $\bm{A}\odot\bm{A}$ in Table 2. Accordingly, we mainly focus on proofs for radial kernels.

In this subsection, we consider the inner-product kernel case with $k(\bm{x},\bm{x}^{\prime})=h\left(\langle\bm{x},\bm{x}^{\prime}\rangle/d\right)$. We briefly introduce our results that can be derived from proofs of Theorem 2 in [LR20] for completeness.

and $k^{\operatorname{lin}}(\bm{X},\bm{x})$ is the transpose of $k^{\operatorname{lin}}(\bm{x},\bm{X})$. Note that $\gamma$ in $\widetilde{\bm{K}^{\operatorname{lin}}}$ corresponds to the implicit regularization and $n\lambda$ corresponds to the explicit regularization. Now we prove Theorem 4 for inner-product kernels. Proof [Proof of Theorem 4 for inner-product kernels] According to the definition of V, we have

where the first inequality comes from Assumption 2. To bound the terms in Eq. (E.2), we need

Combining the above results, with probability at least $1-\delta-d^{-2}$, for any given $\varepsilon>0$, The error bound for the variance in Eq. (E.2) can be further given by

E.2 Radial kernel matrices

In this subsection, we consider the radial kernel case with $k(\bm{x},\bm{x}^{\prime})=h\left(\frac{1}{d}\|\bm{x}-\bm{x}^{\prime}\|_{2}^{2}\right)$. Since the linearization of radial kernel matrices incurs in two additionally terms $\bm{A}$ and $\bm{A}\odot\bm{A}$, estimation for radial kernels is more technical than that of inner-product kernels. Accordingly, to prove Theorem 4 for radial kernels, we need to introduce the following notations and auxiliary results.

Recall $\tau:=\operatorname{tr}(\bm{\Sigma}_{d})/d$, define

where $\bm{A}(\bm{x},\bm{X}):=\psi_{\bm{x}}+[\psi_{1},\psi_{2},\cdots,\psi_{n}]^{\!\top}$ with $\psi_{\bm{x}}=\|\bm{x}\|^{2}_{2}/d-\tau$ and $\psi_{i}=\|\bm{x}_{i}\|^{2}_{2}/d-\tau$ for $i=1,2,\ldots,n$. As discussed in Appendix B, we conclude that $\widetilde{\bm{K}^{\operatorname{lin}}}$ admits the same eigenvalue decay as $\widetilde{\bm{X}}$ since $\bm{A}$ is a rank-2 matrix. Accordingly, we have the following results.

Note that, the matrix $\mathop{\rm diag}(\bm{\Sigma}_{d})\bm{1}_{n}^{\!\top}$ is a rank-one matrix, which implies $\operatorname{rank}(\bm{X}\bm{\Sigma}_{d}^{1/2}\mathop{\rm diag}(\bm{\Sigma}_{d})\bm{1}_{n}^{\!\top})\leq 1$. Accordingly, its non-zero eigenvalue $\lambda_{1}(\bm{X}\bm{\Sigma}_{d}^{1/2}\mathop{\rm diag}(\bm{\Sigma}_{d})\bm{1}_{n}^{\!\top})$ admits

Proof [Proof of Proposition E.2] By virtue of the following results [EK10]

Given a radial kernel, under Assumption 3, for $\theta=\frac{1}{2}-\frac{2}{8+m}$, we have with probability at least $1-d^{-2}$ with respect to the draw of $\bm{X}$, for $d$ large enough, for any given $\varepsilon>0$, we have

where $\widetilde{C_{1}}$ is some constant independent of $n$ and $d$.

Remark: In fact, we only need the $(5+m)$-moment in Assumption 3 but we still follow with it for simplicity.

Proof [Proof of Lemma E.3] We start with the entry-wise Taylor expansion for the smooth kernel at $2\tau$ with $\tau:=\operatorname{tr}(\bm{\Sigma}_{d})/d$

where $\psi_{j}=\|\bm{x}_{j}\|^{2}_{2}/d-\tau$ for $j=1,2,\dots,n$ as defined before. Accordingly, by virtue of $k^{\operatorname{lin}}(\bm{x},\bm{x}_{j})=\frac{\beta\bm{x}^{\!\top}\bm{x}_{j}}{d}-\frac{\beta}{2}(\psi_{\bm{x}}+\psi_{j})$ and Corollary 2 in [EK10], with probability at least $1-d^{-2}$, for any $\varepsilon>0$, we have

where we only need $(5+m)$-moment. Therefore, with probability at least $1-d^{-2}$, for any given $\varepsilon>0$, we have

where $\widetilde{C_{1}}$ and $\widetilde{C_{2}}$ are some constant independent of $n$ and $d$. ∎

E.2.2 Proofs of Theorem 4 for radial kernels

In [EK10], the approximation error between radial kernel matrices and their linearization can be decomposed into three parts: the first-order term $A_{1}$, the second-order term $A_{2}$, and the third-order term $A_{3}$

where $A_{1}$ and $A_{3}$ admit $\|A_{1}\|_{2}\leq d^{-\theta}\log^{2+4\varepsilon}d$ and $\|A_{3}\|_{2}\leq d^{-\theta}\log^{2+4\varepsilon}d$. The second-order term $A_{2}$ admits $\operatorname{Pr}(\|A_{2}\|\leq d^{-\theta}\delta^{-1/2})\leq\delta$ by Proposition A.2 in [LR20] and [EK10]. Accordingly, with probability at least $1-\delta-d^{-2}$, for $\theta=\frac{1}{2}-\frac{2}{8+m}$ and any given $\varepsilon>0$, we have

According to Proposition E.1 and E.2, we have

It can be found that, the above error bounds are the same as that of inner-product kernels, except two additional terms due to the considered $\bm{A}$ and $\bm{A}\odot\bm{A}$ in the linearization, which can be shown small in the large $n,d$ regime.

By virtue of $\left\|(\bm{K}+n\lambda\bm{I})^{-1}\right\|_{2}\leq\frac{2}{n\lambda+\gamma}$ and $\left\|(\bm{K}+n\lambda\bm{I})^{-1}\widetilde{\bm{K}^{\operatorname{lin}}}\right\|_{2}\leq 2$ in [LR20], Lemma E.3, and the above equations, with probability at least $1-\delta-d^{-2}$, for any given $\varepsilon>0$, we have

where the second inequality admits by Lemma E.3, and the last inequality follows by Eq. (E.5). Finally, we conclude the proof. ∎

Appendix F Proof of Proposition 4.1

In this section, we discuss $\mathcal{N}^{n\lambda+\gamma}_{\widetilde{\bm{X}}}$ based on three eigenvalue decays: harmonic decay, polynomial decay, and exponential decay under two regimes $n<d$ and $n>d$.

Recall $b:=n\lambda+\gamma>0$, and $\mathcal{N}^{b}_{\widetilde{\bm{X}}}:=\sum_{i=1}^{n}\frac{\lambda_{i}(\widetilde{\bm{X}})}{\left[b+\lambda_{i}(\widetilde{\bm{X}})\right]^{2}}$, define $F(\lambda_{i}):=\frac{\lambda_{i}}{(b+\lambda_{i})^{2}}$ where $\lambda_{i}$ is short for $\lambda_{i}(\widetilde{\bm{X}})$. We notice that, when $\lambda_{i}\leq b$, $F(\lambda_{i})$ is an increasing function of $\lambda_{i}$, and thus a decreasing function of $i$ when the above three eigenvalue decays are considered. Likewise, when $\lambda_{i}\geq b$, $F(\lambda_{i})$ is a decreasing function of $\lambda_{i}$, and thus an increasing function of $i$. Without loss of generality, we assume that the first $q$ eigenvalues satisfy $\lambda_{i}\geq b$ with $i=1,2,\cdots,q$ and the remaining $n-q$ eigenvalues satisfy $\lambda_{i}\leq b$ with $i=m+1,m+2\cdots,n$. Clearly, the integer $q$ can be chosen from to $n$. Accordingly, denote $r_{*}:=\operatorname{rank}(\widetilde{\bm{X}})$ which includes the rank-deficient case, $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ can be upper bounded by the Riemann sum as follows.

Harmonic decay $\lambda_{i}(\widetilde{\bm{X}})\propto n/i$ for $i\in\{1,2,\dots,r_{*}\}$ and $\lambda_{i}(\widetilde{\bm{X}})=0$ for $i\in\{r_{*}+1,\dots,n\}$

Polynomial decay: $\lambda_{i}(\widetilde{\bm{X}})\propto ni^{-2a}$ with $a>1/2$ for $i\in\{1,2,\dots,r_{*}\}$ and $\lambda_{i}(\widetilde{\bm{X}})=0$ for $i\in\{r_{*}+1,\dots,n\}$. Hence, we actually aim to bound

Exponential decay: $\lambda_{i}(\widetilde{\bm{X}})\propto ne^{-ai}$ with $a>0$ for $i\in\{1,2,\dots,r_{*}\}$ and $\lambda_{i}(\widetilde{\bm{X}})=0$ for $i\in\{r_{*}+1,\dots,n\}$.

Note that, the monotonicity of $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ (also ${\tt{V}}_{1}$) with respect to $n$ is relatively clear for harmonic decay and polynomial decay but is unclear in the case of exponential decay. Here we study the monotonicity in the exponential decay. Denote the function G(n):=\big{(}\frac{1}{b+ne^{-a(r_{*}+1)}}\!-\!\frac{1}{b+ne^{-a}}\big{)} with $b:=n\lambda+\gamma$, taking $\lambda:=\bar{c}n^{-\vartheta}$, its derivation is

It can be found that both $H_{1}(n)$ and $H_{2}(n)$ are decreasing functions with $n$. More specifically, their maximum and minimum can be achieved with

Accordingly, if $H_{1}(1)<H_{2}(1)$, we obtain a decreasing function $G(n)$ of $n$, which implies that $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ will decrease with $n$. Here the condition $H_{1}(1)<H_{2}(1)$ indicates

Accordingly, if the above inequality holds, $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ will decrease with $n$. In Section G.2, we will experimentally check whether this condition holds or not.

F.2 n>d𝑛𝑑n>d case and the large n𝑛n limit

In this section, we consider the $n>d$ case, and further study the trend of ${\tt{V}}_{1}$ as $n\rightarrow\infty$. Note that, in this case, ${\bm{X}\bm{X}^{\!\top}}/{d}$ has at most $r_{*}\leq d$ non-zero eigenvalues. Accordingly, the Riemann sum is counted to $r_{*}$ instead of $n$. Similar to the above description, we also consider the following three eigenvalue decays.

Harmonic decay $\lambda_{i}(\widetilde{\bm{X}})\propto n/i$, $i\in\{1,2,\cdots,d\}$

In particular, taking the limit of $n\rightarrow\infty$, we have

Accordingly, by the squeeze theorem, we can conclude, given $d$, $\mathcal{N}^{b}_{\widetilde{\bm{X}}}$ tends to zero when $n\rightarrow\infty$.

Polynomial decay: $\lambda_{i}(\widetilde{\bm{X}})\propto ni^{-2a}$ with $a>1/2$, $i\in\{1,2,\cdots,d\}$

Exponential decay: $\lambda_{i}(\widetilde{\bm{X}})\propto ne^{-ai}$ with $a>0$, $i\in\{1,2,\cdots,d\}$

Taking the limit of $n\rightarrow\infty$, we can directly have $\lim\limits_{n\rightarrow\infty}\frac{1}{d}\mathcal{N}^{b}_{\widetilde{\bm{X}}}=0$.

Appendix G Additional Experiments

In this section, we present additional experiments including the following parts:

In Section G.1, we add the MNIST dataset [LBBH98] to verify the eigenvalue decay equivalence, and evaluate the effect by different orders in polynomial kernel.

In Section G.2, our model works in a polynomial kernel setting under the polynomial decay and exponential decay of $\widetilde{\bm{X}}$ on the synthetic dataset.

Apart from the YearPredictionMSD dataset in the main text, we add the MNIST dataset [LBBH98] to verify the eigenvalue decay equivalence. We also compute eigenvalues of $\widetilde{\bm{X}}:=\beta\bm{X}\bm{X}^{\!\top}/d+\alpha\bm{1}\bm{1}^{\!\top}$ for validation. Here the parameters $\alpha$ depends on the covariate $\bm{\Sigma}_{d}$, which can be empirically estimated by the sample covariance $\frac{1}{n}\sum_{i=1}^{n}(\bm{x}_{i}-\frac{1}{n}\sum_{j=1}^{n}\bm{x}_{j})(\bm{x}_{i}-\frac{1}{n}\sum_{j=1}^{n}\bm{x}_{j})^{\!\top}$.

Results on the polynomial kernel with order 3 and the Gaussian kernel are presented in Figure 6 and 7, respectively. It can be observed that, the nonlinear kernel matrix $\bm{K}$ admits almost the same eigenvalue as $\widetilde{\bm{X}}:=\beta\bm{X}\bm{X}^{\!\top}/d+\alpha\bm{1}\bm{1}^{\!\top}$ with a constant shift $\gamma$, and accordingly exhibits the same eigenvalue decay with $\widetilde{\bm{X}}$ and $\bm{X}\bm{X}^{\!\top}/d$.