Kernel interpolation

The risk splits exactly into variance and bias. Write $\hat f-f_*$ as the interpolant of the noise plus the projection residual $(P_X-I)f_*$, with $P_X$ the projection of $\mathcal H$ onto $\operatorname{span}\{K(\cdot,x_i)\}$. The cross term has conditional mean zero, so the risk is exactly $V+B$, and $V\le\sigma^2\,\mathbb E_x\|K(X,X)^{-1}K(X,x)\|_2^2$: the variance depends only on the training matrix and the vector for a fresh test point, and the next two steps linearize each.

The training kernel matrix is a linear kernel plus an implicit ridge. This is El Karoui's linearization, made non-asymptotic. In high dimension normalized inner products are nearly deterministic: $|x_i^\top x_j/d|\lesssim d^{-\theta}$ for $i\ne j$ and $\|x_i\|^2/d=\tau+O(d^{-\theta})$, with $\theta=\tfrac12-\tfrac2{8+m}$ reflecting that only $8+m$ moments are assumed. Taylor-expand $h$ entrywise. Off the diagonal the argument is small, so the linear term gives $\beta XX^\top/d$ and the quadratic term, whose matrix of squares is close to its mean, gives $\alpha11^\top$. On the diagonal the argument is $\tau$, which is not small, so $K_{ii}\approx h(\tau)$ while the off-diagonal formula would predict $h(0)+h'(0)\tau$, and the gap is $\gamma I$. The centered quadratic part is controlled only through a fourth-moment trace bound, so Markov gives $\|K(X,X)-K^{lin}(X,X)\|\le\eta:=d^{-\theta}(\delta^{-1/2}+\log^{0.51}d)$, the one place the confidence level enters polynomially. Since $K^{lin}\succeq\gamma I$, Weyl gives $\|K^{-1}\|\le2/\gamma$ and $\|K^{-1}K^{lin}\|\le2$ once $\eta\le\gamma/2$. Without curvature there is no such control: a linear kernel has $\gamma=0$, and the bound on $K^{-1}$ fails.

The test vector is linear and carries no ridge. For a fresh $x$ every argument $x^\top x_j/d$ is small, including the one that would have been a diagonal entry, so $K(x,x_j)=h(0)+\beta x^\top x_j/d+O((x^\top x_j/d)^2)$ with no $\gamma$. Bernstein bounds each remainder by $d^{-2\theta}$ up to logs, and summing $n$ of them gives a squared error of order $d^{-(4\theta-1)}\log^{4.1}d$. This vanishes iff $\theta>1/4$, that is $m>0$, so the $8+m$ moments are exactly what makes this summed error vanish.

With both surrogates in place, the variance is that of ridge regression on linear features. Insert $A^{-1}A$ with $A=K^{lin}(X,X)$ and use Steps 2 and 3 to get $V\le8\sigma^2\,\mathbb E_\mu\|A^{-1}\beta Xx/d\|^2$ plus the linearization error. Averaging over the test point bounds $X\Sigma X^\top$ by $\|\Sigma\|XX^\top$, and adding the PSD term $\alpha\beta11^\top$ turns the middle matrix into $\beta(A-\gamma I)$, which commutes with $A$, so the trace becomes a function of $A$'s spectrum: $$\mathbb E_\mu\|A^{-1}\beta Xx/d\|^2 \le \frac{\|\Sigma\|}{d}\operatorname{Tr}\big(A^{-1}\beta(A-\gamma I)A^{-1}\big) = \frac{\|\Sigma\|}{d}\sum_j\frac{\lambda_j}{(\gamma/\beta + \lambda_j)^2}.$$ This is the variance of ridge regression with ridge $r=\gamma/\beta$.

The bias is a projection residual, bounded by the tail of the kernel spectrum. In the feature coordinates of the integral operator, Cauchy-Schwarz gives $B\le\|f_*\|_{\mathcal H}^2\int\|P^\perp\phi(x)\|^2d\mu(x)$, with $P$ the projection onto $\operatorname{span}\{\phi(x_i)\}$. That span contains the top-$k$ eigenspace $\hat U_k$ of the empirical feature covariance for every $k\le n$, so the residual is at most the kernel-PCA residual of $\hat U_k$, which on the sample is exactly $\frac1n\sum_{j>k}\lambda_j(K(X,X))$. Nothing here uses the dimension or the linearization.

The population residual tracks the sample residual uniformly over rank-$k$ subspaces. Because $\hat U_k$ is fitted to the same points, bound the gap uniformly over all $U$ with $U^\top U=I_k$, as in Shawe-Taylor and Cristianini's kernel-PCA analysis. By Pythagoras $\|P_U^\perp\phi\|^2=K(x,x)-\langle UU^\top,\phi\phi^\top\rangle_F$, so after McDiarmid and symmetrization the class is linear in $UU^\top$, and Cauchy-Schwarz with $\|UU^\top\|_F=\sqrt k$ bounds its Rademacher complexity by $M\sqrt{k/n}$. A union bound over $k=0,\dots,n$ (at level $\delta/(2n)$) lets $k$ be chosen after seeing the data, and adding this bias bound to Step 4's variance bound gives the theorem.

For the exponential kernel the bound becomes explicit. For the bias, Weyl linearizes the kernel spectrum too, $\lambda_j(K)=\lambda_j(K^{lin})\pm\eta$. The block $\alpha11^\top$ contributes one eigenvalue of order $\alpha n$, which one projected direction removes, adding $2M/\sqrt n$. What remains averages to about $\gamma+\beta\operatorname{Tr}(XX^\top)/(nd)\approx\gamma+\beta\tau$, so $B\lesssim r+\tau$, and the implicit ridge adds bias just as an explicit one would. For the variance with $d>n$ there are at most $n$ nonzero summands, each at most $1/(4r)$; with $n>d$, $t/(c+t)^2\le1/t$ gives $V\le\operatorname{Tr}((X^\top X)^{-1})\le\operatorname{Tr}(\Sigma^{-1})/\lambda_{\min}(Z^\top Z)$ for $Z=X\Sigma^{-1/2}$, and Bai-Yin bounds the last factor.

Localize, discard the empty ball, and split conditional on the design. Take $x_0=0$. The estimator fails to be a kernel average only if no point lands within $h$, which has probability $\le e^{-cnh^d}$ because $\bar p=P(\|X_1\|\le h)\ge c_0p_{\min}h^d$ (this needs a fixed fraction of the $h$-ball inside the support). Otherwise $f_n(0)=\sum_iw_iY_i$ with normalized weights $w_i$, and conditioning on the design separates the noise term $\sigma_\xi^2\sigma_X^2$, $\sigma_X^2=E[\sum_iw_i^2\mathbf 1_E]$, from the bias $b^2(0)=E[(\sum_iw_i(f(X_i)-f(0)))^2\mathbf 1_E]$. By exchangeability $\sigma_X^2=nE[w_1^2\mathbf 1_E]$, so the variance reduces to the second moment of a single weight.

The variance is $O(1/(nh^d))$, and this is the only step that uses $a<d/2$. The singularity only increases the denominator: $K_h\ge1$ on the ball, so $\sum_iK_h(X_i)$ is at least the count $N$ of points in the ball, and $N\le n\bar p/2$ is exponentially rare (there use $w_1\le1$). Off that event the denominator exceeds $n\bar p/2$, and the numerator needs $$E[K_h^2(X_1)]\le p_{\max}h^d\int K^2=\frac{p_{\max}|S^{d-1}|}{d-2a}\,h^d,$$ finite iff $a<d/2$. So $\sigma_X^2\lesssim1/(nh^d)$, with a constant that diverges as $a\uparrow d/2$.

The bias is $O(h^\beta)$, deterministically for $\beta\le1$ and by symmetry for $\beta\in(1,2]$. On $E$ the weights are nonnegative, sum to one and vanish outside the $h$-ball, so for $\beta\le1$ the error is a convex combination of local increments, each at most $L_fh^\beta$; a weight, however large, sits on a point within distance $h$, so the spikes are harmless for the bias, and nothing about $p$ is used. For $\beta>1$ that bound leaves a linear term of size $\|\nabla f(0)\|h$, so it must cancel on average. In each off-diagonal pair term the weight product depends on $x_i$ only through $\|x_i\|$ because $K$ is radial, and the linear part is odd, so its integral against $p(0)$ over the full ball vanishes; what is left is controlled by $|p(x)-p(0)|\le L_ph^{\beta-1}$, and this bound on the whole ball, with $p\ge p_{\min}$ for turning $dx$ integrals into $p$-integrals, is why the ball must lie inside the support. Completing the square, each pair term is at most $$h^{2\beta}\Big(L_f+\frac{\|\nabla f(0)\|L_p}{p_{\min}}\Big)^2E[w_iw_j\mathbf 1_E],$$ and the $n(n-1)$ pairs total $O(h^{2\beta})$ because $\sum_{i\ne j}E[w_iw_j\mathbf 1_E]\le1$, and the diagonal terms are bounded by $\sigma_X^2$.

Balancing the bandwidth gives the minimax rate. The risk is $\lesssim(\sigma_\xi^2+1)/(nh^d)+h^{2\beta}$ plus terms of order $ne^{-cnh^d}$, and $h=n^{-1/(2\beta+d)}$ balances the first two at $n^{-2\beta/(2\beta+d)}$, the textbook kernel-smoothing trade-off, unaffected by interpolation. The lower bound is cited (Stone 1982, Tsybakov 2009), not proved; the textbook argument is Le Cam's two-point method, comparing $f_0=0$ with a bump of height $Lh^\beta$ and width $h$, whose KL divergence over $n$ samples is $O(nh^{2\beta+d})=O(1)$ at this bandwidth.

Consequences and limits. Exact fit is not evidence of overfitting. In LR's high-dimensional regime interpolation with a curved kernel behaves like ridge regression, and the data spectrum produces a U-shaped trade-off with no tuning parameter: fast decay (small $\tau$) makes the implicit ridge tiny and the variance large, slow decay makes the bias $r+\tau$ large, and intermediate decay is best, as their synthetic experiments confirm. Isotropic data is the bad case, since $\tau=1$ makes the bias term order one, and the authors assert, without proof, that the true variance is then large too. In BRT's low-dimensional regime $h$ sets the trade-off and the singularity sets the fit, and the two do not interact: the spikes occupy a set of vanishing measure and enter the risk only through $\int K^2$, and against 1-NN the difference is the effective neighborhood, one point versus $nh^d\to\infty$. Limits: LR is only an upper bound, so it can show that interpolation works but not that it fails; its moment assumption is strong, the linearization error decays only like $d^{-m/(8+m)}$, and with $d\asymp n$ the nonlinear part of $K$ reduces to $\gamma I$, so it explains benign interpolation of the linear component and says nothing about learning nonlinear structure; the corollaries also step outside $d\asymp n$. BRT need $\beta\le2$, a bandwidth tuned to $\beta$, and $a<d/2$, and for $a\ge d/2$ it is the proof that breaks down, with no failure of the estimator shown; this is low-dimensional local smoothing, unlike the global minimum-norm interpolation studied for networks. A consequence not in the paper, from the standard plug-in inequality: for $Y\in\{0,1\}$ and $\beta\le1$, the interpolating plug-in classifier has excess risk $O(n^{-\beta/(2\beta+d)})$; BRT themselves prove no classification result.