Parameter-free Statistically Consistent Interpolation: Dimension-independent Convergence Rates for Hilbert kernel regression
Year: 2021

Abstract:
Previously, statistical textbook wisdom has held that interpolation of noisy training 1 data will lead to poor generalization. However, recent work has shown that this 2 is not true and that good generalization can be obtained with function fits that 3 interpolate training data. This could explain why overparameterized deep nets with 4 zero or small training error do not necessarily overfit and could generalize well. 5 Data interpolation schemes have been exhibited that are provably Bayes optimal in 6 the large sample limit and achieve the theoretical lower bounds for excess risk (Sta7 tistically Consistent Interpolation) in any dimension. These interpolation schemes 8 are non-parametric Nadaraya-Watson style estimators with singular kernels, which 9 exhibit statistical consistency in any data dimension for large sample sizes. The 10 recently proposed weighted interpolating nearest neighbors scheme (wiNN) is in 11 this class, as is the previously studied Hilbert kernel interpolation scheme. In 12 the Hilbert scheme, the regression function estimator for a set of labelled data 13 pairs, (xi, yi) ∈ R × R, i = 0, ..., n, has the form f̂(x) = ∑ i yiwi(x), where 14 wi(x) = ‖x − xi‖/ ∑ j ‖x − xj‖. This interpolating function estimator is 15 unique in being entirely free of parameters and does not require bandwidth se16 lection. While statistical consistency was previously proven for this scheme, the 17 precise convergence rates for the finite sample risk were not established. Here, we 18 carry out a comprehensive study of the asymptotic finite sample behavior of the 19 Hilbert kernel regression scheme and prove a number of relevant theorems. We 20 prove under broad conditions that the excess risk of the Hilbert regression estimator 21 is asymptotically equivalent pointwise to σ(x)/ ln(n) where σ(x) is the noise 22 variance. We also show that the excess risk of the plugin classifier is upper bounded 23 by 2|f(x) − 1/2|1−α (1 + ε)ασα(x)(ln(n))−α2 , for any 0 < α < 1, where f is 24 the regression function x 7→ E[y|x]. Our proofs proceed by deriving asymptotic 25 equivalents of the moments of the weight functions wi(x) for large n, for instance 26 for β > 1, E[w i (x)] ∼n→∞ ((β− 1)n ln(n))−1. We further derive an asymptotic 27 equivalent for the Lagrange function and explicitly exhibit the nontrivial extrapola28 tion properties of this estimator. Notably, the convergence rates are independent of 29 data dimension and the excess risk is dominated by the noise variance. The bias 30 term, for which we also give precise asymptotic estimates, is always subleading 31 when the density of data at the considered point is strictly positive. If this local 32 density is zero, we show that the bias term does not vanish in the limit of a large 33 data set and we compute its limit explicitly. Finally, we present heuristic arguments 34 for a universal w−2 power-law behavior of the probability density of the weights 35 in the large n limit. 36 Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.

text: Previously, statistical textbook wisdom has held that interpolation of noisy training1 data will lead to poor generalization. However, recent work has shown that this2 is not true and that good generalization can be obtained with function fits that3 interpolate training data. This could explain why overparameterized deep nets with4 zero or small training error do not necessarily overfit and could generalize well.5 Data interpolation schemes have been exhibited that are provably Bayes optimal in6 the large sample limit and achieve the theoretical lower bounds for excess risk (Sta-7 tistically Consistent Interpolation) in any dimension. These interpolation schemes8 are non-parametric Nadaraya-Watson style estimators with singular kernels, which9 exhibit statistical consistency in any data dimension for large sample sizes. The10 recently proposed weighted interpolating nearest neighbors scheme (wiNN) is in11 this class, as is the previously studied Hilbert kernel interpolation scheme. In12 the Hilbert scheme, the regression function estimator for a set of labelled data13 pairs, (xi, yi) ∈ Rd × R, i = 0, ..., n, has the form f̂(x) = ∑ i yiwi(x), where14
wi(x) = ‖x − xi‖−d/ ∑ j ‖x − xj‖−d. This interpolating function estimator is15 unique in being entirely free of parameters and does not require bandwidth se-16 lection. While statistical consistency was previously proven for this scheme, the17 precise convergence rates for the finite sample risk were not established. Here, we18 carry out a comprehensive study of the asymptotic finite sample behavior of the19 Hilbert kernel regression scheme and prove a number of relevant theorems. We20 prove under broad conditions that the excess risk of the Hilbert regression estimator21 is asymptotically equivalent pointwise to σ2(x)/ ln(n) where σ2(x) is the noise22 variance. We also show that the excess risk of the plugin classifier is upper bounded23 by 2|f(x) − 1/2|1−α (1 + ε)ασα(x)(ln(n))−α2 , for any 0 < α < 1, where f is24 the regression function x 7→ E[y|x]. Our proofs proceed by deriving asymptotic25 equivalents of the moments of the weight functions wi(x) for large n, for instance26 for β > 1, E[wβi (x)] ∼n→∞ ((β− 1)n ln(n))−1. We further derive an asymptotic27 equivalent for the Lagrange function and explicitly exhibit the nontrivial extrapola-28 tion properties of this estimator. Notably, the convergence rates are independent of29 data dimension and the excess risk is dominated by the noise variance. The bias30 term, for which we also give precise asymptotic estimates, is always subleading31 when the density of data at the considered point is strictly positive. If this local32 density is zero, we show that the bias term does not vanish in the limit of a large33 data set and we compute its limit explicitly. Finally, we present heuristic arguments34 for a universal w−2 power-law behavior of the probability density of the weights35 in the large n limit.36
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.
1 Introduction37
Data interpolation and statistical regression of noisy data are both classical subjects but their domain38 of application have been disjoint until recently. Scattered data interpolation techniques [1] are39 generally used for clean data. On the other hand, when supervised learning or statistical regression40 techniques are applied to noisy data, in general smoothing or regularization methods are applied41 to prevent training data interpolation, as the latter is believed to lead to poor generalization [2].42 However, accumulating empirical evidence from overparameterized deep networks has shown that43 data interpolation (equivalently, zero error on the training set) does not automatically imply poor44 generalization [3, 4]. This has in turn given rise to a rapidly growing body of theoretical work to45 understand how and why noisy data interpolation can still lead to good generalization [5, 6, 7, 8, 9,46 10, 11, 12, 13, 14, 15].47
A key observations in this regard is the phenomenon of Statistically Consistent Interpolation [16],48 i.e., regression function estimation that interpolates training data but also generalizes as well as49 possible by achieving the Bayes limit for expected generalization error (risk) when the sample size50 becomes large. This hints at a rich set of theoretical questions at the interface between the disciplines51 of scattered data interpolation and supervised learning, that have only begun to be addressed. In52 particular, there has been comparatively little study of the generalization error or risk of interpolating53 learners. Computation of generalization error bounds in machine learning often relies on the capacity54 of the class of fitting functions [17], however such model complexity based bounds are not tight55 enough to be useful for interpolating learners [4]. For nonparametric interpolation approaches such as56 that considered here, it is also not clear what model complexity means. Thus, there is a need for other57 approaches to understanding the generalization behavior of nonparametric interpolating learners,58 including more direct treatments of the generalization error for specific interpolation schemes so as59 to gain better theoretical understanding. The current paper addresses this need.60
We present a detailed analysis of the finite-sample risk of an interpolating learner with intriguing61 theoretical properties, the Hilbert kernel estimator (Devroye et. al. [18]). A unique property of62 this Nadaraya-Watson (NW) style estimator [19, 20] is that it is fully parameter-free and does not63 have any bandwidth or scale parameter. It is global and uses all data points for each estimate: the64 associated kernel is a power law and thus scale-free. Although statistical consistency of this estimator65 was proven [18] when it was proposed, there has been no systematic analysis of the associated66 convergence rates and asymptotic finite sample behavior. We provide this analysis in the present67 study.68
Related work The only other interpolation scheme we are aware of, that is proven to be statistically69 consistent in arbitrary dimensions under general conditions, is the recently proposed weighted70 interpolating nearest neighbors method (wiNN) [7], which is also a NW estimator utilizing a singular71 power law kernel of a very similar form but with two important differences: a finite number of72 neighbors k is utilized (rather than all data points), and the power law exponent δ of the NW kernel73 satisfies 0 < δ < d/2 rather than δ = d. To achieve consistency k has to scale appropriately with74 sample size. Despite the superficial resemblance, the wiNN and Hilbert Kernel estimators have quite75 different convergence rates, as we will see from the results of this paper. Also worth mentioning is the76 Shepard interpolation scheme [21] originally proposed for interpolation of 2D geospatial data sets,77 also a NW style interpolating estimator, though used in the context of scattered data interpolation.78 In scattered data interpolation [1], the focus is generally on the approximation error (corresponding79 to the “bias” term in our analysis below). The approximation error of the Shepard scheme has been80 analyzed [22] but as we will see below the risk for Hilbert kernel interpolation is dominated by the81 noise or “variance” term. In contrast with wiNN or Hilbert kernel interpolation, other interpolating82 learning methods such as simplex interpolation [7] or ridgeless kernel regression [11] are generally83 not statistically consistent in fixed finite dimension [8].84
Summary of results of this paper Notation and assumptions pertaining to this summary are defined85 in the problem setup section below. We prove under broad conditions that the excess risk of the86 Hilbert regression estimator is asymptotically equivalent pointwise to σ2(x)/ ln(n) where σ2(x) is87 the noise variance. We also show that the excess risk of the plugin classifier is upper bounded by88 2|f(x)− 1/2|1−α (1 + ε)ασα(x)(ln(n))−α2 , for any 0 < α < 1, where f is the regression function89 x 7→ E[y|x]. Our proofs proceed by deriving asymptotic equivalents of the moments of the weight90 functions wi(x) for large n, for instance for β > 1, E[wβi (x)] ∼n→∞ ((β − 1)n ln(n))−1. We91 further derive an asymptotic equivalent for the Lagrange function and explicitly exhibit the nontrivial92
extrapolation properties of this estimator. Notably, the convergence rates are independent of data93 dimension and the excess risk is dominated by the noise variance. The bias term, for which we also94 give precise asymptotic estimates, is always subleading when the density of data at the considered95 point is strictly positive. If this local density is zero, we show that the bias term does not vanish in the96 limit of a large data set and we compute its limit explicitly. Finally, we present heuristic arguments97 for a universal w−2 power-law behavior of the probability density of the weights in the large n limit.98
2 Problem setup99
Notation, Definitions, Statistical Model We model the labelled training data set100 (x0, y0), . . . , (xn, yn) as n + 1 i.i.d. observations of a random vector (X,Y ) with values in101 Rd × R for regression, and with values in Rd × {0, 1} for binary classification. Due to the indepen-102 dence property, the collection X0, . . . , Xn has the product density ∏n i=0 ρ(xi). We will denote by E103 an expectation over the collection of n+ 1 random vectors and by EX the expectation over the col-104 lection X0, . . . , Xn. An expectation over the same collection while holding Xi = xi will be denoted105 EX|xi . The regression function f : Rd → R is defined as the conditional mean of Y given X = x,106 f(x) := E[Y | X = x] and the conditional variance function is σ2(x) := E[|Y − f(X)|2|X = x].107 f minimizes the expected value of the mean squared prediction error (risk under squared loss),108 f = arg minRsq(h) where Rsq(h) := E[(h(X)− Y )2]. Given any regression estimator f̂(x) the109 corresponding risk can be decomposed as E[Rsq(f̂(X))] = Rsq(f) + E[(f̂(X) − f(X)2]. The110 excess risk is given by Rsq(f̂) − Rsq(f) = E[(f̂(X) − f(X))2]. For a consistent estimator this111 excess risk goes to zero as n→∞ and we are interested in characterizing the rate at which it goes to112 zero with increasing n (note our sample size is n + 1 for notational simplicity but for large n this113 does not change the rate).114
In the case of binary classification, Y ∈ {0, 1} and f(x) = P[Y = 1 | X = x)]. Let F : Rd → {0, 1}115 denote the Bayes optimal classifier, defined by F (x) := θ(f(x)− 1/2) where θ(·) is the Heaviside116 theta function. This classifier minimizes the riskR0/1(h) := E[1{h(X)6=Y }] = P(h(X) 6= Y ) under117 zero-one loss. Given the regression estimator f̂ , we consider the plugin classifier F̂ (x) = θ(f̂(x)− 12 ).118 The classification risk for the plugin classifier F̂ is bounded as E[R0/1(F̂ (x))] − R0/1(F (x)) ≤119
2E[|f̂(x)− f(x)|] ≤ 2 √ E[(f̂(x)− f(x))2].120
Finally, we define two sequences an, bn > 0, n ∈ N, to be asymptotically equivalent for n→ +∞,121 denoted an ∼n→+∞ bn, if the limit of their ratio exists and limn→∞ an/bn = 1.122
In summary, our work will focus on the estimation of asymptotic equivalents for E[(f̂(x)− f(x))2]123 and other relevant quantities as this determines the rate at which the excess risk goes to zero for124 regression, and bounds the rate at which the excess risk goes to zero for classification.125
Assumptions. We define the support Ω of the density ρ as Ω = {x ∈ Rd/ρ(x) > 0}, the closed126 support Ω̄ as the closure of Ω, and Ω◦ as the interior of Ω. Our results will not assume any compactness127 condition on Ω or Ω̄. The boundary of Ω is then defined as ∂Ω = Ω̄ \ Ω◦. We assume that ρ has a128 finite variance σ2ρ. In addition, we will most of the time assume that the density ρ is continuous at the129 considered point x ∈ Ω◦, and in some cases, x ∈ ∂Ω ∩ Ω.130 For the regression function f , we will obtain results assuming either of the following conditions131
• CfCont: f is continuous at the considered x,132
• CfHolder: for all x ∈ Ω◦, there exist αx > 0, Kx > 0, and δx > 0, such that133 x′ ∈ Ω and ‖x− x′‖ ≤ δx =⇒ |f(x)− f(x′)| ≤ Kx ‖x− x′‖αx134 (local Hölder smoothness condition),135
where condition CfHolder is obviously stronger than C f Cont. In addition, we will always assume a136 growth condition for the regression function f :137
• CfGrowth: ∫ ρ(y) f 2(y) 1+‖y‖2d d dy <∞.138
As for the variance function σ, we will obtain results assuming either that σ is bounded or satisfies a139 growth condition similar to the one above140
• CσBound: there exists σ 2 0 ≥ 0, such that, for all x ∈ Ω, we have σ2(x) ≤ σ20 ,141 • CσGrowth: ∫ ρ(y) σ 2(y) 1+‖y‖2d d dy <∞.142
When we will assume condition CσGrowth (obviously satisfied when σ 2 is bounded), we will also143 assume a continuity condition CσCont for σ at the considered x.144
Note that all our results can be readily extended in the case where x ∈ ∂Ω = Ω̄ \ Ω◦ but keeping145 the condition ρ(x) > 0 (i.e., x ∈ ∂Ω ∩ Ω), and assuming the continuity at x of ρ as seen as a146 function restricted to Ω, i.e., limy∈Ω→x ρ(y) = ρ(x). Useful examples are when the support Ω of147 ρ is a d-dimensional sphere or hypercube and x is on the surface of Ω (but still with ρ(x) > 0). To148 guarantee these results for x ∈ ∂Ω ∩ Ω, we need also to assume the continuity at x of f , and assume149 that Ω is smooth enough near x, so that there exists a strictly positive local solid angle ωx defined by150
ωx = lim r→0
1
Vdρ(x)rd ∫ ‖x−y‖≤r ρ(y) ddy = lim r→0 1 Vdrd ∫ y∈Ω/‖x−y‖≤r ddy, (1)
where Vd = Sd/d = πd/2/Γ(d/2 + 1) is the volume of the unit ball in d dimensions, and the second151 inequality results from the continuity of ρ at x. If x ∈ Ω◦, we have ωx = 1, while for x ∈ ∂Ω, we152 have 0 ≤ ωx ≤ 1. For instance, if x is on the surface of a sphere or on the interior of a face of a153 hypercube (and in general, when the boundary near x is locally an hyperplane), we have ωx = 12 . If x154 is a corner of a hypercube, we have ωx = 12d . From our methods of proof presented in the appendix,155 it should be clear that all our results for x ∈ Ω◦ perfectly generalize to any x ∈ ∂Ω ∩ Ω for which156 ωx > 0, by simply replacing Vd whenever it appears in our different results by ωxVd.157
Hilbert kernel interpolating estimator and Bias-Variance decomposition. The Hilbert kernel158 regression estimator f̂(x) is a Nadaraya-Watson style estimator employing a singular kernel:159
wi(x) = ‖x− xi‖−d∑n j=0 ‖x− xj‖−d , (2)
f̂(x) = n∑ i=0 wi(x)yi. (3)
The weights wi(x) are also called Lagrange functions in the interpolation literature and satisfy the160 interpolation property wi(xj) = δij , where δij = 1, if i = j, and 0 otherwise. At any given point161 x, they provide a partition of unity so that ∑n i=0 wi(x) = 1. The mean squared error between the162 Hilbert estimator and the true regression function has a bias-variance decomposition (using the i.i.d163 condition and the earlier definitions)164
f̂(x)− f(x) = n∑ i=0 wi(x)[f(xi)− f(x)] + n∑ i=0 wi(x)[yi − f(xi)], (4)
E[(f̂(x)− f(x))2] = B(x) + V(x), (5) (Bias) B(x) = EX [( n∑
i=0
wi(x)[f(xi)− f(x)] )2] , (6)
(V ariance) V(x) = E [ n∑ i=0 w2i (x)[yi − f(xi)]2 ] = EX [ n∑ i=0 w2i (x)σ 2(xi) ] . (7)
The present work derives asymptotic behaviors and bounds for the regression and classification risk165 of the Hilbert estimator for large sample size n. These results are derived by analyzing the large n166 behaviors of the bias and variance terms, which in turn depend on the behavior of the moments of the167 weights or the Lagrange functions wi(x). For all these quantities, asymptotically equivalent forms168 are derived. The proofs exploit a simple integral form of the weight function and details are provided169 in the appendix, while the body of the paper provides the results and associated discussions.170
3 Results171
3.1 The weights, variance and bias terms172
3.1.1 Moments of the weights: large n behavior173
In this section, we consider the moments and the distribution of the weights wi(x) at a given point x.174 The first moment is simple to compute. Since the weights sum to 1 and Xi are i.i.d, it follows that175 EX|xi [wi(x)] are all equal and thus EX|xi [wi(x)] = (n+ 1)−1. The other moments are much less176 trivial to compute and we prove the following theorem in the appendix A.2:177 Theorem 3.1. For x ∈ Ω◦ (so that ρ(x) > 0), we assume ρ continuous at x. Then, the moments of178 the weight w0(x) satisfy the following properties:179
• For β > 1:180
E [ wβ0 (x) ] ∼
n→+∞
1
(β − 1)n ln(n) . (8)
• For 0 < β < 1: defining κβ(x) := ∫ ρ(x+y) ||y||βd d dy <∞, we have181
E [ wβ0 (x) ] ∼
n→+∞
κβ(x)
(Vdρ(x)n ln(n))β . (9)
• For β < 0: all moments for β ≤ −1 are infinite, and the moments of order −1 < β < 0182 satisfy183 E [ wβ0 (x) ] ≤ 1 + nκ|β|(x)κβ(x), (10)
so that a sufficient condition for its existence is κβ(x) = ∫ ρ(x+ y)||y|||β|d ddy <∞.184
Heuristically, the behavior of these moments are consistent with the random variable W = w0(x)185 having a probability distribution satisfying a scaling relation P (W ) = 1Wn p ( W Wn ) , with the scaling186
function p having the universal tail (i.e., independent of x and ρ), p(w) ∼ w→+∞ w−2, and a scale Wn187 expected to vanish with n, when n→ +∞. With this assumption, we can determine the scale Wn by188 imposing the exact condition E[W ] = 1/(n+ 1) ∼ 1/n:189
E[W ] = 1
Wn ∫ 1 0 p ( W Wn ) W dW = Wn ∫ 1 Wn 0 p(w)w dw (11)
∼ Wn ∫ 1 Wn
1
dw
w ∼ −Wn ln(Wn) ∼
1 n , (12)
leading to Wn ∼ 1n ln(n) . Then, the moment of order β > 1 is given by190
E[W β ] = 1
Wn ∫ 1 0 p ( W Wn ) W β dW ∼Wn ∫ 1 0 W β−2 dW ∼ n→+∞
1
(β − 1)n ln(n) , (13)
which indeed coincides with the first result of Theorem 3.1. Our heuristic argument also suggests191 that in the case 0 < β < 1, we have192
E[W ] = 1
Wn ∫ 1 0 p ( W Wn ) W β dW ∼ n→+∞ ∫ +∞ 0 p (w)wβ dw (n ln(n))β , (14)
where the last integral converges since p(w) ∼ w→+∞ w−2 and β < 1. This result is perfectly consistent193 with Eq. (9) in Theorem 3.1, and suggests that ∫ +∞
0 p (w)wβ dw = κβ(x) (Vdρ(x))β . Interestingly, for194 0 < β < 1, and contrary to the case β > 1, we find that the large n equivalent of the moment is not195 universal and depends explicitly on x and the density ρ. As for moments of order −1 < β < 0, we196 conjecture that they are still given by Eq. (9) (and equivalently, by Eq. (14)) provided they exist, and197 that the sufficient condition for their existence κβ(x) <∞ is hence also necessary, since κβ(x) also198 appears in Eq. (9). The fact that moments for β ≤ −1 do not exist strongly suggests that p(0) > 0.199
In fact, Eq. (14)) also suggests that all moments for −1 < β < 0 exist if and only if 0 < p(0) <∞.200 In the Fig. 2 of the appendix, we present numerical simulations confirming our scaling ansatz, the201 fact that p(w) ∼
w→+∞ w−2, and the quantitative prediction for Wn.202
It is shown in Devroye et al. [18] that the Hilbert kernel regression estimate does not converge203 almost surely (a.s.) by giving a specific example. Insight can be gained into this lack of almost sure204 convergence by considering the weight function w0(x), for a sequence of independent training sample205 sets of increasing size n+ 1. Let the corresponding sequence of weights be denoted as ωn ∈ [0, 1].206 From Theorem 3.1, it is clear that ωn converges to zero in probability, since the following Chebyshev207 bound holds (analogous to the bound on the regression risk):208
P(ωn > ε) ≤ 1 + δ
ε2n ln(n) , (15)
for arbitrary ε > 0 and δ > 0, and for n larger than some constant Nx,δ. Alternatively, one can209 exploit the fact that E[ωn] = 1n+1 , leading to P(ωn > ε) ≤ 1 ε n , which is less stringent than Eq. (15)210 as far as the n-dependence is concerned, but is more stringent for the ε-dependence of the bounds.211
Let us show heuristically that ωn does not converge a.s. to zero. Consider the infinite sequence of212 events En ≡ {ωn > ε}, n ∈ N, and the corresponding infinite sum ∑ n P(En) = ∑ n P(ωn > ε).213 Exploiting our previous heuristic argument for the scaling form of the distribution of weights, we214 obtain215
P(ωn > ε) = ∫ 1 ε 1 Wn p ( W Wn ) dW ∼ ∫ n lnn εn lnn dw w2 ∼ 1− ε ε n ln(n) . (16)
Since ∑N n=2 1 n ln(n) ∼ ln(ln(N)) is a divergent series, a Borel-Cantelli argument suggests that an216 infinite number of the events En (i.e., ωn > ε) must occur, which implies that ωn does not converge217 a.s. to 0. Note that the weights are equal to 1 at the data points due to the interpolation condition, so218 that large weights occasionally occur, causing the lack of a.s. convergence.219
3.1.2 Lagrange function: scaling limit220
The expected value of the Lagrange functions wi(x) have a simple form in the large n limit. Due221 to the i.i.d. condition the indices i are exchangeable and we set i = 0 for the computation of the222 expected Lagrange function L0(x) = EX|x0 [w0(x)]. Thus, one of the sample points (denoted x0) is223 held fixed and the other ones are averaged over in computing the expected Lagrange function. For224 x0 6= x kept fixed, we have limn→∞ L0(x) = 0. However, we show in the appendix A.3 that L0(x)225 takes a very simple form when taking a specific scaling limit:226 Theorem 3.2. For x ∈ Ω◦, we assume ρ continuous at x. Then, in the limit (denoted by limZ ), n→227 +∞, ‖x−x0‖−d → +∞ (i.e., x0 → x), and such that zx(n, x0) = Vdρ(x)‖x−x0‖dn log(n)→ Z,228 the Lagrange function L0(x) = EX|x0 [w0(x)] converges to a proper limit,229
lim Z L0(x) =
1
1 + Z . (17)
The proof of this theorem shows that the relative error between L0(x) and 11+Z for finite but large n230 and large ‖x− x0‖−d, such that zx(n, x0) remains close to Z, is O(1/ ln(n)).231 Exploiting Theorem 3.2, we can use a simple heuristic argument to estimate the tail of the distribution232 of the random variable W = w0(x). Indeed, approximating L0(x) for finite but large n by its233 asymptotic form 11+zx(n,x0) , with zx(n, x0) = Vdρ(x)n log(n)‖x− x0‖
d, we obtain234 ∫ 1 W P (W ′) dW ′ ∼ ∫ ρ(x0) θ ( 1 1 + Vdρ(x)n log(n)‖x− x0‖d −W ) ddx0, (18)
∼ Vdρ(x) ∫ +∞
0
θ
( 1
1 + Vdρ(x)n log(n)u −W
) du, (19)
∼ 1 n ln(n)W
=⇒ P (W ) ∼ 1 n ln(n)W 2 , (20)
where θ(.) is the Heaviside function. This heuristic result is again perfectly consistent with our guess235 of the previous section that P (W ) = 1Wn p ( W Wn ) , with the scaling function p having the universal236
tail, p(w) ∼ w→+∞ w−2, and a scale Wn ∼ 1n ln(n) . Indeed, in this case and in the limit n→ +∞, we237 obtain that P (W ) ∼ 1Wn ( Wn W )2 ∼ WnW 2 ∼ 1 n ln(n)W 2 , which is identical to the result of Eq. (20).238
3.1.3 The variance term239
A simple application of the result of Theorem 3.1 for β = 2 (see appendix A.4) allows us to bound240 the variance term V(x) = E [∑n i=0 w 2 i (x)[yi − f(xi)]2 ] for a bounded variance function σ2:241
Theorem 3.3. For x ∈ Ω◦, ρ continuous at x, σ2 ≤ σ20 , and for any ε > 0, there exists a constant242 Nx,ε such that for n ≥ Nx,ε, we have243
V(x) ≤ (1 + ε) σ 2 0
ln(n) . (21)
Relaxing the boundedness condition for σ, but assuming the continuity of σ2 at x along with a growth244 condition, allows us to obtain a precise asymptotic equivalent of V(x), when n→ +∞:245 Theorem 3.4. For x ∈ Ω◦, σ(x) > 0, ρσ2 continuous at x, and assuming the condition CσGrowth,246 i.e., ∫ ρ(y) σ
2(y) 1+‖y‖2d d dy <∞, we have247
V(x) ∼ n→+∞
σ2(x) ln(n) . (22)
Note that if the mean variance ∫ ρ(y)σ2(y) ddy < ∞, which is in particular the case when σ2 is248 bounded over Ω, then the condition CσGrowth is in fact automatically satisfied.249
3.1.4 The bias term250
In appendix A.5, we prove the following three theorems for the bias term.251
Theorem 3.5. For x ∈ Ω◦ (so that ρ(x) > 0), we assume that ρ is continuous at x, and the conditions252
• CfGrowth: ∫ ρ(y) f 2(y) 1+‖y‖2d d dy <∞,253
• CfHolder: there exist αx > 0, Kx > 0, and δx > 0, such that254 x′ ∈ Ω and ‖x− x′‖ ≤ δx =⇒ |f(x)− f(x′)| ≤ Kx ‖x− x′‖αx255 (local Hölder condition for f ).256
Moreover, we define κ(x) = ∫ ρ(x+ y) f(x+y)−f(x)||y||d d dy, where we have |κ(x)| <∞.257
Then, for κ(x) 6= 0, the bias term B(x) = EX [(∑n i=0 wi(x)[f(xi)− f(x)] )2] satisfies258
B(x) ∼ n→+∞
( E [ f̂(x) ] − f(x) )2 , with E [ f̂(x) ] − f(x) ∼
n→+∞
κ(x)
Vdρ(x) ln(n) . (23)
In the non generic case κ(x) = 0, we have the weaker result259
B(x) =  O ( n− 2αx d (ln(n))−1− 2αx d ) , for d > 2αx O ( n−1(ln(n))−1 ) , for d = 2αx O ( n−1(ln(n))−2 ) , for d < 2αx
(24)
Note that κ(x) = 0 is non generic but can still happen, even if f is not constant. For instance, if Ω is a260 sphere centered at x or Ω = Rd, if ρ(x+y) = ρ̂(||y||) is isotropic around x, and if fx : y 7→ f(x+y)261 is an odd function of y, then we indeed have κ(x) = 0 at this symmetric point x.262
Interestingly, for κ(x) 6= 0, Eq. (23) shows that the bias B(x) is asymptotically dominated by263 the square of E [ f̂(x) ] − f(x), showing that the fluctuations of E [ f̂(x) ] − ∑n i=0 wi(x)f(xi) are264
negligible compared to E [ f̂(x) ] − f(x), in the limit n→ +∞ and for κ(x) 6= 0.265
One can relax the local Hölder condition, but at the price of a weaker estimate for B(x) which will266 however be enough to obtain strong results for the regression and classification risks (see below):267
Theorem 3.6. For x ∈ Ω◦, we assume ρ and f continuous at x, and the growth condition CfGrowth:268 ∫ ρ(y) f
2(y) 1+‖y‖2d d dy <∞. Then, the bias term satisfies269
B(x) = o ( 1
ln(n)
) , (25)
or equivalently, for any ε > 0, there exists Nx,ε, such that for n ≥ Nx,ε270
B(x) ≤ ε ln(n) . (26)
Let us now consider a point x ∈ ∂Ω for which we have ρ(x) = 0 (note that x ∈ ∂Ω does not271 necessarily imply ρ(x) = 0). In appendix A.5, we show the following theorem for the expectation272 value of the estimator f̂(x) in the limit n→ +∞:273 Theorem 3.7. For x ∈ ∂Ω such that ρ(x) = 0, we assume that f and ρ satisfy the conditions274
• CfGrowth: ∫ ρ(y) |f(y)| 1+‖y‖d d dy <∞,275
• CρHolder: there exist αx > 0, Kx > 0, and δx > 0, such that276 x′ ∈ Ω and ‖x− x′‖ ≤ δx =⇒ |ρ(x′)| ≤ Kx ‖x− x′‖αx277 (local Hölder condition for ρ).278
Moreover, we define κ(x) = ∫ ρ(x + y) f(x+y)−f(x)||y||d d
dy (|κ(x)| < ∞ under condition CfGrowth),279 and λ(x) = ∫ ρ(x+y) ||y||d d dy (0 < λ(x) <∞ under condition CσHolder). Then,280
lim n→+∞ E[f̂(x)]− f(x) = κ(x) λ(x) . (27)
Hence, in the generic case κ(x) 6= 0 (see Theorem 3.5 and the discussion below it) and under281 condition CρHolder, we find that the bias does not vanish when ρ(x) = 0, and that the estimator f̂(x)282 does not converge to f(x). When ρ(x) = 0, the scarcity of data near the point x indeed prevents the283 estimator to converge to the actual value of f(x). In appendix A.5, we show an example of a density284 ρ continuous at x and such that ρ(x) = 0, but not satisfying the condition CρHolder, and for which285 limn→+∞ E[f̂(x)] = f(x), even if κ(x) 6= 0.286
3.2 Asymptotic equivalent for the regression risk287
In appendix A.6, we prove the following theorem establishing the asymptotic rate at which the excess288 risk goes to zero with large sample size n for Hilbert kernel regression, under mild conditions that do289 not require f or σ to be bounded, but only to satisfy some growth conditions:290
Theorem 3.8. For x ∈ Ω◦, we assume σ(x) > 0, ρ, σ, and f continuous at x, and the growth291 conditions CσGrowth: ∫ ρ(y) σ 2(y) 1+‖y‖2d d dy <∞ and CfGrowth: ∫ ρ(y) f 2(y) 1+‖y‖2d d dy <∞.292
Then the following statements are true:293
• The excess regression risk at the point x satisfies294
E[(f̂(x)− f(x))2] ∼ n→+∞
σ2(x) ln(n) . (28)
• The Hilbert kernel estimate converges pointwise to the regression function in probability.295 More specifically, for any δ > 0, there exists a constant Nx,δ, such that for any ε > 0, we296 have the following Chebyshev bound, valid for n ≥ Nx,δ297
P[|f̂(x)− f(x)| ≥ ε] ≤ 1 + δ ε2
σ2(x) ln(n) . (29)
This theorem is a consequence of the corresponding asymptotically equivalent forms of the variance298 and bias terms presented above. Note that as long as ρ(x) > 0, the variance term dominates over the299 bias term and the regression risk has the same form as the variance term.300
3.3 Rates for the plugin classifier301
In appendix A.7, we prove the following theorem establishing the asymptotic rate at which the302 classification risk goes to zero with large sample size n for Hilbert kernel regression:303 Theorem 3.9. For x ∈ Ω◦, we assume σ(x) > 0, ρ, σ, and f continuous at x. Then, the classification304 risk E[R0/1(F̂ (x))]−R0/1(F (x)) vanishes for n→ +∞.305 More precisely, for any ε > 0, there exists Nx,ε, such that for any n ≥ Nx,ε,306
0 ≤ E[R0/1(F̂ (x))]−R0/1(F (x)) ≤ 2(1 + ε) σ(x)√ ln(n) , (30)
In addition, for any 0 < α < 1, the general inequality307 E[R0/1(F̂ (x))]−R0/1(F (x)) ≤ 2|f(x)− 1/2|1−α E [ |f̂(x)− f(x)|2 ]α 2 , (31)
holds unconditionally and, for n ≥ Nx,ε, leads to308
0 ≤ E[R0/1(F̂ (x))]−R0/1(F (x)) ≤ 2|f(x)− 1/2|1−α (1 + ε)α σα(x)
(ln(n)) α 2 . (32)
For 0 < α < 1, Eq. (32) is weaker than Eq. (30) in terms of its dependence on n, but explicitly309 shows that the classification risk vanishes for f(x) = 1/2. This theorem does not require any growth310 condition for f or σ, since both functions takes values in [0, 1] in the classification context.311

heading: 3.4 Extrapolation behavior outside the support of ρ312
text: We now take the point x outside the closed support Ω̄ of the distribution ρ (which excludes the case313 Ω = Rd). We are interested in the behavior of E [ f̂(x) ] as n→ +∞. In appendix A.8 we prove:314
Theorem 3.10. For x /∈ Ω̄, we assume the growth condition ∫ ρ(y) |f(y)|
1+‖y‖d d dy <∞. Then,315
f̂∞(x) := lim n→+∞
E [ f̂(x) ] = ∫ ρ(y)f(y)‖x− y‖−d ddy∫ ρ(y)‖x− y‖−d ddy , (33)
and f̂∞ is continuous at all x /∈ Ω̄.316 In addition, if ∫ ρ(y)|f(y)| ddy <∞, and defining d(x,Ω) > 0 as the distance between x and Ω, we317 have318 lim
d(x,Ω)→+∞ f̂∞(x) =
∫ ρ(y)f(y) ddy. (34)
Finally, we consider x0 ∈ ∂Ω such that ρ(x0) > 0 (i.e., x0 ∈ ∂Ω ∩ Ω), and assume319 that f and ρ seen as functions restricted to Ω are continuous at x0, i.e. limy∈Ω→x0 ρ(y) =320 ρ(x0) and limy∈Ω→x0 f(y) = f(x0). We also assume that the local solid angle ω0 =321 limr→0
1 Vdρ(x0)rd ∫ ‖x0−y‖≤r ρ(y) d dy exists and satisfies ω0 > 0. Then,322
lim x/∈Ω̄→x0 f̂∞(x) = f(x0). (35)
Eq. (34) shows that far away from Ω (which is possible to realize, for instance, when Ω is bounded),323 f̂∞(x) goes smoothly to the ρ-mean of f . Moreover, Eq. (35) establishes a continuity property for324 the extrapolation f̂∞ at x0 ∈ ∂Ω ∩ Ω under the stated conditions (remember that for x ∈ Ω◦, we325 have limn→+∞ E [ f̂(x) ] = f(x); see Theorem 3.5, and in particular Eq. (23)).326
References327 [1] Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.328
[2] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani. An introduction to329 statistical learning, volume 112. Springer, 2013.330
[3] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals. Understanding331 deep learning requires rethinking generalization. In International Conference on Learning332 Representations, 2017.333
[4] Mikhail Belkin, Siyuan Ma, and Soumik Mandal. To understand deep learning we need to334 understand kernel learning. In Proceedings of the 35th International Conference on Machine335 Learning, pages 541–549, 2018.336
[5] Adele Cutler and Guohua Zhao. Pert-perfect random tree ensembles. Computing Science and337 Statistics, 33:490–497, 2001.338
[6] Abraham J Wyner, Matthew Olson, Justin Bleich, and David Mease. Explaining the success of339 adaboost and random forests as interpolating classifiers. Journal of Machine Learning Research,340 18(48):1–33, 2017.341
[7] Mikhail Belkin, Daniel Hsu, and Partha Mitra. Overfitting or perfect fitting? risk bounds for342 classification and regression rules that interpolate. arXiv preprint arXiv:1806.05161, 2018.343
[8] Alexander Rakhlin and Xiyu Zhai. Consistency of interpolation with laplace kernels is a344 high-dimensional phenomenon. arXiv preprint arXiv:1812.11167, 2018.345
[9] Greg Ongie, Rebecca Willett, Daniel Soudry, and Nathan Srebro. A function space view of346 bounded norm infinite width relu nets: The multivariate case. arXiv preprint arXiv:1910.01635,347 2019.348
[10] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machine-349 learning practice and the classical bias–variance trade-off. Proceedings of the National Academy350 of Sciences, 116(32):15849–15854, 2019.351
[11] Tengyuan Liang and Alexander Rakhlin. Just interpolate: Kernel" ridgeless" regression can352 generalize. arXiv preprint arXiv:1808.00387, 2018.353
[12] Peter L Bartlett, Philip M Long, Gábor Lugosi, and Alexander Tsigler. Benign overfitting in354 linear regression. arXiv preprint arXiv:1906.11300, 2019.355
[13] Andrea Montanari, Feng Ruan, Youngtak Sohn, and Jun Yan. The generalization error of356 max-margin linear classifiers: High-dimensional asymptotics in the overparametrized regime.357 arXiv preprint arXiv:1911.01544, 2019.358
[14] Mina Karzand and Robert D Nowak. Active learning in the overparameterized and interpolating359 regime. arXiv preprint arXiv:1905.12782, 2019.360
[15] Yue Xing, Qifan Song, and Guang Cheng. Statistical optimality of interpolated nearest neighbor361 algorithms. arXiv preprint arXiv:1810.02814, 2018.362
[16] Partha P. Mitra. Fitting elephants in modern machine learning by statistically consistent363 interpolation. Nature Machine Intelligence, 3(5):378–386, May 2021.364
[17] Martin Anthony and Peter L Bartlett. Neural Network Learning: Theoretical Foundations.365 Cambridge University Press, 1999.366
[18] Luc Devroye, Laszlo Györfi, and Adam Krzyżak. The hilbert kernel regression estimate. Journal367 of Multivariate Analysis, 65(2):209–227, 1998.368
[19] Elizbar A Nadaraya. On estimating regression. Theory of Probability & Its Applications,369 9(1):141–142, 1964.370
[20] GS Watson. Smooth regression analysis. Sankhya A: 26: 359-372, (50), 1964.371
[21] Donald Shepard. A two-dimensional interpolation function for irregularly-spaced data. In372 Proceedings of the 1968 23rd ACM national conference, pages 517–524. ACM, 1968.373
[22] Reinhard Farwig. Rate of convergence of shepard’s global interpolation formula. Mathematics374 of Computation, 46(174):577–590, 1986.375
Checklist376
1. For all authors...377 (a) Do the main claims made in the abstract and introduction accurately reflect the paper’s378 contributions and scope? [Yes] Brief statement of the primary theorems are included in the379 abstract together with verbal descriptions of the principal results.380
(b) Did you describe the limitations of your work? [Yes] Theorems include limiting conditions.381 (c) Did you discuss any potential negative societal impacts of your work? [No] The paper382 consists of theorems and proofs about the convergence rates of an interpolation scheme and it383 is hard to conceive any negative social impact of this mathematical exercise. If anything it384 might have a positive impact by shedding theoretical light on the interpolation of noisy data385 in modern ML.386
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them?387 [Yes] We have read the guidelines.388
2. If you are including theoretical results...389 (a) Did you state the full set of assumptions of all theoretical results? [Yes] We state the390 conditions for the theorems proven.391 (b) Did you include complete proofs of all theoretical results? [Yes] Proofs of all theorems392 are provided in the appendix. The relevant appendix section is noted in the corresponding393 sections containing the theorems.394
3. If you ran experiments...395 (a) Did you include the code, data, and instructions needed to reproduce the main experimental396 results (either in the supplemental material or as a URL)? [N/A] The paper is directed to397 proving a number of theorems.398
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were399 chosen)? [N/A] We do not rely on numerical experiments for the results presented in the400 paper.401
(c) Did you report error bars (e.g., with respect to the random seed after running experiments402 multiple times)? [N/A]403
(d) Did you include the total amount of compute and the type of resources used (e.g., type of404 GPUs, internal cluster, or cloud provider)? [N/A]405
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...406 (a) If your work uses existing assets, did you cite the creators? [N/A]407 (b) Did you mention the license of the assets? [N/A]408 (c) Did you include any new assets either in the supplemental material or as a URL? [N/A]409 (d) Did you discuss whether and how consent was obtained from people whose data you’re410 using/curating? [N/A]411 (e) Did you discuss whether the data you are using/curating contains personally identifiable412 information or offensive content? [N/A]413 5. If you used crowdsourcing or conducted research with human subjects...414
(a) Did you include the full text of instructions given to participants and screenshots, if applica-415 ble? [N/A] There are no human subjects.416
(b) Did you describe any potential participant risks, with links to Institutional Review Board417 (IRB) approvals, if applicable? [N/A]418
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on419 participant compensation? [N/A]420

Review summaries:
This article investigates the generalization guarantees of the Nadaraya-Watson estimator. The focus was put on the study of the finite sample risk of the estimator for two learning tasks: regression and classification. The contributions are theoretical: the excess risk at a point $x$ (in the interior of the domain) converges to $0$ at the rate $\mathcal{O}(\sigma(x)^2/\log^{\alpha}(n))$, where $\sigma(x)$ is the noise variance at the point $x$, $\alpha$ is a positive real number that depends on the learning task, and n is the number of nodes used for the interpolation. Most importantly, the convergence rate is independent of the dimension $d$. In the process, many interesting theoretical results were proven: the study of the asymptotic equivalents of the moments of the weight functions, along with a scaling limit of the associated Lagrange functions.

The paper studies convergence rates for point-wise kernel-based regression in $R^d$ with the Hilbert kernel. This kernel is parameter-free and singular, yielding an interpolator of the data. Various results are given for this estimator, in particular, it is shown that under some continuity and growth conditions on the regression function $f(x) = E[y|X=x]$ and the density $\rho$, the excess error follows the asymptotic rate of $\Theta(\sigma(x)^2/ \text{log}(n))$, where $\sigma(x)^2$ is the variance of $y$ at $x$.

This paper provides an asymptotic analysis of the Hilbert kernel interpolating estimator, at a fixed covariate point $x$. First, it is noted that, since the estimator is a Nadaraya-Watson-type estimator, the bias and variance of the estimator depend primarily on the distribution of the weights $w_i(x)$. The paper thus provides a fairly precise asymptotic characterization of the moments of the weights. This characterization is then used to bound the variance (by $O(\sigma^2(x)/\log(n)$) and the bias (by $O(1/\log(n)$). This also leads to a comparable bound on classification risk for a plug-in classifier based on Hilbert kernel regression.