How rotational invariance of common kernels prevents generalization in high dimensions

Konstantin Donhauser, Mingqi Wu, Fanny Yang

Introduction

Traditional analysis establishes good generalization properties of kernel ridge regression when the dimension dd is relatively small compared to the number of samples nn. These minimax optimal and consistency results however become less powerful for modern data sets with large dd close to nn. High-dimensional asymptotic theory aims to fill this gap by providing bounds that assume d,n→∞d,n\to\infty and are often much more predictive of practical observations even for finite dd.

While recent work establishes explicit asymptotic upper bounds for the bias and variance for high-dimensional linear regression, the results for kernel regression are less conclusive in the regime d/nβ→cd/n^{\beta}\to c with β∈(0,1)\beta\in(0,1). In particular, even though several papers show that the variance decreases with the dimensionality of the data, the bounds on the bias are somewhat inconclusive. On the one hand, Liang et al. prove asymptotic consistency for ground truth functions with asymptotically bounded Hilbert norms for neural tangent kernels (NTK) and inner product (IP) kernels. In contrast, Ghorbani et al. show that for uniform distributions on the product of two spheres, consistency cannot be achieved unless the ground truth is a low-degree polynomial. This polynomial approximation barrier can also be observed for random feature and neural tangent regression .

Notably, the two seemingly contradictory consistency results hold for different distributional settings and are based on vastly different proof techniques. While proves consistency for general input distributions including isotropic Gaussians, the lower bounds in the papers are limited to data that is uniformly sampled from the product of two spheres. Hence, it is a natural question to ask whether the polynomial approximation barrier is a more general phenomenon or restricted to the explicit settings studied in . Concretely, this paper addresses the following question

Can we overcome the polynomial approximation barrier when considering different high-dimensional input distributions, eigenvalue decay rates or scalings of the kernel function?

We unify previous distributional assumptions in one proof framework and thereby characterize how the rotational invariance property of common kernels induces a bias towards low-degree polynomials. Specifically, we show that the polynomial approximation barrier persists for

a broad range of common rotationally invariant kernels such as radial basis functions (RBF) with vastly different eigenvalue decay rates, inner product kernels and NTK of any depth .

general input distributions including anisotropic Gaussians where the degree of the polynomial depends only on the growth of deff:=tr⁡(Σd)/∣∣Σd∣∣opd_{\text{eff}}:=\operatorname{tr}(\Sigma_{d})/\left|\left|\Sigma_{d}\right|\right|_{\textrm{op}} and not on the specific structure of Σd\Sigma_{d}. In particular, we cover the distributions studied in previous related works .

different scalings τ\tau with kernel function kτ(x,x′)=k(xτ,x′τ)k_{\tau}(x,x^{\prime})=k(\frac{x}{\sqrt{\tau}},\frac{x^{\prime}}{\sqrt{\tau}}) beyond the classical choice of τ≍deff\tau\asymp d_{\text{eff}}.

As a result, this paper demonstrates that the polynomial approximation barrier is a general high-dimensional phenomenon for rotationally invariant kernels, restricting the set of functions for which consistency can at most be reached to low degree polynomials.

Rotationally invariant kernels are a natural choice if no prior information on the structure of the ground truth is available, as they treat all dimensions equally. Since our analysis covers a broad range of distributions, eigenvalue decays and different scalings, our results motivate future work to focus on the symmetries respectively asymmetries of the kernel incorporating prior knowledge on the structure of the high-dimensional problem, e.g. .

This paper is organized as follows. First of all, we show in Section 2.2 that the bounded norm assumption that previous consistency results rely on is violated as d→∞d\to\infty even for simple functions such as f⋆(x)=e1⊤xf^{\star}(x)=e_{1}^{\top}x. We then introduce our generalized setting in Section 2 and present our main results in Section 3 where we show a lower bound on the bias that increases with the dimensionality of the data. Finally, in Section 4 we empirically illustrate how the bias dominates the risk in high dimensions and therefore limits the performance of kernel regression. As a result, we argue that it is crucial to incorporate prior knowledge of the ground truth function (such as sparsity) in high-dimensional kernel learning even in the noiseless setting and empirically verify this on real-world data.

Problem setting

In this section, we briefly introduce kernel regression estimators in reproducing kernel Hilbert spaces and subsequently introduce our assumptions on the kernel, data distribution and high-dimensional regime.

with λ>0\lambda>0 and the minimum norm interpolator (also called the kernel \sayridgeless estimate)

that can be obtained as the limit of the ridge estimate f^0=lim⁡λ→0f^λ\hat{f}_{0}=\lim_{\lambda\to 0}\hat{f}_{\lambda} for fixed nn. It is well-known that the ridge estimator can attain consistency as n→∞n\to\infty for some sequence of λ\lambda such that λn→0\frac{\lambda}{n}\to 0. Recently, some works have also analyzed the consistency behavior of ridgeless estimates motivated by the curiously good generalization properties of neural networks with zero training error.

Note that consistency in terms of R(f^λ)→0\textbf{R}(\hat{f}_{\lambda})\to 0 as n→∞n\rightarrow\infty can only be reached if the bias vanishes. In this paper, we lower bound the bias B which, in turn, implies a lower bound on the risk and the inconsistency of the estimator. The theoretical results in Section 3 hold for both ridge regression and minimum norm interpolation. However, it is well known that the ridge penalty controls the bias-variance trade-off, and hence we are primarily interested in lower bounds of the bias for the minimum norm interpolant.

2 Prior work on the consistency of kernel regression

For ridge regression estimates in RKHS, a rich body of work shows consistency and rate optimality when appropriately choosing the ridge parameter λ\lambda both in the non-asymptotic setting, e.g. , and the classical asymptotic setting, e.g. , as n→∞n\to\infty.

In words, for simple sequences of sparse product functions, the Hilbert norm diverges as the dimension d→∞d\to\infty. The precise conditions on the kernel and sequence Hd\mathcal{H}_{d} of induced Hilbert spaces can be found in Appendix B. Figure 1 illustrates this phenomenon for f⋆(x)=e1⊤xf^{\star}(x)=e_{1}^{\top}x for the Laplace and exponential inner product kernel.

The discussion so far implies that generalization upper bounds that rely on the bounded Hilbert norm assumption become void even for simple ground truth functions. A natural follow-up question is hence: Do we actually fail to learn sparsely parameterized functions consistently or is it simply a loose upper bound? A recent line of work by Ghorbani et al. shows that kernel regression estimates can indeed only consistently learn polynomials of degree at most log⁡nlog⁡deff\frac{\log n}{\log d_{\text{eff}}} as d,n→∞d,n\to\infty (which we refer to as the polynomial approximation barrier). While the results provide some intuition for the behavior of kernel regression, the proofs heavily rely on significant simplifications that hold for the specific distributional assumptions on the sphere. It is a priori unclear whether they apply to more general settings including the ones considered in . In the next sections, relying on a different proof technique we show that the polynomial approximation barrier indeed holds for a broad spectrum of data distributions that also capture the distributions studied in the papers and for an entire range of eigenvalue decay rates of the kernel functions (e.g. polynomial and exponential decay rates) and choices of the scaling τ\tau. As a consequence, our results suggest that the polynomial approximation barrier is strongly tied to the rotational invariance of the kernel function and not specific to the settings studied so far.

3 Our problem setting

The framework we study in this paper covers random vectors that are generated from a covariance matrix model, i.e. X=Σd1/2WX=\Sigma_{d}^{1/2}W with vector WW consisting of i.i.d entries and their projections onto the d−1d-1-dimensional unit sphere. We now specify all relevant assumptions on the kernel and data distribution.

Throughout this paper, we focus on continuous and rotationally invariant kernels. They include a vast majority of the commonly used kernels such as fully connected NTK, RBF and inner product kernels as they only depend on the squared Euclidean norm of xx and x′x^{\prime} and the inner product x⊤x′x^{\top}x^{\prime}. We first present one set of assumptions on the kernel functions that are sufficient for our main results involving local regularity of kk around the sphere where the data concentrates in high dimensions.

Rotational invariance and local power series expansion: The kernel function kk is rotationally invariant and there is a function gg such that k(x.x′)=g(∥x∥22,∥x′∥22,x⊤x′)k(x.x^{\prime})=g(\|x\|_{2}^{2},\|x^{\prime}\|_{2}^{2},x^{\top}x^{\prime}). Furthermore, gg can be expanded as a power series of the form

We show in Corollary 3.2 that the kernels for which our main results hold cover a broad range of commonly studied kernels in practice. In particular, Theorem 3.1 also holds for α\alpha-exponential kernels defined as k(x,x′)=exp⁡(−∥x−x′∥2α)k(x,x^{\prime})=\exp(-\|x-x^{\prime}\|_{2}^{\alpha}) for α∈(0,2)\alpha\in(0,2) even though we could not yet show that they satisfy Assumptions A.1-A.2. In Appendix C, we show how the proof of Theorem 3.1 crucially relies on the rotational invariance Assumption C.1 and the fact that the eigenvalues of the kernel matrix KK are asymptotically lower bounded by a positive constant. Both conditions are also satisfied by the α\alpha-exponential kernels and the separate treatment is purely due to different proof technique used to lower bound the kernel matrix eigenvalues.

We impose the following assumptions on the data distribution.

Covariance model: We assume that the input data distribution is from one of the following sets

High dimensional regime: We assume that the effective dimension grows with the sample size nn s.t. deff/nβ→cd_{\text{eff}}/n^{\beta}\to c for some β,c>0\beta,c>0.

and parameterize the scaling by sequence of parameters τ\tau dependent on nn. In Section 3.1 we study the standard scaling τdeff→c>0\frac{\tau}{d_{\text{eff}}}\to c>0, before discussing τdeff→0\frac{\tau}{d_{\text{eff}}}\to 0 and τdeff→∞\frac{\tau}{d_{\text{eff}}}\to\infty respectively in Section 3.2, where we show that the polynomial approximation barrier is not a consequence of the standard scaling.

Main Results

We now present our main results that hold for a wide range of distributions and kernels and show that kernel methods can at most consistently learn low-degree polynomials. Section 3.1 considers the case τ≍deff\tau\asymp d_{\text{eff}} while Section 3.2 provides lower bounds for the regimes τdeff→0\frac{\tau}{d_{\text{eff}}}\to 0 and τdeff→∞\frac{\tau}{d_{\text{eff}}}\to\infty.

The bias of the kernel estimators f^λ\hat{f}_{\lambda} is asymptotically almost surely lower bounded for any ϵ>0\epsilon>0,

Corollary 3.2 gives examples for kernels that satisfy the assumptions of the theorem. The almost sure statements refers to the sequence of matricies X of random vectors xix_{i} as n→∞n\to\infty, but also hold true with probability ≥1−n2exp⁡(−Cϵ′log⁡(n)(1+ϵ′))\geq 1-n^{2}\exp(-C_{\epsilon^{\prime}}\log(n)^{(1+\epsilon^{\prime})}) over the draws of X (see Lemma C.1 for further details).

The attentive reader might notice that Ghorbani et al. achieve a lower barrier m=⌊1/β⌋m=\lfloor 1/\beta\rfloor for their specific setting which implies that our results are not tight. However, this is not the main focus in this work as we primarily intend to demonstrate that the polynomial approximation barrier persists for general covariance model data distributions (Ass. B.1-B.2) and hence asymptotic consistency in high-dimensional regimes is at most reachable for commonly used rotationally invariant kernels if the ground truth function is a low degree-polynomial. We leave tighter bounds as interesting future work. Finally, we remark that as β→0\beta\to 0 we enter again a classical asymptotic regime where nn is much larger compared to dd and hence m→∞m\to\infty. This underlines the difference between classical and high-dimensional asymptotics and shows that our results are only meaningful in the latter.

We now present a short proof sketch to provide intuition for why the polynomial approximation barrier holds. The full proof can be found in Appendix C.

The proof of the main theorem is primarily based on the concentration of Lipschitz continuous functions of vectors with i.i.d entries. In particular, we show in Lemma C.1 that

where we use tr⁡(Σd)≍nβ\operatorname{tr}(\Sigma_{d})\asymp n^{\beta}. Furthermore, Assumption A.1and hence the rotational invariance of the kernel function kk, implies that for inner product kernels where gjg_{j} are constants,

with θ\theta some constant such that 1<θ<(m+1)β21<\theta<(m+1)\frac{\beta}{2} that exists because m≥⌊2/β⌋m\geq\lfloor 2/\beta\rfloor. Hence, as n→∞n\to\infty, kτ(xi,X)k_{\tau}(x_{i},X) converge to low-degree polynomials. Using the closed form solution of f^λ\hat{f}_{\lambda} based on the representer theorem we can hence conclude the first statement in Theorem 3.1 if K+λI≻cIK+\lambda I\succ cI for some constant c>0c>0. The result follows naturally for ridge regression with non vanishing λ>0\lambda>0. However for the minimum norm interpolator, we need to show that the eigenvalues of the kernel matrix KK themselves are asymptotically lower bounded by a positive non-zero constant. This follows from the additional assumption in Theorem 3.1 and the observation that (X⊤X)∘j′→In(\textbf{X}^{\top}\textbf{X})^{\circ j^{\prime}}\to I_{n} in operator norm with ∘\circ being the Hadamard product. Finally, the case where gjg_{j} depend on ∥xi∥22\|x_{i}\|_{2}^{2} and ∥xj∥22\|x_{j}\|_{2}^{2} requires a more careful analysis and constitute the major bulk of the proof in Appendix C. ∎

The assumptions in Theorem 3.1 cover a broad range of commonly used kernels, including the ones in previous works. The following corollary summarizes some relevant special cases

The α\alpha-exponential kernel for α∈(0,2]\alpha\in(0,2], including Laplace (α=1\alpha=1) and the Gaussian (α=2\alpha=2) kernels

The fully-connected NTK of any depth with regular activation functions including the ReLU activation σ(x)=max⁡(0,x)\sigma(x)=\max(0,x)

The precise regularity conditions of the activation functions for the NTK and the proof of the corollary can be found in Appendix C.3.

In particular, the assumptions hold for instance for the α\alpha-exponential kernels with α∈(0,2)\alpha\in(0,2) (see ) and the popular Matern RBF kernels with ν<2\nu<2. The proof of the theorem can be found in Appendix D.1 and is based on the flat limit literature on RBFs . Finally, we remark that Theorem 3.3 applies only for the interpolating estimator f^0\hat{f}_{0}. However, because the ridge penalty is well known to regulate the bias-variance tradeoff, we argue that the bias increases with the ridge penalty and hence attains its minimum at λ=0\lambda=0.

3 Discussion of theoretical results

As a consequence of our unifying treatment, we cover most of the settings studied in the existing literature and show that the polynomial approximation barrier (5) is neither restricted to a few specific distributions nor to a particular choice of the scaling or the eigenvalue decay of the kernel function. We summarize our setting alongside those of previous works in Table 1.

The Assumptions B.1-B.2 allow very general distributions and include the ones in the current literature. In particular, we also cover the settings studied in the papers and hence put their optimistic conclusions into perspective. Besides, our results hold true for arbitrary covariance matrices and only depend on the growth rate of the effective dimension, but are independent of the explicit structure of the covariance matrix. This stands in contrast to linear regression where consistency can (only) be guaranteed for spiky covariance matrices .

Our results do not only apply for the standard choice of the scaling τ/deff→c>0{\tau/d_{\text{eff}}\to c>0}, but also apply to general RBF kernels in the flat limit scaling, i.e. where τ→∞\tau\to\infty. This case is particularly important since this is where we empirically find the bias to attain its minimum in Figure 3(c). We therefore conjecture that the polynomial approximation barrier cannot be overcome with different choices of the scaling.

Furthermore, by explicitly showing that the polynomial barrier persists for all α\alpha-exponential kernels with α∈(0,2]\alpha\in(0,2] (that have vastly different eigenvalue decay rates), we provide a counterpoint to previous work that suggests consistency for α→0\alpha\to 0. In particular, the paper proves minimax optimal rates for the Nadaraya-Watson estimators with singular kernels for fixed dimensions and empirical work by suggests that spikier kernels have more favorable performances. Our results suggest that in high dimensions, the effect of eigenvalue decay (and hence “spikiness”) may be dominated by asymptotic effects of rotationally invariant kernels. We discuss possible follow-up questions in Section 5. As a result, we can conclude that the polynomial approximation barrier is a rather general phenomenon that occurs for commonly used rotationally invariant kernels in high dimensional regimes. For ground truths that are inherently higher-degree polynomials that depend on all dimensions, our theory predicts that consistency of kernel learning with fully-connected NTK, standard RBF or inner product kernels is out of reach if the data is high-dimensional. In practice however, it is possible that not all dimensions carry equally relevant information. In Section 4.3 we show how feature selection can be used in such settings to circumvent the bias lower bound. On the other hand, for image datasets like CIFAR-10 where the ground truth is a complex function of all input dimensions, kernels that incorporate convolutional structures (such as CNTK or compositional kernels ) and hence break the rotational symmetry can perform quite well.

Experiments

In this section we describe our synthetic and real-world experiments to further illustrate our theoretical results and underline the importance of feature selection in high dimensional kernel learning.

In Figure 1, we demonstrate how the Hilbert norm of the simple sparse linear function f⋆(x)=x(1)f^{\star}(x)=x_{(1)} grows with dimension dd. We choose the scaling τ=d\tau=d and consider the Hilbert space induced by the scaled Gaussian kτ(x,x′)=exp⁡(−∥x−x′∥2τ)k_{\tau}(x,x^{\prime})=\exp(-\frac{\|x-x^{\prime}\|^{2}}{\tau}), Laplace kτ(x,x′)=exp⁡(−∥x−x′∥2τ)k_{\tau}(x,x^{\prime})=\exp(-\frac{\|x-x^{\prime}\|_{2}}{\sqrt{\tau}}) and exponential inner product kτ(x,x′)=exp⁡(−xTx′τ)k_{\tau}(x,x^{\prime})=\exp(-\frac{x^{T}x^{\prime}}{\tau}) kernels. To estimate the norm, we draw 75007500 i.i.d. random samples with noiseless observations from the uniform distribution on X=d\mathcal{X}=^{d}.

2 Illustration of the polynomial approximation barrier

We now provide details for the numerical experiments in Figure 2,3 that illustrate the lower bounds on the bias in Theorem 3.1 and Theorem 3.3. For this purpose, we consider the following three different distributions that satisfy the assumptions of the theorems and are covered in previous works

While our theoretical results only show lower bounds on the bias, Figure 2(a) and 2(b) suggest that as β\beta increases, the bias in fact stepwise aligns with the best polynomial of lower and lower order. This has also been shown in for the uniform distribution from the product of two spheres. Indeed, with decreasing β\beta, for the cubic polynomial in Figure 2(b) we first learn linear functions (first descent in the curve). Since the best degree 2 polynomial approximation of f2∗f^{*}_{2} around is a linear function, the curve then enters a plateau before descending to zero indicating that we successfully learn the ground truth.

3 Feature selection for high-dimensional kernel learning

In this section, we demonstrate how the polynomial approximation barrier limits the performance in real world datasets and how one may overcome this issue using feature selection. How our theory motivates feature selection can be most cleanly illustrated for complex ground truth functions that only depend on a number of covariates that is much smaller than the input dimension (which we refer to as sparse)for more general functions that are not sparse we still expect a U-curve, albeit potentially less pronounced, whenever the function is not a low-degree polynomial of order 2/β2/\beta:

Based on our theoretical results we expect that for sparse ground truths, the bias follows a U-shape as dimension increases: until all relevant features are included, the bias first decreases before it then starts to decrease due to the polynomial approximation barrier that holds for large dd when asymptotics start to kick in. Since recent work shows that the variance vanishes in high-dimensional regimes (see e.g. ), we expect the risk to follow a U-shaped curve as well. Hence, performing feature selection could effectively yield much better generalization for sparse ground truth functions. We would like to emphasize that although the described curve may ring a familiar bell, this behavior is not due to the classical bias-variance trade-off, since the U-shaped curve can be observed even in the noiseless case where we have zero variance. We now present experiments that demonstrate the U-shape of the risk curve for both synthetic experiments on sparse ground truths and real-world data. We vary the dimensionality dd by performing feature selectionwe expect other approaches to incorporate sparsity such as automatic relevance determination to yield a similar effect using the algorithm proposed in the paper . In order to study the impact of high-dimensionality on the variance, we add different levels of noise to the observations.

For the real-world experiments we are not able to decompose the risk to observe separate trends of the bias and the variance. However, we can argue that the ridge estimator should perform significantly better than the minimum norm interpolator whenever the variance dominates. Vice versa, if their generalization errors are close, the effect of the variance vanishes. This is due to the fact that the ridge penalty decreases the variance and increases the bias. Hence, for real-world experiments, we use the comparison between the risks of the ridge estimator and minimum norm interpolator to deduce that the bias is in fact dominating the risk for high dimensions.

We now explore the applicability of our results on real-world data where the assumptions of the theorems are not necessarily satisfied. For this purpose we select datasets where the number of features is large compared to the number of samples. In this section we show results on the regression dataset residential housing (RH) with n=372n=372 and d=107d=107 to predict sales prices from the UCI website . Further datasets can be found in Appendix A.3. In order to study the effect of noise, we generate an additional dataset (RH-2) where we add synthetic i.i.d. noise drawn from the uniform distribution on [−1/2,1/2][-1/2,1/2] to the observations. The plots in Figure 4 are then generated as follows: we increase the number of features using a greedy forward selection procedure (see Appendix A.3 for further details ). We then plot the risk achieved by the kernel ridge and ridgeless estimate using the Laplace kernel on the new subset of features.

Figure 4(b) shows that the risks of the minimum norm interpolator and the ridge estimator are identical, indicating that the risk is essentially equivalent to the bias. Hence our first conclusion is that, similar to the synthetic experiment, the bias follows a U-curve. For the dataset RB-2 in Figure 4(c), we further observe that even with additional observational noise, the ridge and ridgeless estimator converge, i.e. the bias dominates for large dd. We observe both trends in other high-dimensional datasets discussed in Appendix A.3 as well. As a consequence, we can conclude that even for some real-world datasets that do not necessarily satisfy the conditions of our bias lower bound, feature selection is crucial for kernel learning for noisy and noiseless observations alike. We would like to note that this conclusion does not necessarily contradict empirical work that demonstrates good test performance of RBFs on other high-dimensional data such as MNIST. In fact, the latter only suggests that linear or polynomial fitting would do just as well for these datasets which has indeed been suggested in .

Conclusion and future work

Kernel regression encourages estimators to have a certain structure by means of the RKHS norm induced by the kernel. For example, the eigenvalue decay of α\alpha-exponential kernels results in estimators that tend to be smooth (i.e. Gaussian kernel) or more spiky (i.e. small α<1\alpha<1). A far less discussed fact is that many kernels implicitly incorporate additional structural assumptions. For example, rotational invariant kernels are invariant under permutations and hence treat all dimensions equally Even though rotational invariance is a natural choice when no prior information on the structure of the ground truth is available, this paper shows that the corresponding inductive bias in high dimensions is in fact restricting the average estimator to a polynomial. In particular, we show in Theorems 3.1 and 3.3 that the lower bound on the bias is simply the projection error of the ground truth function onto the space of polynomials of degree at most 2⌊2/β⌋2\lfloor 2/\beta\rfloor respectively ⌊2/β⌋\lfloor 2/\beta\rfloor. Apart from novel technical insights that result from our unified analysis (discussed in Sec. 3.3), our result also opens up new avenues for future research.

Modern datasets which require sophisticated methods like deep neural networks to obtain good predictions are usually inherently non-polynomial and high-dimensional. Hence, our theory predicts that commonly used rotationally invariant kernels cannot perform well for these problems due to a high bias. In particular, our bounds are independent of properties like the smoothness of the kernel function and cannot be overcome by carefully choosing the eigenvalue decay. Therefore, in order to understand why certain highly overparameterized methods generalize well in these settings, our results suggest that it is at least equally important to understand how prior information can be incorporated to break the rotational symmetry of the kernel function. Examples for recent contributions in this direction are kernels relying on convolution structures (such as CNTK or compositional kernels ) for image datasets.

Another relevant future research direction is to present a tighter non-asymptotic analysis that allows a more accurate characterization of the estimator in practice. The presented results in this paper are asymptotic statements, meaning that they do not provide explicit bounds for fixed n,dn,d. Therefore, for given finite n,dn,d it is unclear which high-dimensional regime provides the most accurate characterization of the estimator’s statistical properties. For instance, our current results do not provide any evidence whether the estimator follows the bias lower bounds for n=dβn=d^{\beta} with β=log⁡(n)/log⁡(d)\beta=\log(n)/\log(d) or n=γdn=\gamma d. We remark that the methodology used to prove the statements in this paper could also be used to derive non-asymptotic bounds, allowing us to further investigate this problem. However, we omitted such results in this paper for the sake of clarity and compactness of our theorem statements and proofs and leave this for future work.

Acknowledgments

We would like to thank the anonymous reviewers for their helpful feedback, Xiao Yang for insightful discussions on the flat limit kernel and Armeen Taeb for his comments on the manuscript.

References

Appendix A Experiments

This section contains additional experiments not shown in the main text. Our code is publicly available at https://www.github.com/DonhauserK/High-dim-kernel-paper/

In this section, we provide additional experiments that discuss Theorem 3.1. In particular, we investigate kernels beyond the Laplace kernel and study the behaviour of the bias with respect to β\beta when dd is fixed and nn varies. The experimental setting is the same as the one in Section 4.2.

Instead of comparing the bias curves for different input distribution as in Figure 2(a), Figure 5 shows the bias with respect to β\beta for the α\alpha-exponential kernel, i.e. k(x,x′)=exp⁡(−∥x−x′∥2α)k(x,x^{\prime})=\exp(-\|x-x^{\prime}\|_{2}^{\alpha}), for different choices of α\alpha and hence for kernels with distinct eigenvalue decays (α=2\alpha=2 results in an exponential eigenvalue decay while α<2\alpha<2 in a polynomial eigenvalue decay). Clearly, we can see that the curves transition at a similar value for β\beta which confirms the the discussion of Theorem 3.1 in Section 3.3 where we argue that the polynomial approximation barrier occurs independently of the eigenvalue decay.

Figure 6 shows the bias of the minimum norm interpolant B(f^0)\textbf{B}(\hat{f}_{0}) normalized by B(0)\textbf{B}(0) for the ground truth function f⋆(x)=2x(1)3f^{\star}(x)=2x_{(1)}^{3} and the Laplace kernel as in Section 4.2 with τ=deff\tau=d_{\text{eff}}. We observe that the asymptotics already kick in for d≈40d\approx 40 since all curves for d≥40d\geq 40 resemble each other. This confirms the the trend in Figure 2(b).

A.2 Feature selection - Synthetic

The goal of this experiment is to compare the bias variance trade-off of ridge regression and minimum norm interpolation. We use the same experimental setting as the ones used for Figure 4(a) (see Section 4.3). We set the bandwidth to τ=deff\tau=d_{\text{eff}} and choose the ridge parameter λ\lambda using 55-fold cross validation. While for small dimensions dd, ridge regularization is crucial to achieve good performance, the bias becomes dominant as the dimension grows and the difference of the risks of both methods shrinks. This aligns well with Theorem 3.1 which predicts that the bias starts to increase with dd for fixed nn once we enter the asymptotic regime.

A.3 Feature selection - Real world

We now present details for our real world experiments to emphasize the relevance of feature selection when using kernel regression for practical applications, as discussed in Section 4.3.

The residential housing regression data set from the UCI website where we predict the sales prices and construction costs given a variety of features including the floor area of the building or the duration of construction.

The ALLAML classification data set from the ASU feature selection website where we classify patients with acute myeloid leukemia (AML) and acute lymphoblastic leukemia (ALL) based on features gained from gene expression monitoring via DNA microarrays.

The CLL_SUB_111 classification dataset from the ASU feature selection web-page where we classify genetically and clinically distinct subgroups of B-cell chronic lymphocytic leukemia (B-CLL) based on features consisting of gene expressions from high density oligonucleotide arrays. While the original dataset contains three different classes, we only use the classes 2 and 3 for our experiment to obtain a binary classification problem.

Because the number of features in the ALLAML and CLL_SUB_111 datasets massively exceed the number of samples, we run the feature selection algorithm in and pre-select the best 100100 features chosen by the algorithm. In order to reduce the computational expenses, we run the algorithm in batches of 20002000 features and iteratively remove all features except for the best 200200 features chosen by the algorithm. We do this until we reduce the total number of features to 20002000 and then select in a last round the final 100100 features used for the further procedure. Reducing the amount of features to 100100 is important for the computational feasibility of greedy forward features selection in our experiments. The properties of the datasets are summarized in Table 2.

Results of the experiments: The following figures present the results of our experiments on all datasets except for the ones predicting the sales prices in the residential housing dataset, which we presented in Figure 4(b),4(c) in the main text. Similar to the observations made in Section 4.3, Figure 8,9,10 show that the risk reaches its minimum around d≈25d\approx 25, with significant differences to the right at d≈100d\approx 100. In particular, this holds for both ridge regression and interpolation, which again shows that the bias becomes dominant as the dimension increases. Surprisingly, we also note that the relevance of ridge regularization seems to be much smaller for classification tasks than regression tasks.

Appendix B Bounded Hilbert norm assumption

This section gives a formal statement of Lemma 2.1. We begin with the conditions under which the lemma holds. We consider tensor product kernels of the form

For any j>mj>m, define fj=1f_{j}=1. First, we note that the proof follows trivially if any of the fjf_{j} is not contained in the RKHS induced by qq since this implies that the Hilbert norm ∥f∥Hk=∞\|f\|_{\mathcal{H}_{k}}=\infty. Hence, we can assume that for all jj, fjf_{j} is contained in the RKHS for all dd. Furthermore, because kk is a product kernel, we can write ∥f∥Hk=∏j=1d∥fj∥Hq\|f\|_{\mathcal{H}_{k}}=\prod_{j=1}^{d}\|f_{j}\|_{\mathcal{H}_{q}} where ∥.∥Hq\|.\|_{\mathcal{H}_{q}} is the Hilbert norm induced by qq on X\mathcal{X}. Because we are only interested to see whether the sequence of Hilbert norms diverge, without loss of generality we can assume that m=1m=1, and hence,

Next, by Mercer’s theorem there exists an orthonormal eigenbasis {ϕi}i=1∞\{\phi_{i}\}_{i=1}^{\infty} in L2(X,μ)\mathcal{L}_{2}(\mathcal{X},\mu) with corresponding eigenvalues {λi}i=1∞\{\lambda_{i}\}_{i=1}^{\infty} such that for any g∈Hqg\in\mathcal{H}_{q}, ∥g∥Hq=∑i=1∞(⟨f,ϕi⟩)2λi\|g\|_{\mathcal{H}_{q}}=\sum_{i=1}^{\infty}\frac{(\langle f,\phi_{i}\rangle)^{2}}{\lambda_{i}}, where ⟨f,ϕi⟩=∫ϕi(x)f(x)dμ(x)\langle f,\phi_{i}\rangle=\int\phi_{i}(x)f(x)d\mu(x). Note that because the kernel qq depends on dd, λi\lambda_{i} and ϕi\phi_{i} also depend on dd. Next, because by assumption f(x)=1f(x)=1 is contained in the RKHS, there exists αi\alpha_{i} such that for every x∈Xx\in\mathcal{X}, 1=∑i=1∞αiϕi(x)1=\sum_{i=1}^{\infty}\alpha_{i}\phi_{i}(x) and ∑i=1∞αi2=1\sum_{i=1}^{\infty}\alpha_{i}^{2}=1. Furthermore,

Furthermore, there exists βi\beta_{i} such that f1(x)=∑i=1∞βiϕi(x)f_{1}(x)=\sum_{i=1}^{\infty}\beta_{i}\phi_{i}(x). Again, because we are only interested to see whether the sequence of Hilbert norms diverge, without loss of generality we can assume that ∑i=1∞βi2=1\sum_{i=1}^{\infty}\beta_{i}^{2}=1 and hence also ∥fj∥Hq≥1\|f_{j}\|_{\mathcal{H}_{q}}\geq 1.

Together with the fact that αjd2→1\alpha_{j_{d}}^{2}\to 1 it then follows that ∑i≠jdβi2\sum_{i\neq j_{d}}\beta_{i}^{2} has to be asymptotically lower bounded by some positive non-zero constant c2c_{2} and hence

This contradicts the assumption that ∥f1∥Hq\|f_{1}\|_{\mathcal{H}_{q}} is upper bounded by some constant for every dd. Hence, we are only left with the case where ∥1∥Hq≥c>1\|1\|_{\mathcal{H}_{q}}\geq c>1, however, this case diverges due to Equation 9. Hence, the proof is complete. ∎

Appendix C Proof of Theorem 3.1

Before presenting the proof of the (generalized) theorem, we first state the key concentration inequalities used throughout the proof. It is an extension of Lemma A.2 in the paper , which iteself is a consequence of the concentration of Lipschitz continuous functions of i.i.d random vectors.

Then, there exists some constant C>0C>0 such that for nn sufficiently large,

In particular, the event EX\mathcal{E}_{\textbf{X}} holds almost surely with respect to the sequence of data sets X as n→∞n\to\infty, that is the probability that for infinitely many nn, EX\mathcal{E}_{\textbf{X}} does not hold, is zero.

The proof of the lemma can be found in Section E.1.

The proof of the Theorem is primarily separated into two parts

We first state Theorem C.2 which shows under the weaker Assumption C.1 that the results of 3.1 hold for the ridge estimate f^λ\hat{f}_{\lambda} for non-vanishing λ>0\lambda>0 or the ridgeless estimate whenever the eigenvalues of KK are asymptotically lower bounded.

We finish the proof for the ridgeless estimate by invoking Theorem C.2 and showing that KK indeed has asymptotically lower bounded eigenvalues under the stricter assumptions A.1-A.3 imposed in Theorem 3.1.

For the clarity we denote with A.3 the β\beta-dependent assumptions in Theorem 3.1

β\beta-dependent assumptions: gig_{i} is (⌊2/β⌋+1−i)(\lfloor 2/\beta\rfloor+1-i)-times continuously differentiable in a neighborhood of (1,1)(1,1) and there exists j′>⌊2/β⌋j^{\prime}>\lfloor 2/\beta\rfloor such that gj′(1,1)>0g_{j^{\prime}}(1,1)>0.

We start by introducing the following weaker assumptions that allows us to jointly treat α\alpha-exponential kernels and kernels satisfying Assumption A.1-A.3 when the kernel eigenvalues are lower bounded in Theorem C.2. Note that this assumption implies that the kernel is rotationally invariant.

The kernel function kk is rotationally invariant and there exists a function gg such that k(x.x′)=g(∥x∥22,∥x′∥22,x⊤x′)k(x.x^{\prime})=g(\|x\|_{2}^{2},\|x^{\prime}\|_{2}^{2},x^{\top}x^{\prime}). Furthermore, gg can be expanded as a power series of the form

with m=⌊2/β⌋m=\lfloor 2/\beta\rfloor that converges in a neighborhood N(δ,δ′)N(\delta,\delta^{\prime}) of the sphere for some δ,δ′>0\delta,\delta^{\prime}>0 and where gig_{i} is (⌊2/β⌋+1−i)(\lfloor 2/\beta\rfloor+1-i)-times continuously differentiable in an neighborhood of (1,1)(1,1) and the remainder term rr is a continuous function around the point (1,1,0)(1,1,0).

The bias of the kernel estimators f^λ\hat{f}_{\lambda} is asymptotically lower bounded, for any ϵ>0\epsilon>0,

We can find a polynomial pp such that for any ϵ,ϵ′>0\epsilon,\epsilon^{\prime}>0, there exists C>0C>0 such that asymptotically with probability ≥1−n2exp⁡(−C(log⁡(n))1+ϵ′)\geq 1-n^{2}\exp(-C(\log(n))^{1+\epsilon^{\prime}}) over the draws of XX,

The proof of this theorem can be found in Section C.1. Theorem C.2 states Theorem 3.1 under the assumption that (K+λI)(K+\lambda I) has asymptotically lower bounded eigenvalues and the weaker Assumption C.1. For the proof of Theorem 3.1, it remains to show that Assumptions A.1-A.3 of Theorem 3.1 and the α\alpha-exponential kernel both

induce kernel matrices with almost surely asymptotically positive lower bounded eigenvalues

Point (a) is relatively simple to prove and deferred to Section E.6. The bulk of the work in fact lies in showing (b) separately for the case for A.1-A.3 and α\alpha-exponential kernels with α∈(0,2)\alpha\in(0,2) in the following two propositions, as these two cases require two different proof techniques.

where λmin⁡(K)\lambda_{\min}(K) is the minimum eigenvalue of the kernel matrix KK.

Assume that the Assumptions B.1-B.2 hold true. Then, the minimum eigenvalue of the kernel matrix of the α\alpha-exponential kernel with α∈(0,2)\alpha\in(0,2) is lower bounded by some positive constant almost surely as n→∞n\to\infty.

The proof of the Propositions C.3 and C.4 can be found in the Sections C.2.1 and C.2.2 respectively which concludes the proof of the theorem.

The almost sure statement in Proposition C.4 can also be replaced with an in probability statement as in Lemma C.1 and hence also the statements in Theorem 3.1.

As a result of Lemma C.1 it is sufficient to condition throughout the rest of this proof on the intersection of the events EX\mathcal{E}_{\textbf{X}} and the event where the eigenvalues of the kernel matrix KK are lower bounded by a positive constant.

The idea of the proof is to decompose the analysis into the term emerging from the error in the high probability region EZ∣Z\mathcal{E}_{Z|\mathbf{Z}} and the error emerging from the low probability region EZ∣Zc\mathcal{E}_{Z|\mathbf{Z}}^{\mathsf{c}}. The proof essentially relies on the following lemma.

We can construct a polynomial pp of degree ≤m\leq m such that for n→∞n\to\infty,

The proof of the lemma can be found in Section E.2. As a result, Equation 16 follows immediately and Equation 17 is a consequence of

The first two terms vanish due to Lemma C.6. To see that the third term vanishes, note that for nn sufficiently large,

Next, the lower bound for the bias. Due to Lemma C.8, we have that

and hence, for any γ1>0\gamma_{1}>0 and nn sufficiently large,

Thus, the result follows from the definition of the infimum. ∎

C.2 Proofs for the lower bound of the eigenvalues

We use the same notaiton as used in the proof of Theorem C.2. As a result of Lemma C.1 it is sufficent to condition on EX\mathcal{E}_{\textbf{X}} throughout the rest of this proof. The proof follows straight forwardly from the following Lemma C.7 which gives an asymptotic description of the kernel matrix KK based on a similar analysis as the one used in the proof of Theorem 2.1 and 2.2 in the paper . In essence, it is again a consequence of the concentration inequality from Lemma C.1 and the stronger Assumption A.1-A.3 and in particular the power series expansion of gg. We denote with ∘i\circ i the ii-times Hadamard product.

Given that the assumption in Proposition C.3 hold. For m=⌊2/β⌋m=\lfloor 2/\beta\rfloor,

where GgqG_{g_{q}} is the positive semi-definite matrix with entries (Ggq)i,j=gl(∥zi∥22,∥zj∥22)(G_{g_{q}})_{i,j}=g_{l}(\|z_{i}\|_{2}^{2},\|z_{j}\|_{2}^{2}).

The proof of the lemma can be found in Section E.3. The proof of Proposition C.3 then follows straight forwardly when using Schur’s product theorem which shows that

where we use that gig_{i} is positive semi-definite by Assumption A.1. To see that the eigenvalues are lower bounded, we thus simply need to show that g(1,1,1)−∑q=0mgq(1,1)>0g(1,1,1)-\sum_{q=0}^{m}g_{q}(1,1)>0. This holds because the positive semi-definiteness of gqg_{q} implies that gq(1,1)≥0g_{q}(1,1)\geq 0 and hence g(1,1,1)=∑q=0∞gq(1,1)g(1,1,1)=\sum_{q=0}^{\infty}g_{q}(1,1) is a sum of positive coefficients and because by Assumption A.3 there exists j′>⌊2/β⌋j^{\prime}>\lfloor 2/\beta\rfloor such that gj′(1,1)>0g_{j^{\prime}}(1,1)>0. Hence, there exists a positive constant c>0c>0 such that λmin⁡(M)≥c\lambda_{\min}(M)\geq c. We can conclude the proof when applying Lemma C.7, which implies that λmin⁡(K)→λmin⁡(M)\lambda_{\min}(K)\to\lambda_{\min}(M) as n→∞n\to\infty. ∎

C.2.2 Proof of Proposition C.4

We use the same notaiton as used in the proof of Theorem C.2 and define DαD_{\alpha} to be the n×nn\times n matrix with entries (Dα)i,j=dα(zi,zj):=∥zi−zj∥2α(D_{\alpha})_{i,j}=d_{\alpha}(z_{i},z_{j}):=\|z_{i}-z_{j}\|_{2}^{\alpha}. We separate the proof into two steps. In a first step, we decompose kk in the terms

The proof of the lemma can be found in Section E.4. In particular, note that we can use the same argument as used in Lemma C.1 to show that there exists almost surely over the draws of Z as n→∞n\to\infty an additional vector z0z_{0}, such that for any two vectors z,z′∈Z∪{z0}z,z^{\prime}\in\textbf{Z}\cup\{z_{0}\},

Throughout the rest of this proof, we conditioned on the event EX\mathcal{E}_{\textbf{X}} and the additional event that Equation (21) holds, and remark that the intersection of these two events holds true almost surely as n→∞n\to\infty. It is then straight forward to show that the eigenvalues of the matrix

where we have used that 1T(vT−1⊤v)=01^{T}\begin{pmatrix}v^{T}\\ -1^{\top}v\end{pmatrix}=0. As a result, we can see that

Next, using Lemma C.1 and the fact that dα(x,x′)=2α/2+O(∣∣x−x′∣∣222−1)d_{\alpha}(x,x^{\prime})=2^{\alpha/2}+O(\frac{\lvert\lvert x-x^{\prime}\rvert\rvert_{2}^{2}}{2}-1), we can see that γ≥exp⁡(2α/2/2)>1{\gamma\geq\exp(2^{\alpha/2}/2)>1}. Hence, it is sufficient to show that ψψ⊤−12e2 11⊤\psi\psi^{\top}-\frac{1}{2e^{2}}~{}11^{\top} is positive semi-definite. This is true if and only if 1⊤ψψ⊤1≥12e21⊤11⊤11^{\top}\psi\psi^{\top}1\geq\frac{1}{2e^{2}}1^{\top}11^{\top}1, which is equivalent to saying that (∑i=1nexp⁡(−dα(zi,z0)))2≥n22e2\left(\sum_{i=1}^{n}\exp(-d_{\alpha}(z_{i},z_{0}))\right)^{2}\geq\frac{n^{2}}{2e^{2}}. Using again the same argument as for γ\gamma, we can see that max⁡i∣2α/2−dα(zi,z0)∣→0\underset{i}{\max}\left|2^{\alpha/2}-d_{\alpha}(z_{i},z_{0})\right|\to 0 for any ii, which completes the proof. ∎

C.3 Proof of Corollary 3.2

First, note that the Assumption A.1-A.3 straight forwardly hold true for the exponential inner product kernel with k(x,x′)=exp⁡(x⊤x′)=∑j=0∞1j!(x⊤x′)jk(x,x^{\prime})=\exp(x^{\top}x^{\prime})=\sum_{j=0}^{\infty}\frac{1}{j!}(x^{\top}x^{\prime})^{j} and for the Gaussian kernel with

Next, note that the α\alpha-exponential kernel with α<2\alpha<2 is already explicitly covered in Theorem 3.1. Hence, the only thing left to show is that Theorem 3.1 also applies to ReLU-NTK.

Assume that the activation function σ\sigma is kk-homogeneous and both the activation function and its derivative possess a Hermite-polynomial series expansion (see ) where there exits j′≥⌊2/β⌋j^{\prime}\geq\lfloor 2/\beta\rfloor such that the j′j^{\prime}-th coefficient aj′≠0a_{j^{\prime}}\neq 0. Then, the NTK satisfies the Assumption A.1-A.3 and hence Theorem 3.1 applies.

In fact, we can easily see that any non linear activation function which is homogeneous and both the activation function and its derivative possesses a Hermite polynomial extension satisfies the assumptions in Proposition C.9. In particular, this includes the popular ReLU activation function σ(x)=max⁡(x,0)\sigma(x)=\max(x,0) where the explicit expression for the Hermite polynomial extension can be found in . ∎

The only thing left to show is Assumption A.3. While we have already shown that gjg_{j} are smooth in a neighborhood of (1,1)(1,1), we still need to show that there exists j′>⌊2/β⌋j^{\prime}>\lfloor 2/\beta\rfloor such that gj′(1,1)>0g_{j^{\prime}}(1,1)>0. However, this follows from the fact that by assumption there exists j′>⌊2/β⌋j^{\prime}>\lfloor 2/\beta\rfloor such that aj′≠0a_{j^{\prime}}\neq 0 where aja_{j} are the Hermite coefficients of the activation function σ\sigma. ∎

Appendix D Different scalings τ𝜏\tau

In this section, we present results for different choices of the scaling beyond the standard choice τ≍deff\tau\asymp d_{\text{eff}}. In Subsection D.1, we give a proof of Theorem 3.3 describing the flat limit, i.e. the limit of the interpolant where for any fixed n,dn,d, τ→∞\tau\to\infty. Furthermore, in order to get a more comprehensive picture, we additionally present straight forward results for other choices of τ\tau in Section D.2.

First, although the limit lim⁡τ→∞K−1\lim_{\tau\to\infty}K^{-1} does not exists, we can apply Theorem 3.12 in to show that the flat limit interpolator fFL:=lim⁡τ→∞f^0f_{\textrm{FL}}:=\lim_{\tau\to\infty}\hat{f}_{0} of any kernel satisfying the assumption in Theorem 3.3 exists and has the form

Furthermore, for the α\alpha-exponential kernel, we use Theorem 2.1 in to show that it satisfies the assumptions imposed on the eigenvalue decay in Theorem 3.3.

The estimator fFLf_{\textrm{FL}} is also called the polyharmonic spline interpolator. This estimator is invariant under rescalings of the input data which is also the reason why we can rescale the input data by deff\sqrt{d_{\text{eff}}}, i.e. consider zi=xi/deffz_{i}=x_{i}/\sqrt{d_{\text{eff}}} as input data points.

and 1TDα−11>01^{T}D_{\alpha}^{-1}1>0. In particular, this allows us to use the block matrix inverse to show that

Furthermore, by Lemma C.1 and the fact that ∥zi−Z∥22=zi⊤zi+Z⊤Z−2zi⊤Z\|z_{i}-Z\|_{2}^{2}=z_{i}^{\top}z_{i}+Z^{\top}Z-2z_{i}^{\top}Z, we can see that for q=⌊2/β⌋q=\lfloor 2/\beta\rfloor, nO((12∥zi−Z∥22−1)q+1)→0nO\left(\left(\frac{1}{2}\|z_{i}-Z\|_{2}^{2}-1\right)^{q+1}\right)\to 0. Hence, assuming that the absolute eigenvalues ∣λi(A)∣|\lambda_{i}(A)| of AA are all upper bounded by a non-zero positive constant, we can use exactly the same argument as used in the proof of Theorem C.2 to conclude the proof.

Thus, we only need to show that the eigenvalues ∣λi(A)∣|\lambda_{i}(A)| are upper bounded. We already know from Lemma C.8 that there exists some constant c>0c>0 independent of nn, such that ∣∣Dα−1∣∣op≤c\left|\left|D_{\alpha}^{-1}\right|\right|_{\textrm{op}}\leq c. Thus, we only need to show that ∣∣Dα−111TDα−11TDα−11∣∣op\left|\left|\frac{D_{\alpha}^{-1}11^{T}D_{\alpha}^{-1}}{1^{T}D_{\alpha}^{-1}1}\right|\right|_{\textrm{op}} is almost surely upper bounded. Because Dα−111TDα−11TDα−11\frac{D_{\alpha}^{-1}11^{T}D_{\alpha}^{-1}}{1^{T}D_{\alpha}^{-1}1} is a rank one matrix, we know that

where we use that 1TDα−11>01^{T}D_{\alpha}^{-1}1>0 which we already know from the discussion above. Because by the binomial expansion, dα(zi,zj)=2α/2+O(12∥zi−zj∥22−1)d_{\alpha}(z_{i},z_{j})=2^{\alpha/2}+O(\frac{1}{2}\|z_{i}-z_{j}\|_{2}^{2}-1), Lemma C.1 implies that max⁡i≠jdα(zi,zj)→2α/2\underset{i\neq j}{\max}d_{\alpha}(z_{i},z_{j})\to 2^{\alpha/2} and hence, 1n1TDα1≳n\frac{1}{n}1^{T}D_{\alpha}1\gtrsim n, for any nn sufficiently large. Therefore, λ1≥n\lambda_{1}\geq n. Hence, there exists some constant c>0c>0 independent of nn, such that α121λ1≤c\alpha_{1}^{2}\frac{1}{\lambda_{1}}\leq c. As a consequence,

Next, note that our assumption imply that for nn sufficiently large, v⊤v−γ2≥1/2v⊤v  v^{\top}v-\gamma^{2}\geq 1/2v^{\top}v~{}~{}. Therefore,

D.2 Additional results

In this section, we present some additional results for different choices of the scaling. The results presented in this section are straight forward but provide a more complete picture for different choices of the scaling τ\tau. We use again the same notation as used in Appendix C.1.

First, we show the case where τ→0\tau\to 0. We assume that kk is the α\alpha-exponential kernel with α∈(0,2]\alpha\in(0,2], i.e. k(x,x′)=exp⁡(−∥x−x′∥2α)k(x,x^{\prime})=\exp(-\|x-x^{\prime}\|_{2}^{\alpha}).

and with probability ≥1−(n+1)2exp⁡(−C(log⁡(n))1+ϵ)\geq 1-(n+1)^{2}\exp(-C(\log(n))^{1+\epsilon}) over the draws of XX,

Hence, the result follows immediately from f^λ(X)=yT(K+λI)−1kτ(X,X)\hat{f}_{\lambda}(X)=y^{T}(K+\lambda I)^{-1}k_{\tau}(\textbf{X},X). ∎

We can also show a similar result for the case where τ→∞\tau\to\infty and λ\lambda does not vanish.

with c=f(X)T(11T+λIn)−11c=f(\textbf{X})^{T}(11^{T}+\lambda I_{n})^{-1}1

We use the same notation as in Lemma D.3. Again due to Lemma C.1, we find that

where we have used that θα/2>1\theta\alpha/2>1.

The only thing left to show is that f(X)T(11T+λIn)−11f(\textbf{X})^{T}(11^{T}+\lambda I_{n})^{-1}1 does not diverge. For this, let aIn+b11Ta\textrm{I}_{n}+b11^{T} be the inverse of 11T+λIn11^{T}+\lambda\textrm{I}_{n}. As a result of a simple computation we find that a=λ(λ+1+(n−1))+(n−1)λ(λ+1+(n−1))a=\frac{\lambda(\lambda+1+(n-1))+(n-1)}{\lambda(\lambda+1+(n-1))} and b=−1λ(λ+1+(n−1))b=\frac{-1}{\lambda(\lambda+1+(n-1))}. Hence,

Appendix E Technical lemmas

We begin with the following Lemma, which is a direct consequence of the results in Appendix A in .

i.i.d entries x(i)x_{(i)} almost surely bounded ∣x(i)∣≤c|x_{(i)}|\leq c by some constant c>0c>0 and with zero mean and unit variance.

standard normal distributed i.i.d entries.

Let MM be any symmetric matrix with ∣∣M∣∣op=1\left|\left|M\right|\right|_{\textrm{op}}=1 and let M=M+−M−M=M_{+}-M_{-} be the decomposition of MM into two positive semi-definite matrices M+M_{+} and M−M_{-} with ∣∣M+∣∣op,∣∣M−∣∣op≤1\left|\left|M_{+}\right|\right|_{\textrm{op}},\left|\left|M_{-}\right|\right|_{\textrm{op}}\leq 1. Then, there exists some positive constants C1,C2,C3C_{1},C_{2},C_{3} independent of MM such that for any r>ζ=C1/tr⁡(M+)r>\zeta=C_{1}/\operatorname{tr}(M_{+}),

Following the same argument as the one used in Corollary A.2 in , we can use Lemma E.1 to show that there exists constants C2,C3>0C_{2},C_{3}>0 such that for n→∞n\to\infty,

We now make us of the Borel-Cantelli Lemma. For any ϵ>0\epsilon>0, let r(n)=12n−β/2 (log⁡(n))(1+ϵ)/2r(n)=\frac{1}{\sqrt{2}}n^{-\beta/2}\ (\log(n))^{(1+\epsilon)/2} and note that because tr⁡(Σd)≍nβ\operatorname{tr}(\Sigma_{d})\asymp n^{\beta}, ζ\zeta decays at rate n−βn^{-\beta} and in particular, for any nn sufficiently large, r(n)/2>ζ→0r(n)/2>\zeta\to 0. Hence, we can see that there exists some constant C>0C>0 such that for any nn sufficiently large

which allows us to apply the Borel-Cantelli Lemma. Hence,

which concludes the first step of the proof.

Next, we already know from the previous discussion that for any nn sufficiently large, P(EX)≥1−n2exp⁡(−C (log⁡(n))(1+ϵ)/2)P(\mathcal{E}_{\textbf{X}})\geq 1-n^{2}\exp(-C\ (\log(n))^{(1+\epsilon)/2}). Furthermore, because XX is independently drawn from the same distribution as xix_{i}, P(EX∪EX∣X)≥1−(n+1)2exp⁡(−C (log⁡(n))(1+ϵ)/2)P(\mathcal{E}_{\textbf{X}}\cup\mathcal{E}_{X|\textbf{X}})\geq 1-(n+1)^{2}\exp(-C\ (\log(n))^{(1+\epsilon)/2}). Hence, for any nn sufficiently large,

First, note that the case where i=ji=j is clear. Let si=xi∥xi∥2s_{i}=\frac{x_{i}}{\|x_{i}\|_{2}} and zi=xitr⁡(Σd)z_{i}=\frac{x_{i}}{\sqrt{\operatorname{tr}(\Sigma_{d})}}. Since we are in the Euclidean space, the inner product is given by

Due to Equation (24), we have that ∣ ∥zi∥22−1∣≤n−β/2(log⁡(n))(1+ϵ)/2 a.s. as n→∞|~{}\|z_{i}\|_{2}^{2}-1|\leq n^{-\beta/2}(\log(n))^{(1+\epsilon)/2}~{}a.s.~{}\text{as}~{}n\to\infty and (ziTzj)2≤(n−β/2(log⁡(n))(1+ϵ)/2)2a.s. as n→∞(z_{i}^{T}z_{j})^{2}\leq(n^{-\beta/2}(\log(n))^{(1+\epsilon)/2})^{2}a.s.~{}\text{as}~{}n\to\infty. Therefore,

The rest of the proof then follows straight forwardly. ∎

We know that ff is λmax⁡(M+)/tr⁡(M+)\lambda_{\max}(M_{+})/\sqrt{\operatorname{tr}(M_{+})}-Lipschitz continuous and hence also a ∣∣M∣∣op/tr⁡(M+)\sqrt{\left|\left|M\right|\right|_{\textrm{op}}}/\sqrt{\operatorname{tr}(M_{+})}-Lipschitz continuous. For the case where the entries xx are bounded i.i.d. random variables, we can use the simple fact that the norm is convex in order to apply Corollary 4.10 in and Proposition 1.8 in in . For the case where the entries are normally distributed we can apply Theorem \@slowromancapv@.\@slowromancapi@ in . As a result, we can see that there exists a constant C4>0C_{4}>0 independent of MM, such that

The proof then follows straight forwardly following line by line the proof of Lemma A.2 in .

E.2 Proof of Lemma C.6

For any j≤m+1,j\leq m+1, α=(i1,i2)\alpha=(i_{1},i_{2}), let gj(α)g_{j}^{(\mathbf{\alpha})} denote the partial derivatives gj(α)(x,y)=∂∣α∣∂t1i1t2i2gj(t1,t2)∣x,yg_{j}^{(\mathbf{\alpha})}(x,y)=\frac{\partial^{|\alpha|}}{\partial t_{1}^{i_{1}}t_{2}^{i_{2}}}g_{j}(t_{1},t_{2})|_{x,y}. Define s=⌊2/β⌋s=\lfloor 2/\beta\rfloor First of all, note that due to Lemma C.1, for any δ,δ′>0\delta,\delta^{\prime}>0 and nn sufficiently large, for any Z∈EZ∣ZZ\in\mathcal{E}_{Z|\mathbf{Z}}

As a result, we can make use of Assumption C.1. We are heavily going to make use of this fact throughout the proof. The proof is separated into two steps where we first show 1. and then 2. using the expression for pp from the first step. Proof of the first statement We construct a polynomial p(Z)p(Z) using the power series expansion of kk from Assumption A.1 and in addition the Taylor series approximation of glg_{l} around the point (1,1)(1,1). For any nn sufficiently large, we can write

where ηl1,l2l,i∈Br(1,1)\eta^{l,i}_{l_{1},l_{2}}\in B_{r}(1,1) are points contained in the closed ball around the point (1,1)(1,1) with radius r2=(∥zi∥22−1)2+(∥Z∥22−1)2→0r^{2}=(\|z_{i}\|_{2}^{2}-1)^{2}+(\|Z\|_{2}^{2}-1)^{2}\to 0. Hence, using the fact that gig_{i} is s+1−is+1-i-times continuously differentiable, we can see that any ∥gl(l1,l2)(ηl1,l2l,i)∥\|g_{l}^{(l_{1},l_{2})}(\eta^{l,i}_{l_{1},l_{2}})\| is almost surely upper bounded by some constant.

Let vZv_{Z} be the vector defined in Equation (25) We define the polynomial pp as

Note that pp is a linear combination of the terms (Z⊤Z)p1(zi⊤Z)p2(Z^{\top}Z)^{p_{1}}(z_{i}^{\top}Z)^{p_{2}} with p1+p2≤sp_{1}+p_{2}\leq s, and hence a polynomial of ZZ of degree at most 2s2s. If glg_{l} are constant, i.e. the kernel is an inner product kernel, vZv_{Z} contains only the terms (zi⊤Z)l(z_{i}^{\top}Z)^{l} and hence p(Z)p(Z) is a polynomial of ZZ of degree at most ss.

In order to conclude the first step of the proof, we only need to show that both terms go to zero. First, we show that the term nmax⁡iB1i→0n\max\limits_{i}B^{i}_{1}\to 0. Recall that ∣gq(l1,l2)(ηl1,l2q,i)∣\left|g_{q}^{(l_{1},l_{2})}(\eta^{q,i}_{l_{1},l_{2}})\right| is upper bounded as n→∞n\to\infty independent of ii. Hence, we can apply Lemma C.1 which shows that for any integers q,l1q,l_{1} and l2l_{2} such that q+l1+l2=s+1q+l_{1}+l_{2}=s+1,

which holds true for any positive constant ϵ>0\epsilon>0. Finally, because s=⌊2/β⌋s=\lfloor 2/\beta\rfloor, (β/2)(s+1)>1(\beta/2)(s+1)>1, and hence

Furthermore, because rr is a continuous function and zi,Zz_{i},Z are contained in a closed neighborhood around of (1,1,0)(1,1,0), r(∥zi∥,∥Z∥,zi⊤Z)r(\|z_{i}\|,\|Z\|,z_{i}^{\top}Z) is upper bounded by some constant independent of ii as n→∞n\to\infty. Therefore, we also have that

where we have again used Lemma C.1. Hence, we can conclude the first step of the proof when observing that we have only assumed that Z∈EZ∣ZZ\in\mathcal{E}_{Z|\mathbf{Z}} and hence the convergence is uniformly. Proof of the second statement We can see from the definition of pp and the subsequent discussion that

We can decompose for nn sufficiently large,

where we have used Lemma C.1 in the second inequality. The first term vanishes trivially from the concentration inequality. Indeed, for nn sufficiently large, we have that

For the second term, note that we can see from the proof of Lemma E.1 that there exists some constant c>0c>0, such that for nn sufficiently large,

We can now apply integration by parts to show that

Hence, combining these terms, we get the desired result

E.3 Proof of Lemma C.7

As in , we separately analyze the off and on diagonal terms of KK. Let AA be the off-diagonal matrix of KK, with diagonal entries Ai,j=(1−δi,j)Ki,jA_{i,j}=(1-\delta_{i,j})K_{i,j} and let DD be the diagonal matrix of KK with entries Di,j=δi,jKi,jD_{i,j}=\delta_{i,j}K_{i,j}. We have

Similarily, decompose MM into its off-diagonal, MAM_{A}, and its diagonal MDM_{D}. We have

We begin with the first term. Note that MAM_{A} has off-diangoal entries (MA)i,j:=∑q=0m(zi⊤zj)qgq(∥zi∥22,∥zj∥22)(M_{A})_{i,j}:=\sum_{q=0}^{m}(z_{i}^{\top}z_{j})^{q}g_{q}(\|z_{i}\|_{2}^{2},\|z_{j}\|_{2}^{2}), and hence,

where we have the same argument as used in the proof of Lemma C.6 and the fact that the Assumptions A.1-A.3 imply Assumption C.1, as shown in Lemma E.2.

Because D−MDD-M_{D} is a diagonal matrix, for nn sufficiently large,

where we have used that by assumption gg is δL\delta_{L}-Lipschitz continuous on the restriction {(x,x,x)∣x∈[1−δL,1+δL]}⊂Ω\{(x,x,x)|x\in[1-\delta_{L},1+\delta_{L}]\}\subset\Omega for some δL>0\delta_{L}>0. Clearly T1→0T_{1}\to 0 due to Lemma C.1. Furthermore, by Assumption C.1, for any q≤mq\leq m, gqg_{q} is continuously differentiable and hence also Lipschitz continuous in a closed ball around (1,1)(1,1). Thus, T3→0T_{3}\to 0. Hence, it is only left to show that T2→0T_{2}\to 0, which is a consequence of the following claim. Claim: For any ϵ>0\epsilon>0 and any q>0q>0,

where cqc_{q} is a constant only depending on qq. Proof of the claim: In order to prove the claim, recall that due to Lemma C.1, for every q>0q>0

We prove the claim by induction. The case where q=1q=1 holds trivially with c1=1c_{1}=1. For q>1q>1,

Next, by induction, max⁡i ∣(xi⊤xi)j(tr⁡(Σd))j−1∣≤c1max⁡[n−β/2(log⁡(n))(1+ϵ)/2,n−jβ/2(log⁡(n))j((1+ϵ)/2)]\underset{i}{\max}~{}\left|\frac{(x_{i}^{\top}x_{i})^{j}}{(\operatorname{tr}(\Sigma_{d}))^{j}}-1\right|\leq c_{1}\max\left[n^{-\beta/2}(\log(n))^{(1+\epsilon)/2},n^{-j\beta/2}(\log(n))^{j((1+\epsilon)/2)}\right] for any j<qj<q. Furthermore, ∑j=1q(−1)j(qj)=−1\sum_{j=1}^{q}(-1)^{j}{q\choose j}=-1, which shows that

which completes the induction and thus the proof ∎

E.4 Proof of Lemma C.8

which holds for all t≥0t\geq 0. Hence, for t:=∥zi−zj∥2α⩾0t:=\|z_{i}-z_{j}\|_{2}^{\alpha}\geqslant 0, we can write

E.5 Proof of Lemma C.10

We start the proof with a discussion of existing results in the literature. As shown in Appendix E.1 in , the homogeneity of σ\sigma allows us to write

The function tσt_{\sigma} is continuous in $andsmoothinand smooth in(-1,1)$

Based on this discussion, we now prove the lemma. In a first step, we derive a closed form expression for Σ(i)(x,x)\Sigma^{(i)}(x,x). Based on the discussion above and particularly Equation (30) we can see that

The goal is now to show that whenever x,x′≠0x,x^{\prime}\neq 0, Σ(i)(x,x′)\Sigma^{(i)}(x,x^{\prime}) can be expressed as a sum of the form

with ηj,l(i)≥0\eta_{j,l}^{(i)}\geq 0. We prove by induction. In a first step, note that the case where i=0i=0 holds trivially true. Next, assume that Equation (32) holds true for Σ(i−1)\Sigma^{(i-1)}. Due to the above discussion, tσt_{\sigma} can be expressed as a Taylor series around with positive coefficients ai2a_{i}^{2}. Thus,

Next, note that any of the properties 1-4 from the above discussion also hold true for tσ˙t_{\dot{\sigma}}. Therefore, we can use exactly the same argument for Σ˙(i)\dot{\Sigma}^{(i)} to show that for any x,x′≠0x,x^{\prime}\neq 0,

with η˙j,l(i+1)≥0\dot{\eta}_{j,l}^{(i+1)}\geq 0. Finally, we can conclude the proof because

and when using the same argument as used in the induction step above to show that the resulting multi-sum converges absolutely. ∎

E.6 Additional lemmas

Any kernel which satisfies Assumption A.1 and A.3 also satisfies Assumption C.1.

The only point which does not follow immediately is to show that rr is a continuous function. For this, write gg as a function of the variables x,y,zx,y,z, i.e.

For every x,yx,y, define the function gx,y(z)=g(x,y,z)g_{x,y}(z)=g(x,y,z). Due to the series expansion, we can make use of the theory on the Taylor expansion which implies that for any x,yx,y, gx,yg_{x,y} is a smooth function in the interior of N(δ)N(\delta) (using the definition from Assumption A.1). Hence, we can conclude that there exists a function rx,yr_{x,y} such that

In particular, the smoothness of gx,y(z)g_{x,y}(z) implies that rx,yr_{x,y} is continuous in the interior of {z:(x,y,z)∈N(δ)}\{z:(x,y,z)\in N(\delta)\}. Next, define r(x,y,z)=rx,y(z)r(x,y,z)=r_{x,y}(z) and note that that the continuity of gg implies that rr is continuous everywhere except for the plane z=0z=0. Finally, because rx,y(0)r_{x,y}(0) exists point wise, we conclude that rr exists and is a continuous function in the interior of N(δ)N(\delta). ∎

Any RBF kernel k(x,x′)=h(∥x−x′∥22)k(x,x^{\prime})=h(\|x-x^{\prime}\|_{2}^{2}) with hh locally analytic around 22 satisfies Assumption C.1.

Because by assumption hh has a local Taylor series around 22, we can write

is a sum of non negative summands. Hence, we get that the sum

converges absolutely. Thus, we can arbitrarily reorder the summands: