Just Interpolate: Kernel "Ridgeless" Regression Can Generalize

Introduction

According to conventional wisdom, explicit regularization should be added to the least-squares objective when the Hilbert space $\mathcal{H}$ is high- or infinite-dimensional (Golub et al., 1979; Wahba, 1990; Smola and Schölkopf, 1998; Shawe-Taylor and Cristianini, 2004; Evgeniou et al., 2000; De Vito et al., 2005; Alvarez et al., 2012):

The regularization term is introduced to avoid “overfitting” since kernels provide enough flexibility to fit training data exactly (i.e. interpolate it). From the theoretical point of view, the regularization parameter $\lambda$ is a knob for balancing bias and variance, and should be chosen judiciously. Yet, as noted by a number of researchers in the last few years,In particular, we thank M. Belkin, B. Recht, L. Rosasco, and N. Srebro for highlighting this phenomenon. the best out-of-sample performance, empirically, is often attained by setting the regularization parameter to zero and finding the minimum-norm solution among those that interpolate the training data. The mechanism for good out-of-sample performance of this interpolation method has been largely unclear (Zhang et al., 2016; Belkin et al., 2018b).

As a concrete motivating example, consider the prediction performance of Kernel Ridge Regression for various valuesWe take $\lambda\in\{0,0.01,0.02,0.04,0.08,0.16,0.32,0.64,1.28\}$. of the regularization parameter $\lambda$ on subsets of the MNIST dataset. For virtually all pairs of digits, the best out-of-sample mean squared error is achieved at $\lambda=0$. Contrary to the standard bias-variance-tradeoffs picture we have in mind, the test error is monotonically decreasing as we decrease $\lambda$ (see Figure 1 and further details in Section 6).

We isolate what appears to be a new phenomenon of implicit regularization for interpolated minimum-norm solutions in Kernel “Ridgeless” Regression. This regularization is due to the curvature of the kernel function and “kicks in” only for high-dimensional data and for “favorable” data geometry. We provide out-of-sample statistical guarantees in terms of spectral decay of the empirical kernel matrix and the empirical covariance matrix, under additional technical assumptions.

Our analysis rests on the recent work in random matrix theory. In particular, we use a suitable adaptation of the argument of (El Karoui, 2010) who showed that high-dimensional random kernel matrices can be approximated in spectral norm by linear kernel matrices plus a scaled identity. While the message of (El Karoui, 2010) is often taken as “kernels do not help in high dimensions,” we show that such a random matrix analysis helps in explaining the good performance of interpolation in Kernel “Ridgeless” Regression.

Grace Wahba (Wahba, 1990) pioneered the study of nonparametric regression in reproducing kernel Hilbert spaces (RKHS) from the computational and statistical perspectives. One of the key aspects in that work is the role of the decay of eigenvalues of the kernel (at the population level) in rates of convergence. The analysis relies on explicit regularization (ridge parameter $\lambda$) for the bias-variance trade-off. The parameter is either chosen to reflect the knowledge of the spectral decay at the population level (De Vito et al., 2005) (typically unknown to statistician), or by the means of cross-validation (Golub et al., 1979). Interestingly, the explicit formula of Kernel Ridge Regression has been introduced as “kriging” in the literature before, and was widely used in Bayesian statistics (Cressie, 1990; Wahba, 1990).

In the learning theory community, Kernel Ridge Regression is known as a special case of Support Vector Regression (Vapnik, 1998; Shawe-Taylor and Cristianini, 2004; Vovk, 2013). Notions like metric entropy (Cucker and Smale, 2002) or “effective dimension” (Caponnetto and De Vito, 2007) were employed to analyze the guarantees on the excess loss of Kernel Ridge Regression, even when the model is misspecified. We refer the readers to Györfi et al. (2006) for more details. Again, the analysis leans crucially on the explicit regularization, as given by a careful choice of $\lambda$, for the model complexity and approximation trade-off, and mostly focusing on the fixed dimension and large sample size setting. However, to the best of our knowledge, the literature stays relatively quiet in terms of what happens to the minimum norm interpolation rules, i.e., $\lambda=0$. As pointed out by (Belkin et al., 2018b, a), the existing bounds in nonparametric statistics and learning theory do not apply to interpolated solution either in the regression or the classification setting. In this paper, we aim to answer when and why interpolation in RKHS works, as a starting point for explaining the good empirical performance of interpolation using kernels in practice (Zhang et al., 2016; Belkin et al., 2018b).

Preliminaries

In this paper we study the case when $K(X,X)$ is full rank, taking (2.2) as the starting point. For this interpolating estimator, we provide high-probability (with respect to a draw of $X$) upper bounds on the integrated squared risk of the form

Here the expectation is over $\mathbf{x}\sim\mu$ and $Y|X$, and $\phi_{n,d}$ is a data-dependent upper bound. We remark that upper bounds of the form (2.3) also imply prediction loss bounds for excess square loss with respect to the class $\mathcal{H}$, as $\mathbf{E}(\widehat{f}(\mathbf{x})-f_{*}(\mathbf{x}))^{2}=\mathbf{E}(\widehat{f}(\mathbf{x})-\mathbf{y})^{2}-\mathbf{E}(f_{*}(\mathbf{x})-\mathbf{y})^{2}$.

2 Notation and Background on RKHS

Let us introduce the integral operator $\mathcal{T}_{\mu}\mathrel{\mathop{\mathchar 58\relax}}L^{2}_{\mu}\rightarrow L^{2}_{\mu}$ with respect to the marginal measure $\mu(x)$:

and denote the set of eigenfunctions of this integral operator by $e(x)=\{e_{1}(x),e_{2}(x),\ldots,e_{p}(x)\}$, where $p$ could be $\infty$. We have that

Denote $T={\rm diag}(t_{1},\ldots,t_{p})$ as the collection of non-negative eigenvalues. Adopting the spectral notation,

Via this spectral characterization, the interpolation estimator (2.1) takes the following form

and the associated eigenvalues to be $\lambda_{j}(\widehat{\mathcal{T}})$, indexed by $j$. The eigenvalues are the same as those of $\frac{1}{n}K(X,X)$. It is sometimes convenient to express $\widehat{\mathcal{T}}$ as the linear operator under the basis of eigenfunctions, in the following matrix sense

We write $\mathbf{E}_{\mu}[\cdot]$ to denote the expectation with respect to the marginal $\mathbf{x}\sim\mu$. Furthermore, we denote by

the squared $L^{2}$ norm with respect to the marginal distribution. The expectation $\mathbf{E}_{Y|X}[\cdot]$ denotes the expectation over $y_{1},\ldots,y_{n}$ conditionally on $x_{1},\ldots,x_{n}$.

Main Result

High dimensionality: there exists universal constants $c,C\in(0,\infty)$ such that $c\leq d/n\leq C$. Denote by $\Sigma_{d}=\mathbf{E}_{\mu}[x_{i}x_{i}^{*}]$ the covariance matrix, assume that the operator norm $\|\Sigma_{d}\|_{\rm op}\leq 1$.

Noise condition: there exists a $\sigma>0$ such that $\mathbf{E}[(f_{*}(\mathbf{x})-\mathbf{y})^{2}|\mathbf{x}=x]\leq\sigma^{2}$ for all $x\in\Omega$.

Non-linear kernel: for any $x\in\Omega$, $K(x,x)\leq M$. Furthermore, we consider the inner-product kernels of the form

While we state the main theorem for inner product kernels, the results follow under suitable modificationsWe refer the readers to El Karoui (2010) for explicit extensions to RBF kernels. for Radial Basis Function (RBF) kernels of the form

We postpone the discussion of the assumptions until after the statement of the main theorem.

Let us first define the following quantities related to curvature of $h$:

Under the assumptions (A.1)-(A.4) and for $d$ large enough, with probability at least $1-2\delta-d^{-2}$ (with respect to a draw of design matrix $X$), the interpolation estimator (2.2) satisfies

Here the remainder term $\epsilon(n,d)=O(d^{-\frac{m}{8+m}}\log^{4.1}d)+O(n^{-\frac{1}{2}}\log^{0.5}(n/\delta))$.

A few remarks are in order. First, the upper bound is data-dependent and can serve as a certificate (assuming that an upper bound on $\sigma^{2},\|f_{*}\|_{\mathcal{H}}^{2}$ can be guessed) that interpolation will succeed. The bound also suggests the regimes when the interpolation method should work. The two terms in the estimate of Theorem 1 represent upper bounds on the variance and bias of the interpolation estimator, respectively. Unlike the explicit regularization analysis (e.g. (Caponnetto and De Vito, 2007)), the two terms are not controlled by a tunable parameter $\lambda$. Rather, the choice of the non-linear kernel $K$ itself leads to an implicit control of the two terms through curvature of the kernel function, favorable properties of the data, and high dimensionality. We remark that for the linear kernel ($h(a)=a$), we have $\gamma=0$, and the bound on the variance term can become very large in the presence of small eigenvalues. In contrast, curvature of $h$ introduces regularization through a non-zero value of $\gamma$. We also remark that the bound “kicks in” in the high-dimensional regime: the error term decays with both $d$ and $n$.

We left the upper bound of Theorem 1 in a data-dependent form for two reasons. First, an explicit dependence on the data tells us whether interpolation can be statistically sound on the given dataset. Second, for general spectral decay, current random matrix theory falls short of characterizing the spectral density non-asymptotically except for special cases (Bose et al., 2003; El Karoui, 2010).

The assumption in (A.1) that $c\leq d/n\leq C$ emphasizes that we work in a high-dimensional regime where $d$ scales on the order of $n$. This assumption is used in the proof of (El Karoui, 2010), and the particular dependence on $c,C$ can be traced in that work if desired. Rather than doing so, we “folded” these constants into mild additional power of $\log d$. The same goes for the assumption on the scaling of the trace of the population covariance matrix.

The assumption in (A.2) that $Z_{i}(k)$ are i.i.d. across $k=1,\ldots,d$ is a strong assumption that is required to ensure the favorable high-dimensional effect. Relaxing this assumption is left for future work.

The existence of $(8+m)$-moments for $|z_{i}(k)|$ is enough to ensure $|z_{i}(k)|\leq C\cdot d^{\frac{2}{8+m}}$ for $1\leq i\leq n,1\leq k\leq d$ almost surely (see, Lemma 2.2 in Yin et al. (1988)). Remark that the assumption of existence of $(8+m)$-moments in (A.2) is relatively weak. In particular, for bounded or subgaussian variables, $m=\infty$ and the error term $\epsilon(n,d)$ scales as $d^{-1}+n^{-1/2}$, up to log factors. See Lemma B.1 for an explicit calculation in the Gaussian case.

Finally, as already mentioned, the main result is stated for the inner product kernel, but can be extended to the RBF kernel using an adaptation of the analysis in (El Karoui, 2010).

Behavior of the Data-dependent Bound

We can further bound the variance and the bias, with the choice $k=0$, as

We first illustrate numerically the bias-variance trade-off by varying the geometric properties of the data in terms of the population spectral decay of $\mathbf{x}$. We shall parametrize the eigenvalues of the covariance, for $0<\kappa<\infty$, as

The parameter $\kappa$ controls approximate “low-rankness” of the data: the closer $\kappa$ is to , the faster does the spectrum of the data decay. This is illustrated in the top row of Figure 6 on page 6. By letting $\kappa\rightarrow 0$, ${\bf r}$ can be arbitrary small, as

We will focus on three cases, $\kappa\in\{e^{-1},e^{0},e^{1}\}$, for the decay parameter, and values $d=100$, $n\in\{500,2000\}$. The data-dependent upper bounds on $\mathbf{V}$ and $\mathbf{B}$ are summarized in Table 1. More detailed plots are postponed to Figure 6 (in this figure, we plot the ordered eigenvalues and the spectral density for both the population and empirical covariances). Table 1 shows that as $\kappa$ increases (a slower spectral decay), the implicit regularization parameter becomes larger, resulting in a decreasing variance and an increasing bias.

We also perform simulations to demonstrate the trade-off between bias and variance in the generalization error. The result is shown in Figure 2. For each choice of $(n,d)$ pair, we vary the spectral decay of the kernel by changing gradually $\kappa\in[e^{-2},e^{2}]$, and plot the generalization error on the log scale. We postpone the experiment details to Section 6.2, but point out important phenomenona observed in Figures 2-3: (1) an extremely fast spectral decay (small $\kappa$) will generate insufficient implicit regularization that would hurt the generalization performance due to a large variance term; (2) a very slow spectral decay (large $\kappa$) will result in a large bias, which can also hurt the generalization performance; (3) certain favorable spectral decay achieves the best trade-off, resulting in the best generalization error.

We now theoretically demonstrate scalings within the $n>d$ regime when both $\bf V$ and $\bf B$ vanish. For simplicity, we consider Gaussian $\mathbf{x}$.

Consider general eigenvalue decay with $\|\Sigma_{d}\|_{\rm op}\leq 1$. Then with high probability,

To illustrate the behavior of the estimates in Corollary 4.1, consider the following assumptions on the population covariance matrix:

Let $\Sigma_{d}={\rm diag}(1,\ldots,1,0,\ldots,0)$ with $\epsilon d$ ones, $\epsilon\in(0,1)$. In this case ${\bf r}=\epsilon^{2}$, and $\lambda_{j}(XX^{*}/n)\geq(1-\sqrt{\epsilon d/n})^{2}$ with high probability by standard results in random matrix theory. Then

Therefore, as $\epsilon\rightarrow 0$, both terms vanish for $n>d$.

Let $\Sigma_{d}={\rm diag}(1,\epsilon,\ldots,\epsilon)$ for small $\epsilon>0$. In this case, ${\bf r}=\epsilon^{2}$ and $\lambda_{j}(XX^{*}/n)\geq\epsilon(1-\sqrt{d/n})^{2}$ with high probability. Then

For instance, for $\epsilon\asymp(d/n)^{1/2}$, both terms vanish for $n\gg d$.

Consider $\lambda_{j}(\Sigma_{d})=j^{-\alpha}$ for $0<\alpha<1$. Then ${\bf r}\asymp d^{-2\alpha}$. One can bound w.h.p. (see (B.4))

Balancing the two terms, one obtains a nonparametric upper bound $n^{-\frac{\alpha}{2\alpha+1}}.$ A similar analysis can be carried out for $\alpha\geq 1$.

Case d>n𝑑𝑛d>n

In this case, we can further bound the variance and the bias, with the choice $k=0$, as

We first numerically illustrate the trade-off between the variance and the bias upper bounds. We consider three cases $\kappa\in\{e^{-1},e^{0},e^{1}\}$, and $d=2000$, $n\in\{400,100\}$. As before, we find a trade-off between $\mathbf{V}$ and $\mathbf{B}$ with varying $\kappa$; the results are summarized in Table 2. Additionally, Figure 7 provides a plot of the ordered eigenvalues, as well as spectral density for both the population and empirical covariances. As one can see, for a general eigenvalue decay, the spectral density of the population and the empirical covariance can be quite distinct. We again plot the generalization error in Figure 3 as a function of $\kappa$.

We now theoretically showcase an example in the $d\gg n$ regime where both $\bf V$ and $\bf B$ vanish. Again consider $\mathbf{x}$ being Gaussian for simplicity.

The variance bound follows from the fact that $\frac{t}{({\bf r}+t)^{2}}\leq\frac{1}{4{\bf r}}$ for all $t$.

Proofs

To prove Theorem 1, we decompose the mean square error into the bias and variance terms (Lemma 5.1), and provide data-dependent bound for each (Sections 5.2 and 5.3).

The following is a standard bias-variance decomposition for an estimator. We remark that it is an equality, and both terms have to be small to ensure the desired convergence.

The following decomposition for the interpolation estimator (2.2) holds

Recall the closed form solution of the interpolation estimator:

Define $E=Y-\mathbf{E}[Y|X]=Y-f_{*}(X)$. Since $\mathbf{E}_{Y|X}E=0$, we have

2 Variance

In this section, we provide upper estimates on the variance part $\mathbf{V}$ in (5.2).

Let $\delta\in(0,1)$. Under the assumptions (A.1)-(A.4), with probability at least $1-\delta-d^{-2}$ with respect to a draw of $X$,

for $\theta=\frac{1}{2}-\frac{2}{8+m}$ and for $d$ large enough.

Let us discuss the first term in Eq. (5.4) and its role in implicit regularization induced by the curvature of the kernel, eigenvalue decay, and high dimensionality. In practice, the data matrix $X$ is typically centered, so $1^{*}X=0$. Therefore the first term is effectively

This formula explains the effect of implicit regularization, and captures the “effective rank” of the training data $X$. We would like to emphasize that this measure of complexity is distinct from the classical notion of effective rank for regularized kernel regression (Caponnetto and De Vito, 2007), where the “effective rank” takes the form $\sum_{j}g_{r}(t_{j})$ with $g_{r}(t)=t/(r+t)$, with $t_{j}$ is the eigenvalue of the population integral operator $\mathcal{T}$.

From the definition of $\mathbf{V}$ and $E[Y|X]=f_{*}(X)$,

Due to the fact that $\mathbf{E}_{Y|X}\left[(Y_{i}-f_{*}(X_{i}))(Y_{j}-f_{*}(X_{j}))\right]=0$ for $i\neq j$, and $\mathbf{E}_{Y|X}\left[(Y_{i}-f_{*}(X_{i}))^{2}\right]\leq\sigma^{2}$, we have that $\|\mathbf{E}_{Y|X}\left[(Y-f_{*}(X))(Y-f_{*}(X))^{*}\right]\|\leq\sigma^{2}$ and thus

Let us introduce two quantities for the ease of derivation. For $\alpha,\beta,\gamma$ defined in (3), let

and $K^{\rm lin}(x,X)$ being the transpose of $K^{\rm lin}(X,x)$. By Proposition A.2, with probability at least $1-\delta-d^{-2}$, for $\theta=\frac{1}{2}-\frac{2}{8+m}$ the following holds

As a direct consequence, one can see that

provided $d$ is large enough, in the sense that

By Lemma B.2 (for Gaussian case, Lemma B.1),

where the the third inequality relies on (5.9) and (5.7), and the fourth inequality follows from (5.8).

We conclude that with probability at least $1-\delta-d^{-2}$,

3 Bias

Let $\delta\in(0,1)$. The bias, under the only assumptions that $K(x,x)\leq M$ for $x\in\Omega$, and $X_{i}$’s are i.i.d. random vectors, is upper bounded as

In this proof, when there is no confusion, we use $f(x)=\sum_{i=1}^{p}e_{i}(x)f_{i}$ where $f_{i}$ denotes the coefficients of $f$ under the basis $e_{i}(x)$. Adopting this notation, we can write $f(x)=e(x)^{*}f$ where $f=[f_{1},f_{2},\ldots,f_{p}]^{T}$ also denotes a possibly infinite vector. For the bias, it is easier to work in the frequency domain using the spectral decomposition. Recalling the spectral characterization in the preliminary section,

Here we use the fact that $T^{-1/2}f_{*}=\sum_{i}t_{i}^{-1/2}f_{*,i}e_{i}$ and $\|T^{-1/2}f_{*}\|^{2}=\sum_{i}f_{*,i}^{2}/t_{i}=\|f_{*}\|_{\mathcal{H}}^{2}$. Next, recall the empirical Kernel operator with its spectral decomposition $\widehat{T}=\widehat{U}\widehat{\Lambda}\widehat{U}^{*}$, with $\widehat{\Lambda}_{jj}=\frac{1}{n}\lambda_{j}\left(K(X,X)\right)$. Denote the top $k$ columns of $\widehat{U}$ to be $\widehat{U}_{k}$, and $P_{\widehat{U}_{k}}^{\perp}$ to be projection to the eigenspace orthogonal to $\widehat{U}_{k}$. By observing that $T^{1/2}e(X)(e(X)^{*}Te(X))^{-1}e(X)^{*}T^{1/2}$ is a projection matrix, it is clear that for all $k\leq n$,

We continue the study of the last quantity using techniques inspired by Shawe-Taylor and Cristianini (2004). Denote the function $g$ indexed by any rank-$k$ projection $U_{k}$ as

Clearly, $\|U_{k}U_{k}^{T}\|_{F}=\sqrt{k}$. Define the function class

It is clear that $g_{\widehat{U}_{k}}\in\mathcal{G}_{k}$. Observe that $g_{\widehat{U}_{k}}$ is a random function that depends on the data $X$, and we will bound the bias term using the empirical process theory. It is straightforward to verify that

Using symmetrization Lemma B.4 with $M=\sup_{x\in\Omega}K(x,x)$, with probability at least $1-2\delta$,

by the Pythagorean theorem. Since $\epsilon_{i}$’s are symmetric and zero-mean and $\left\|T^{1/2}e(x_{i})\right\|^{2}$ does not depend on $U_{k}$, the last expression is equal to

We further bound the Rademacher complexity of the set $\mathcal{G}_{k}$

by the Cauchy-Schwarz inequality and the fact that $\|U_{k}U_{k}^{T}\|_{F}\leq\sqrt{k}$. The last expression is can be further evaluated by the independence of $\epsilon_{i}$’s

Therefore, for all $k\leq n$, with probability at least $1-2n\delta$,

Let us compare the bounds obtained in this paper to those one can obtain for classification with a margin. For classification, Thm. 21 in Bartlett and Mendelson (2002) shows that the misclassification error is upper bounded with probability at least $1-\delta$ as

where $\phi_{\gamma}(t)\mathrel{\mathop{\mathchar 58\relax}}=\max(0,1-t/\gamma)\wedge 1$ is the margin loss surrogate for the indicator loss ${\bf 1}(t<0)$. By tuning the margin $\gamma$, one obtains a family of upper bounds.

Now consider the noiseless regression scenario (i.e. $\sigma=0$ in (A.1)). In this case, the variance contribution to the risk is zero, and

where $P_{k}$ is the best-rank $k$ projection (based on $X$) and $P_{k}^{\perp}$ denotes its orthogonal projection. By tuning the parameter $k$ (similar as the $1/\gamma$ in classification), one can balance the RHS to obtain the optimal trade-off.

However, classification is easier than regression in the following sense: $\widehat{f}$ can present a non-vanishing bias in estimating $f_{*}$, but as long as the bias is below the empirical margin level, it plays no effect in the margin loss $\phi_{\gamma}(\cdot)$. In fact, for classification, under certain conditions, one can prove exponential convergence for the generalization error (Koltchinskii and Beznosova, 2005).

Experiments

In this section we provide full details of the experiments on MNIST (LeCun et al., 2010). Our first experiment considers the following problem: for each pair of distinct digits $(i,j)$, $i,j\in\{0,1,\ldots,9\}$, label one digit as $1$ and the other as $-1$, then fit the Kernel Ridge Regression with Gaussian kernel $k(x,x^{\prime})=\exp(-\|x-x^{\prime}\|^{2}/d)$, where $d=784$ is the dimension as analyzed in our theory (also the default choice in Scikit-learn package (Pedregosa et al., 2011)). For each of the $\binom{10}{2}=45$ pairs of experiments, we chose $\lambda=0$ (no regularization, interpolation estimator), $\lambda=0.1$ and $\lambda=1$. We evaluated the performance on the out-of-sample test dataset, with the error metric

Remarkably, among all 45 experiments, no-regularization performs the best. We refer to the table in Section C for a complete list of numerical results. For each experiment, the sample size is roughly $n\approx 10000$.

The second experiment is to perform the similar task on a finer grid of regularization parameter $\lambda\in\{0,0.01,0.02,0.04,0.08,0.16,0.32,0.64,1.28\}$. Again, in all but one pair, the interpolation estimator performs the best in out-of-sample prediction. We refer to Figure 4 for details.

To conclude this experiment, we plot the eigenvalue decay of the empirical kernel matrix and the sample covariance matrix for the 5 experiments shown in the introduction. The two plots are shown in Figure 5. Both plots exhibit a fast decay of eigenvalues, supporting the theoretical finding that interpolation performs well on a test set in such situations.

On the other hand, it is easy to construct examples where the eigenvalues do not decay and interpolation performs poorly. This is the case, for instance, if $X_{i}$ are i.i.d. from spherical Gaussian. One can show that in the high-dimensional regime, the variance term itself (and not just the upper bound on it) is large. Since the bias-variance decomposition is an equality, it is not possible to establish good $L^{2}_{\mu}$ convergence.

2 A Synthetic Example

In this section we provide the details of the synthetic experiments mentioned in Section 4 for Tables 1-2 and Figures 2-3. We choose the RBF kernel as the non-linearity with $h(t)=\exp(-t)$. Again, we consider a family of eigenvalue decays for the covariance matrix parametrized by $\kappa$, with the small $\kappa$ describing fast spectral decay

We set a target non-linear function $f_{*}$ in the RKHS with kernel $K(x,x^{\prime})=h\left(\|x-x^{\prime}\|^{2}/d\right)$ as

For each parameter triplet $(n,d,\kappa)$, we generate data in the following way

for $1\leq i\leq n$ where $\epsilon_{i}\sim N(0,\sigma^{2})$ is independent noise, with $\sigma=0.1$ (Figures 2-3) and $\sigma=0.5$ (Figures 8). Figures 6-7 contrasts the difference between the population and empirical eigenvalues for various parameter triplets $(n,d,\kappa)$.

We now explain Figures 2-3, which illustrate the true generalization error in this synthetic example, by varying the spectral decay $\kappa$, for a particular case of high dimensionality ratio $d/n$. Here we plot the out-of-sample test error for the interpolated min-norm estimator $\widehat{f}$ on fresh new test data $(x_{t},y_{t})$ from the same data generating process, with the error metric

The error plots are shown in Figure 2 (for $n>d$) and 3 (for $d>n$), and Figure 8 for the high noise case. On the x-axis, we plot the $\log(\kappa)$, and on the y-axis the $\log(\text{error})$. Each curve corresponds to the generalization error behavior (and the bias and variance trade-off) as we vary spectral decay from fast to slow (as $\kappa$ increases) for a particular choice of $d/n$ or $n/d$ ratio. Clearly, for a general pair of high dimensionality ratio $d/n$, there is a “sweet spot” of $\kappa$ (favorable geometric structure) such that the trade-off is optimized.

Further Discussion

This paper is motivated by the work of Belkin et al. (2018b) and Zhang et al. (2016), who, among others, observed the good out-of-sample performance of interpolating rules. This paper continues the line of work in (Belkin et al., 2018a, c; Belkin, 2018) on understanding theoretical mechanisms for the good out-of-sample performance of interpolation. We leave further investigations on the connection between kernel ridgeless regression and two-layer neural networks as a future work (Dou and Liang, 2019).

From an algorithmic point of view, the minimum-norm interpolating solution can be found either by inverting the kernel matrix, or by performing gradient descent on the least-squares objective (starting from ). Our analysis can then be viewed in the light of recent work on implicit regularization of optimization procedures (Yao et al., 2007; Neyshabur et al., 2014; Gunasekar et al., 2017; Li et al., 2017).

The paper also highlights a novel type of implicit regularization. In addition, we discover that once we parametrize the geometric properties — the spectral decay — we discover the familiar picture of the bias-variance trade-off, controlled by the implicit regularization that adapts to the favorable geometric property of the data. Moreover, if one explicitly parametrizes the choice of the kernel by, say, the bandwidth, we are likely to see the familiar picture of the bias-variance trade-off, despite the fact that the estimator is always interpolating. Whether one can achieve optimal rates of estimation (under appropriate assumptions) for the right choice of the bandwidth appears to be an interesting and difficult statistical question. Another open question is whether one can characterize situations when the interpolating minimum-norm solution is dominating the regularized solution in terms of expected performance.

References

Appendix A Propositions

We first borrow a technical result for concentration of quadratic forms under a mild moment condition.

The proof follows almost exactly as in Lemma A.3 (El Karoui, 2010). The only point of clarification is that one can assert

and thus for $d$ large enough, say $c(\log d)^{1.02}\geq 4\log d+\log 2/c^{2}$, we have that $2n^{2}\exp(-c(\log d)^{1.02})\leq d^{-2}$. ∎

The following proposition is a non-asymptotic adaptation of Theorem 2.1 in (El Karoui, 2010). Our contribution here is only to carefully spell out the terms and emphasize that the error rate can be very slow (this is why (El Karoui, 2010) only provides a convergence in probability result).

Under the assumptions (A.1), (A.2), and (A.4), for $\theta=\frac{1}{2}-\frac{2}{8+m}$, with probability at least $1-\delta-d^{-2}$,

for $d$ large enough and $\delta$ small enough.

In (El Karoui, 2010), the approximation error can be decomposed into first-order term $E_{1}$ (diagonal approximation), second-order off-diagonal term $E_{2}$, and third-order off-diagonal approximation $E_{3}$,

where $\|E_{1}\|\leq\frac{1}{2}d^{-\theta}\log^{0.51}d$, $\|E_{3}\|\leq\frac{1}{2}d^{-\theta}\log^{0.51}d$ with probability at least $1-d^{-2}$. However, for $E_{2}$ only convergence in probability is obtained. On Page 19 in (El Karoui, 2010), the last line reads

Appendix B Lemmas and Corollaries

Under the assumptions (A.1), (A.4), and that $x_{i}\sim\mathcal{N}(0,\Sigma_{d})$ i.i.d. Then with probability at least $1-d^{-2}$ with respect to a draw of $X$,

Start with entry-wise Taylor expansion for the smooth kernel,

Conditionally on $x_{j},1\leq j\leq n$, with probability at least $1-2\exp(-t^{2}/2)$ on $\mathbf{x}$ drawn from $\mathcal{N}(0,\Sigma)$,

Using standard $\chi^{2}$ concentration bound, we know that with probability at least $1-d^{-2}$ on $X$

Therefore with probability at least $1-2\exp(-t^{2}/2)$ on $\mathbf{x}\sim\mu$, conditionally on $x_{j},1\leq j\leq n$, we have

Define $Z=\|K(\mathbf{x},X)-K^{\rm lin}(\mathbf{x},X)\|$. The above says that, conditioned on $X$, $P\left(Z\geq C\cdot d^{-1/2}(s+\log^{1.02}d)\right)\leq 2\exp(-s/2)$ for all $s>0$. Therefore, by defining change of variables $z=C\cdot d^{-1/2}(s+\log^{1.02}d)$,

with probability at least $1-d^{-2}$ on $X$, for $d$ large enough. ∎

Under the assumptions (A.1), (A.2), and (A.4), for $\theta=\frac{1}{2}-\frac{2}{8+m}$, we have with probability at least $1-d^{-2}$ with respect to the draw of $X$, for $d$ large enough,

We start with entry-wise Taylor expansion for the smooth kernel,

Conditionally on $X_{j},1\leq j\leq n$, by Bernstein’s inequality (Boucheron et al., 2013, p. 38), with probability at least $1-\exp(-t)$ on $\mathbf{x}$, for all $j\in[n]$

Here the second line uses the fact that $\max_{k}\left|[\Sigma^{1/2}x_{j}](k)\cdot[\Sigma^{-1/2}\mathbf{x}](k)\right|\leq\|\Sigma^{1/2}x_{j}\|_{\infty}d^{\frac{2}{8+m}}$ due to the assumption (A.2) for each entry $[\Sigma^{-1/2}\mathbf{x}](k)\leq Cd^{\frac{2}{8+m}}$. Applying Proposition A.1 with the matrix taken to be identity, for all $j$, with probability at least $1-d^{-2}$ on $X$

Therefore with probability at least $1-\exp(-t)$ with respect to $\mathbf{x}\sim\mu$, conditionally on $x_{j},1\leq j\leq n$

Define $Z=\|K(\mathbf{x},X)-K^{\rm lin}(\mathbf{x},X)\|$. The above says that, conditioned on $X$, $P\left(Z\geq C\cdot d^{-2\theta+1/2}(s^{2}+\log^{2.04}d)\right)\leq 2\exp(-s)$ for all $s>0$. Therefore, by change of variables $z=C\cdot d^{-2\theta+1/2}(s^{2}+\log^{2.04}d)$, the expectation satisfies

with probability at least $1-d^{-2}$ on $X$, for $d$ large enough.

If $A\succeq B$, then $\|A^{-1}C\|_{\rm op}\leq\|B^{-1}C\|_{\rm op}$. Since $\Sigma^{-1}\succeq I$, it holds that

For the variance part, with $Z=X\Sigma_{d}^{-1/2}$, we have

where the last step uses (B.4). Therefore, using standard random matrix theory, one can further upper bound the above equation by

by the same argument as in the proof of Corollary 4.1. ∎

where $\mathbf{E}_{\epsilon}$ denotes the conditional expectation with respect to i.i.d. Rademacher random variables $\epsilon_{1},\ldots,\epsilon_{n}$.

The proof is a standard exercise using McDiarmid’s inequality and symmetrization. We include here for completeness. See (Mendelson, 2003, Theorem 2.21, 2.23 and their corollaries). ∎

Appendix C MNIST Result

Here the error is in percentage, so 2.921 corresponds to an error 2.921%.