Overfitting or perfect fitting? Risk bounds for classification and regression rules that interpolate

Mikhail Belkin, Daniel Hsu, Partha Mitra

Introduction

The central problem of supervised inference is to predict labels of unseen data points from a set of labeled training data. The literature on this subject is vast, ranging from classical parametric and non-parametric statistics to more recent machine learning methods, such as kernel machines , boosting , random forests , and deep neural networks . There is a wealth of theoretical analyses for these methods based on a spectrum of techniques including non-parametric estimation , capacity control such as VC-dimension or Rademacher complexity , and regularization theory . In nearly all of these results, theoretical analysis of generalization requires “what you see is what you get” setup, where prediction performance on unseen test data is close to the performance on the training data, achieved by carefully managing the bias-variance trade-off. Furthermore, it is widely accepted in the literature that interpolation has poor statistical properties and should be dismissed out-of-hand. For example, in their book on non-parametric statistics, Györfi et al. [25, page 21] say that a certain procedure “may lead to a function which interpolates the data and hence is not a reasonable estimate”.

Yet, this is not how many modern machine learning methods are used in practice. For instance, the best practice for training deep neural networks is to first perfectly fit the training data . The resulting (zero training loss) neural networks after this first step can already have good performance on test data . Similar observations about models that perfectly fit training data have been made for other machine learning methods, including boosting , random forests , and kernel machines . These methods return good classifiers even when the training data have high levels of label noise .

An important effort to show that fitting the training data exactly can under certain conditions be theoretically justified is the margins theory for boosting and other margin-based methods . However, this theory lacks explanatory power for the performance of classifiers that perfectly fit noisy labels, when it is known that no margin is present in the data . Moreover, margins theory does not apply to regression and to functions (for regression or classification) that interpolate the data in the classical sense .

In this paper, we identify the challenge of providing a rigorous understanding of generalization in machine learning models that interpolate training data. We take first steps towards such a theory by proposing and analyzing interpolating methods for classification and regression with non-trivial risk and consistency guarantees.

Many existing forms of generalization analyses face significant analytical and conceptual barriers to being able to explain the success of interpolating methods.

Existing capacity-based bounds (e.g., VC dimension, fat-shattering dimension, Rademacher complexity) for empirical risk minimization do not give useful risk bounds for functions with zero empirical risk whenever there is non-negligible label noise. This is because function classes rich enough to perfectly fit noisy training labels generally have capacity measures that grow quickly with the number of training data, at least with the existing notions of capacity . Note that since the training risk is zero for the functions of interest, the generalization bound must bound their true risk, as it equals the generalization gap (difference between the true and empirical risk). Whether such capacity-based generalization bounds exist is open for debate.

Generalization analyses based on algorithmic stability control the difference between the true risk and the training risk, assuming bounded sensitivity of an algorithm’s output to small changes in training data. Like standard uses of capacity-based bounds, these approaches are not well-suited to settings when training risk is identically zero but true risk is non-zero.

Many analyses are available for regularization approaches to statistical inverse problems, ranging from Tikhonov regularization to early stopping . To obtain a risk bound, these analyses require the regularization parameter λ\lambda (or some analogous quantity) to approach zero as the number of data nn tends to infinity. However, to get (the minimum norm) interpolation, we need λ→0\lambda\to 0 while nn is fixed, causing the bounds to diverge.

There is an extensive literature on local prediction rules in non-parametric statistics . Nearly all of these analyses require local smoothing (to explicitly balance bias and variance) and thus do not apply to interpolation. (Two exceptions are discussed below.)

Recently, Wyner et al. proposed a thought-provoking explanation for the performance of AdaBoost and random forests in the interpolation regime, based on ideas related to “self-averaging” and localization. However, a theoretical basis for these ideas is not developed in their work.

The analyses of the nearest neighbor rule and Hilbert kernel regression estimate are not based on bounding generalization gap, the difference between the true risk and the empirical risk. Rather, the true risk is analyzed directly by exploiting locality properties of the prediction rules. In particular, the prediction at a point depends primarily or entirely on the values of the function at nearby points. This inductive bias favors functions where local information in a neighborhood can be aggregated to give an accurate representation of the underlying regression function.

What we do.

Our approach to understanding the generalization properties of interpolation methods is to understand and isolate the key properties of local classification, particularly the nearest neighbor rule. First, we construct and analyze an interpolating function based on multivariate triangulation and linear interpolation on each simplex (Section 3), which results in a geometrically intuitive and theoretically tractable prediction rule. Like nearest neighbor, this method is not statistically consistent, but, unlike nearest neighbor, its asymptotic risk approaches the Bayes risk as the dimension becomes large, even when the Bayes risk is far from zero—a kind of “blessing of dimensionality”This does not remove the usual curse of dimensionality, which is similar to the standard analyses of kk-NN and other non-parametric methods.. Moreover, under an additional margin condition the difference between the Bayes risk and our classifier is exponentially small in the dimension.

A similar finding holds for regression, as the method is nearly consistent when the dimension is high.

Next, we propose a weighted & interpolated nearest neighbor (wiNN) scheme based on singular weight functions (Section 4). The resulting function is somewhat less natural than that obtained by simplicial interpolation, but like the Hilbert kernel regression estimate, the prediction rule is statistically consistent in any dimension. Interestingly, conditions on the weights to ensure consistency become less restrictive in higher dimension—another “blessing of dimensionality”. Our analysis provides the first known non-asymptotic rates of convergence to the Bayes risk for an interpolated predictor, as well as tighter bounds under margin conditions for classification. In fact, the rate achieved by wiNN regression is statistically optimal under a standard minimax settingAn earlier version of this article paper contained a bound with a worse rate of convergence based on a loose analysis. The subsequent work found that a different Nadaraya-Watson kernel regression estimate (with a singular kernel) could achieve the optimal convergence rate; this inspired us to seek a tighter analysis of our wiNN scheme..

Our results also suggest an explanation for the phenomenon of adversarial examples , which are seemingly ubiquitous in modern machine learning. In Section 5, we argue that interpolation inevitably results in adversarial examples in the presence of any amount of label noise. When these schemes are consistent or nearly consistent, the set of adversarial examples (where the interpolating classifier disagrees with the Bayes optimal) has small measure but is asymptotically dense. Our analysis is consistent with the empirical observations that such examples are difficult to find by random sampling , but are easily discovered using targeted optimization procedures, such as Projected Gradient Descent .

Finally, we discuss the difference between direct and inverse interpolation schemes; and make some connections to kernel machines, and random forests in (Section 6).

All proofs are given in Appendix A. We informally discuss some connections to graph-based semi-supervised learning in Appendix B.

Preliminaries

In this paper, we analyze two interpolating schemes, one based on triangulating and constructing the simplicial interpolant for the data, and another, based on weighted nearest neighbors with singular weight function.

2 Smoothness, margin, and regularity conditions

Below we list some standard conditions needed for further development.

For all x,x′x,x^{\prime} in the support of μ\mu,

There exist c0>0c_{0}>0 and r0>0r_{0}>0 such that

The regularity condition from Audibert and Tsybakov is not very restrictive. For example, if supp⁡(μ)=B⁡(0,1)\operatorname{supp}(\mu)=\operatorname{B}(0,1), then c0≈1/2c_{0}\approx 1/2 and r0≥1r_{0}\geq 1.

In what follows, we mostly assume uniform marginal distribution μ\mu over a certain domain. This is done for the sake of simplicity and is not an essential condition. For example, in every statement the uniform measure can be substituted (with a potential change of constants) by an arbitrary measure with density bounded from below.

Interpolating scheme based on multivariate triangulation

In this section, we describe and analyze an interpolating scheme based on multivariate triangulation. Our main interest in this scheme is in its natural geometric properties and the risk bounds for regression and classification which compare favorably to those of the original nearest neighbor rule (despite the fact that neither is statistically consistent in general).

The predictions of the plug-in classifier based on simplicial interpolation are qualitatively very different from those of the nearest neighbor rule. This is true even when restricting attention to a single simplex. Suppose, for example, that η(x)<1/2\eta(x)<1/2 for all x∈conv⁡(x1,…,xd+1)x\in\operatorname{conv}(x_{1},\dotsc,x_{d+1}), so the Bayes classifier predicts for all xx in the simplex. On the other hand, due to label noise, we may have some yi=1y_{i}=1. Suppose in fact that only yd+1=1y_{d+1}=1, while yi=0y_{i}=0 for all i=1,…,di=1,\dotsc,d. In this scenario (depicted in Figure 1 for d=2d=2), the nearest neighbor rule (erroneously) predicts 11 on a larger fraction of the simplex than the plug-in classifier based on η^\hat{\eta}. The difference can be striking in high dimensions: 1/d1/d for nearest neighbor versus 1/2d1/2^{d} for simplicial interpolation in dd-dimensional version of Figure 1. This provides an intuition why, in contrast to the nearest neighbor rule, simplicial interpolation can yield to classifiers that are nearly optimal in high dimensions.

One consequence of Proposition 3.1 for η^\hat{\eta} is that if xx is contained in two adjacent simplices (that share a <d{<}d-dimensional face), then it does not matter which simplex is used to define UT(x)U_{T}(x); the value of η^(x)\hat{\eta}(x) is the same in any case. Geometrically, we see that the restriction of the interpolating linear function to a face of the simplex coincides with the interpolating linear function constructed on a sub-simplex formed by that face. Therefore, we deduce that η^\hat{\eta} is a piecewise linear and continuous interpolation of the data (x1,y1),…,(xn,yn)(x_{1},y_{1}),\dotsc,(x_{n},y_{n}) on conv⁡(x1,…,xn)\operatorname{conv}(x_{1},\dotsc,x_{n}).

We note that our prediction rule requires only locating the vertices of the simplex containing a given point, rather than the considerably harder problem of constructing a full triangulation. In fact, locating the containing simplex in a Delaunay triangulation reduces to solving polynomial-size linear programs ; in contrast, computing the full Delaunay triangulation has complexity exponential in the (intrinsic) dimension .

2 Mean squared error

In general, each YiY_{i} may deviate from its conditional mean η(Xi)\eta(X_{i}) by a non-negligible amount, and hence any function that interpolates the training data is “fitting noise”. Nevertheless, in high dimension, the mean squared error of such a function will be quite close to that of the (optimal) conditional mean function.

3 Classification risk

We now analyze the statistical risk of the plug-in classifier based on η^\hat{\eta}, given by

We first state an easy consequence of Corollary 3.3 using known properties of plug-in classifiers.

Under the same conditions as Corollary 3.3,

When the conditional mean function satisfies a margin condition, the 1/d1/\sqrt{d} in Corollary 3.4 can be replaced with a quantity that is exponentially small in dd, as we show next.

Both Corollary 3.4 and Theorem 3.5 show that the risk of f^\hat{f} can be very close to the Bayes risk in high dimensions, thus exhibiting a certain “blessing of dimensionality". This stands in contrast to the nearest neighbor rule, whose asymptotic risk does not diminish with the dimension and is bounded by twice the Bayes risk, 2R0/1(f∗)2\mathcal{R}_{0/1}(f^{*}).

Interpolating nearest neighbor schemes

In this section, we describe a weighted nearest neighbor scheme that, like the 11-nearest neighbor rule, interpolates the training data, but is similar to the classical (unweighted) kk-nearest neighbor rule in terms of other properties, including convergence and consistency. (The classical kk-nearest neighbor rule is not generally an interpolating method except when k=1k=1.)

In what follows, we investigate the properties of interpolating schemes of this type.

We will need two key observations for the analyses of these algorithms.

The second key point is that η^(x)\hat{\eta}(x) is an interpolating scheme, provided that w(x,z)w(x,z) has a singularity when z=xz=x. Indeed, it is easily seen that if lim⁡z→xw(x,z)=∞\lim_{z\to x}w(x,z)=\infty, then lim⁡x→xiη^(x)=yi\lim_{x\to x_{i}}\hat{\eta}(x)=y_{i}. Extending η^\hat{\eta} continuously to the data points yields a weighted & interpolated nearest neighbor (wiNN) scheme.

Concretely, we will consider ϕ\phi that diverge near t=0t=0 as t↦−log⁡(t)t\mapsto-\log(t) or t↦t−δt\mapsto t^{-\delta}, δ>0\delta>0.

The denominator ∥x−x(k+1)∥\|x-x_{(k+1)}\| in the argument of ϕ\phi is not strictly necessary, but it allows for convenient normalization in view of the conditional independence of kk-nearest neighbors given x(k+1)x_{(k+1)}. Note that the weights depend on the sample and are thus data-adaptive.

Although w(x,x(i))w(x,x_{(i)}) are unbounded for singular weight functions, concentration only requires certain bounded moments. Geometrically, the volume of the region around the singularity needs to be small enough. For radial weight functions that we consider, this condition is more easily satisfied in high dimension. Indeed, the volume around the singularity becomes exponentially small in high dimension.

Our wiNN schemes are related to Nadaraya-Watson kernel regression . The use of singular kernels in the context of interpolation was originally proposed by Shepard ; they do not appear to be commonly used in machine learning and statistics, perhaps due to a view that interpolating schemes are unlikely to generalize or even be consistent; the non-adaptive Hilbert kernel regression estimate (essentially, k=nk=n and δ=d\delta=d) is the only exception we know of.

2 Mean squared error

Let η^\hat{\eta} be a wiNN scheme with singular weight function ϕ\phi. Assume the following conditions:

η\eta satisfies the (A,α)(A,\alpha)-smoothness for some A>0A{>}0 and α>0\alpha{>}0.

ϕ(t)=t−δ\phi(t)=t^{-\delta} for some 0<δ<d/20<\delta<d/2.

Let Z0\mathchar58=λ(supp⁡(μ))/λ(B⁡(0,1))Z_{0}\mathrel{\mathop{\mathchar 58\relax}}=\lambda(\operatorname{supp}(\mu))/\lambda(\operatorname{B}(0,1)), and assume n>2Z0k/(c0r0d)n>2Z_{0}k/(c_{0}r_{0}^{d}). For any x0∈supp⁡(μ)x_{0}\in\operatorname{supp}(\mu), let rk+1,n(x0)r_{k+1,n}(x_{0}) be the distance from x0x_{0} to its (k+1)(k+1)st nearest neighbor among X1,…,XnX_{1},\dotsc,X_{n}. Then

The bound in Theorem 4.3 is stated in terms of the expected distance to the (k+1)(k+1)st nearest neighbor raised to the 2α2\alpha power; this is typically bounded by O((k/n)2α/d)O((k/n)^{2\alpha/d}). Choosing k=n2α/(2α+d)k=n^{2\alpha/(2\alpha+d)} leads to a convergence rate of n−2α/(2α+d)n^{-2\alpha/(2\alpha+d)}, which is minimax optimal.

3 Classification risk

We now analyze the statistical risk of the plug-in classifier f^(x)=\mathds1{η^(x)>1/2}\hat{f}(x)=\mathds{1}_{\{\hat{\eta}(x)>1/2\}} based on η^\hat{\eta}.

As in Section 3.3, we obtain the following easy consequence of Theorem 4.3 using known properties of plug-in classifiers.

Under the same conditions as Theorem 4.3,

Choosing k=n2α/(2α+d)k=n^{2\alpha/(2\alpha+d)} leads to a convergence rate of n−α/(2α+d)n^{-\alpha/(2\alpha+d)}.

We now give a more direct analysis, largely based on that of Chaudhuri and Dasgupta for the standard kk-nearest neighbor rule, that leads to improved rates under favorable conditions.

For 0<γ<1/20<\gamma<1/2, define the effective interiors of the two classes by

Points away from the boundary, i.e., in Xp,γ−\mathcal{X}_{p,\gamma}^{-} or Xp,γ+\mathcal{X}_{p,\gamma}^{+} for p≈k/np\approx k/n, are likely to have kk nearest neighbors in Xp,γ−\mathcal{X}_{p,\gamma}^{-} or Xp,γ+\mathcal{X}_{p,\gamma}^{+}, respectively, so that interpolating their labels yields accurate predictions.

Let η^\hat{\eta} be a wiNN scheme with singular weight function ϕ\phi, and let f^\hat{f} be the corresponding plug-in classifier. Fix any 0<γ<1/20<\gamma<1/2 and p>k/np>k/n. Then

While Theorem 4.5 is quite general, the values of quantities involved can be non-trivial to express in terms of nn. The following corollary leads to explicit rates under certain conditions.

η\eta satisfies the (A,α)(A,\alpha)-smoothness and (B,β)(B,\beta)-margin conditions for some A>0A{>}0, α>0\alpha{>}0, B>0B{>}0, β≥0\beta{\geq 0}.

ϕ(t)=t−δ\phi(t)=t^{-\delta} for some 0<δ<d/20<\delta<d/2.

Let Z0\mathchar58=λ(supp⁡(μ))/λ(B⁡(0,1))Z_{0}\mathrel{\mathop{\mathchar 58\relax}}=\lambda(\operatorname{supp}(\mu))/\lambda(\operatorname{B}(0,1)), and assume

For consistency, we set k\mathchar58=n(2+β)α/((2+β)α+d)k\mathrel{\mathop{\mathchar 58\relax}}=n^{(2+\beta)\alpha/((2+\beta)\alpha+d)}, and in the bound, we plug-in p\mathchar58=2k/np\mathrel{\mathop{\mathchar 58\relax}}=2k/n and γ\mathchar58=A(Z0p/c0)α/d\gamma\mathrel{\mathop{\mathchar 58\relax}}=A(Z_{0}p/c_{0})^{\alpha/d}. This leads to a convergence rate of n−αβ/(α(2+β)+d)n^{-\alpha\beta/(\alpha(2+\beta)+d)}.

The factor 1/k1/k in the final term in Corollary 4.6 results from an application of Chebyshev inequality. Under additional moment conditions, which are satisfied for certain functions ϕ\phi (e.g., ϕ(t)=−log⁡(t)\phi(t)=-\log(t)) with better-behaved singularity at zero than t−δt^{-\delta}, it can be replaced by e−Ω(γ2k)e^{-\Omega(\gamma^{2}k)}. Additionally, while the condition ϕ(t)=t−δ\phi(t)=t^{-\delta} is convenient for analysis, it is sufficient to assume that ϕ\phi approaches infinity no faster than t−δt^{-\delta}.

Ubiquity of adversarial examples in interpolated learning

The recently observed phenomenon of adversarial examples in modern machine learning has drawn a significant degree of interest. It turns out that by introducing a small perturbation to the features of a correctly classified example (e.g., by changing an image in a visually imperceptible way or even by modifying a single pixel ) it is nearly always possible to induce neural networks to mis-classify a given input in a seemingly arbitrary and often bewildering way.

We will now discuss how our analyses, showing that Bayes optimality is compatible with interpolating the data, provide a possible mechanism for these adversarial examples to arise. Indeed, such examples are seemingly unavoidable in interpolated learning and, thus, in much of the modern practice. As we show below, any interpolating inferential procedure must have abundant adversarial examples in the presence of any amount of label noise. In particular, in consistent on nearly consistent schemes, like those considered in this paper, while the predictor agrees with the Bayes classifier on the bulk of the probability distribution, every “incorrectly labeled” training example (i.e., an example whose label is different from the output of the Bayes optimal classifier) has a small “basin of attraction” with every point in the basin misclassified by the predictor. The total probability mass of these “adversarial” basins is negligible given enough training data, so that a probability of misclassifying a randomly chosen point is low. However, assuming non-zero label noise, the union of these adversarial basins asymptotically is a dense subset of the support for the underlying probability measure and hence there are misclassified examples in every open set. This is indeed consistent with the extensive empirical evidence for neural networks. While their output is observed to be robust to random feature noise , adversarial examples turn out to be quite difficult to avoid and can be easily found by targeted optimization methods such as PCG . We conjecture that it may be a general property or perhaps a weakness of interpolating methods, as some non-interpolating local classification rules can be robust against certain forms of adversarial examples .

Let An={x∈Ω\mathchar58f^n(x)≠f∗(x)}{\cal A}_{n}=\{x\in\Omega\mathrel{\mathop{\mathchar 58\relax}}\hat{f}_{n}(x)\neq f^{*}(x)\} be the set of points at which f^n\hat{f}_{n} disagrees with the Bayes optimal classifier f∗f^{*}; in other words, An{\cal A}_{n} is the set of “adversarial examples” for f^n\hat{f}_{n}. Consistency of f^\hat{f} implies that, with probability one, lim⁡n→∞μ(An)=0\lim_{n\to\infty}\mu({\cal A}_{n})=0 or, equivalently, lim⁡n→∞∥f^n−f∗∥Lμ2=0\lim_{n\to\infty}\|\hat{f}_{n}-f^{*}\|_{L^{2}_{\mu}}=0. On the other hand, the following result shows that the sets An{\cal A}_{n} are asymptotically dense in Ω\Omega, so that there is an adversarial example arbitrarily close to any xx.

Let (X1,Y1),…,(Xn,Yn)(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n}) be the training data used to construct f^n\hat{f}_{n}. Fix a finite ϵ\epsilon-cover of Ω\Omega with respect to the Euclidean distance. Since f^n\hat{f}_{n} is interpolating and η\eta is never zero nor one, for every ii, there is a non-zero probability (over the outcome of the label YiY_{i}) that f^n(Xi)=Yi≠f∗(Xi)\hat{f}_{n}(X_{i})=Y_{i}\neq f^{*}(X_{i}); in this case, the training point XiX_{i} is an adversarial example for f^n\hat{f}_{n}. By choosing n=n(μ,ϵ,δ)n=n(\mu,\epsilon,\delta) large enough, we can ensure that with probability at least δ\delta over the random draw of the training data, every element of the cover is within distance ϵ\epsilon of at least one adversarial example, upon which every point in Ω\Omega is within distance 2ϵ2\epsilon (by triangle inequality) of the same. ∎

A similar argument for regression shows that while an interpolating η^\hat{\eta} may converge to η\eta in Lμ2L^{2}_{\mu}, it is generally impossible for it to converge in L∞L_{\infty} unless there is no label noise. An even more striking result is that for the Hilbert scheme of Devroye et al., the regression estimator almost surely does not converge at any fixed point, even for the simple case of a constant function corrupted by label noise . This means that with increasing sample size nn, at any given point xx misclassification will occur an infinite number of times with probability one. We expect similar behavior to hold for the interpolation schemes presented in this paper.

Discussion and connections

In this paper, we considered two types of algorithms, one based on simplicial interpolation and another based on interpolation by weighted nearest neighbor schemes. It may be useful to think of nearest neighbor schemes as direct methods, not requiring optimization, while our simplicial scheme is a simple example of an inverse method, using (local) matrix inversion to fit the data. Most popular machine learning methods, such as kernel machines, neural networks, and boosting, are inverse schemes. While nearest neighbor and Nadaraya-Watson methods often show adequate performance, they are rarely best-performing algorithms in practice. We conjecture that the simplicial interpolation scheme may provide insights into the properties of interpolating kernel machines and neural networks.

To provide some evidence for this line of thought, we show that in one dimension simplicial interpolation is indeed a special case of interpolating kernel machine. We will briefly sketch the argument without going into the details. Consider the space H\mathcal{H} of real-valued functions ff with the norm ∥f∥H2=∫(d ⁣⁡f/d ⁣⁡x)2+κ2f2d ⁣⁡x\|f\|_{\mathcal{H}}^{2}=\int(\operatorname{d\!}f/\operatorname{d\!}x)^{2}+\kappa^{2}f^{2}\operatorname{d\!}x. This space is a reproducing kernel Hilbert Space corresponding to the Laplace kernel e−κ∣x−z∣e^{-\kappa|x-z|}. It can be seen that as κ→0\kappa\to 0 the minimum norm interpolant f∗=arg⁡min⁡⁡f∈H,∀if(xi)=yi∥f∥Hf^{*}=\operatorname*{\arg\min}_{f\in\mathcal{H},\forall_{i}f(x_{i})=y_{i}}\|f\|_{\mathcal{H}} is simply linear interpolation between adjacent points on the line. Note that this is the same as our simplicial interpolating method.

Finally, we note that while kernel machines (which can be viewed as two-layer neural networks) are much more theoretically tractable than general neural networks, none of the current theory applies in the interpolated regime in the presence of label noise . We hope that simplicial interpolation can shed light on their properties and lead to better understanding of modern inferential methods.

Acknowledgements

We would like to thank Raef Bassily, Luis Rademacher, Sasha Rakhlin, and Yusu Wang for conversations and valuable comments. We acknowledge funding from NSF. DH acknowledges support from NSF grants CCF-1740833 and DMR-1534910. PPM acknowledges support from the Crick-Clay Professorship (CSHL) and H N Mahabala Chair (IITM). This work grew out of discussions originating at the Simons Institute for the Theory of Computing in 2017, and we thank the Institute for the hospitality. PPM and MB thank ICTS (Bangalore) for their hospitality at the 2017 workshop on Statistical Physics Methods in Machine Learning.

References

Appendix A Proofs

A.2 Proof of Theorem 3.2

Throughout we condition on X1,…,XnX_{1},\dotsc,X_{n}, and write

For the first term, observe that if X∉C^X\notin\widehat{C}, then η^(X)=1/2\hat{\eta}(X)=1/2 and hence (η^(X)−η(X))2≤1/4(\hat{\eta}(X)-\eta(X))^{2}\leq 1/4.

We now consider the second term, conditional on Z\mathchar58=(X1,…,Xn)Z\mathrel{\mathop{\mathchar 58\relax}}=(X_{1},\dotsc,X_{n}) and X∈C^X\in\widehat{C}. Let LT(X)=\mathchar58{(X(1),Y(1)),…,(X(d+1),Y(d+1))}L_{T}(X)=\mathrel{\mathop{\mathchar 58\relax}}\{(X_{(1)},Y_{(1)}),\dotsc,(X_{(d+1)},Y_{(d+1)})\}. Since X∈conv⁡(UT(X))X\in\operatorname{conv}(U_{T}(X)), its barycentric coordinates W\mathchar58=(W1,…,Wd+1)W\mathrel{\mathop{\mathchar 58\relax}}=(W_{1},\dotsc,W_{d+1}) in conv⁡(UT(X))\operatorname{conv}(U_{T}(X)) are distributed as Dirichlet⁡(1,…,1)\operatorname{Dirichlet}(1,\dotsc,1). Let ϵ(i)\mathchar58=Y(i)−η(X(i))\epsilon_{(i)}\mathrel{\mathop{\mathchar 58\relax}}=Y_{(i)}-\eta(X_{(i)}) and b(i)\mathchar58=η(X(i))−η(X)b_{(i)}\mathrel{\mathop{\mathchar 58\relax}}=\eta(X_{(i)})-\eta(X) for i=1,…,d+1i=1,\dotsc,d+1. Also, let v(x)\mathchar58=var⁡(Y∣X=x)v(x)\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{var}(Y\mid X=x) be the conditional variance function. By the smoothness assumptions, we have

by Jensen’s inequality and the bound on ∣b(i)∣|b_{(i)}|. For the second term, we have

by the bound on ∣v(X(i))−v(X)∣|v(X_{(i)})-v(X)|. Therefore

The conclusion follows by taking expectation with respect to ZZ and XX. ∎

A.3 Proof of Corollary 3.3

Recall that μ\mu is supported uniformly on a convex polytope, and that C^\widehat{C} is the convex hull of X1,…,XnX_{1},\dotsc,X_{n}. Consider the probability mass outside of C^\widehat{C}. This quantity has been intensely studied in the context of stochastic geometry [see 37, for a review]. The following exemplifies the kind of result one may expect.

Next we consider δ^T\hat{\delta}_{T}, the maximum diameter of any simplex in the triangulation TT (defined in Theorem 3.2). For many natural triangulation schemes, we expect δ^T→0\hat{\delta}_{T}\to 0 as n→∞n\to\infty. This is indeed the case with Delaunay triangulation, in which the edges of each simplex in TT are obtained by connecting the centroids of neighboring cells in the Voronoi tessellation for the given point set x1,…,xnx_{1},\dotsc,x_{n}.

Consider the Voronoi tessellation corresponding to the set x1,…,xnx_{1},\ldots,x_{n}. The Voronoi cell corresponding to xix_{i} is defined simply as {x∈C\mathchar58∀j≠i⋅∥x−xi∥≤∥x−xj∥}\{x\in C\mathrel{\mathop{\mathchar 58\relax}}\forall j\neq i\centerdot\|x-x_{i}\|\leq\|x-x_{j}\|\}, the set of points closest to xix_{i} than any other xjx_{j}. It is easy to see that each Voronoi cell is a convex set. Moreover, the distance from xix_{i} to any xx in its corresponding cell cannot exceed ϵ\epsilon, as the set x1,…,xnx_{1},\ldots,x_{n} is ϵ\epsilon-dense for CC. The edges of Delaunay triangulation connect the centroids of neighboring elements of Voronoi tessellation and thus are bounded by 2ϵ2\epsilon by the triangle inequality. The diameter of the simplex is the length of the longest edge, so the claim is proved. ∎

vanish as n→∞n\to\infty. This follows by applying Theorem A.1 and Lemma A.2. ∎

A.4 Proof of Theorem 3.5

The proof of Theorem 3.5 relies on the following tail bound.

for some absolute constants c1,c2>0c_{1},c_{2}>0 (which may depend on δ\delta but not pˉ\bar{p} nor kk).

Recall that WW has the same distribution as (G1,…,Gk)/∑i=1kGi(G_{1},\dotsc,G_{k})/\sum_{i=1}^{k}G_{i}, where G1,…,GkG_{1},\dotsc,G_{k} are independent Gamma random variables, each with unit shape and scale parameters. Let S\mathchar58=∑i=1kGi(Yi−1/2)S\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i=1}^{k}G_{i}(Y_{i}-1/2). Then, by Proposition 3.1,

Set λ∗\mathchar58=δ\lambda^{*}\mathrel{\mathop{\mathchar 58\relax}}=\delta, so we obtain

which is of the form c1pˉ⋅e−c2kc_{1}\bar{p}\cdot e^{-c_{2}k} for c1=8c_{1}=8 and c2=δ2/(4−δ2)c_{2}=\delta^{2}/(4-\delta^{2}).

Now, we prove the right tail bound for SS under the assumption that pˉ≤1/8\bar{p}\leq 1/8. Fix any y∈{0,1}ky\in\{0,1\}^{k}, and let t\mathchar58=∑i=1kyit\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i=1}^{k}y_{i} denote be the number of 11’s in yy. If t=0t=0, then we have

If 0<t<k/20<t<k/2, then by the summation property of the Gamma distribution, the conditional distribution of SS given Y=yY=y is the same as that of (Ht−Hk−t)/2(H_{t}-H_{k-t})/2, where HtH_{t} and Hk−tH_{k-t} are independent Gamma random variables with unit scale, HtH_{t} has shape parameter tt, and Hk−tH_{k-t} has shape parameter k−tk-t. The moment generating function for Ht−Hk−tH_{t}-H_{k-t} is

Since 0<t<k/20<t<k/2, the minimizer of the moment generating function is achieved at λ∗\mathchar58=1−2t/k\lambda^{*}\mathrel{\mathop{\mathchar 58\relax}}=1-2t/k. So

where RE⁡(p,q)\mathchar58=pln⁡pq+(1−p)ln⁡1−p1−q\operatorname{RE}(p,q)\mathrel{\mathop{\mathchar 58\relax}}=p\ln\frac{p}{q}+(1-p)\ln\frac{1-p}{1-q} is the binary relative entropy. Therefore, using Equation 1 and Equation 2,

where ∣Y∣\mathchar58=Y1+⋯+Yk|Y|\mathrel{\mathop{\mathchar 58\relax}}=Y_{1}+\dotsb+Y_{k}.

To put Equation 3 into the desired form, first observe that

Moreover, by a standard coupling argument, if BB is a binomial random variable with kk trials and success probability pˉ\bar{p}, then

where the second-to-last inequality uses the assumption pˉ≤1/8\bar{p}\leq 1/8, and the last inequality uses a standard Chernoff bound for binomial random variables. Therefore, combining Equation 3, Equation 4, and Equation 5,

Therefore, it is sufficient to prove our bound for the union of the simplices contained entirely within in the interior of one class. Moreover, since our bound is preserved under taking unions of sets, it is sufficient to prove the bound for the interior of a single simplex.

Let LT(X)=\mathchar58{(X(1),Y(1)),…,(X(d+1),Y(d+1))}L_{T}(X)=\mathrel{\mathop{\mathchar 58\relax}}\{(X_{(1)},Y_{(1)}),\dotsc,(X_{(d+1)},Y_{(d+1)})\} be the training examples defining one such simplex Δ\mathchar58=conv⁡(X(1),…,X(d+1))\Delta\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{conv}(X_{(1)},\dotsc,X_{(d+1)}). Without loss of generality we can assume that η(X(i))≤1/2−h\eta(X_{(i)})\leq 1/2-h for all ii (the analysis for η(X(i))≥1/2+h\eta(X_{(i)})\geq 1/2+h is the same). Conditional on X1,…,XnX_{1},\dotsc,X_{n}, the random vector XX is uniformly distributed in Δ\Delta. Therefore, the barycentric coordinates (W1,…,Wd+1)(W_{1},\dotsc,W_{d+1}) of XX within Δ\Delta are distributed as Dirichlet⁡(1,…,1)\operatorname{Dirichlet}(1,\dotsc,1). Since XX is independent of (X1,Y1),…,(Xn,Yn)(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n}), it follows that (W1,…,Wd+1)(W_{1},\dotsc,W_{d+1}) is independent of Y(1),…,Y(d+1)Y_{(1)},\dotsc,Y_{(d+1)}. Therefore, by Lemma A.3, we have

for some absolute constants c1,c2>0c_{1},c_{2}>0 (depending only on hh). Since η\eta is Lipschitz on the class interior, we have for any x∈Δx\in\Delta,

Since the above argument holds for any simplex Δ\Delta contained entirely within a class interior, we conclude

A.5 Proof of Theorem 4.5

Following Chaudhuri and Dasgupta (and in particular, the proof of their Theorem 5), we bound the probability of the event f^(X)≠f∗(X)\hat{f}(X)\neq f^{*}(X) by the sum of probabilities of three events:

the (k+1)(k+1)st nearest neighbor of XX is more than distance rp(X)r_{p}(X) from XX;

¬(E1∪E2)\neg(E_{1}\cup E_{2}) and yet f^(X)≠f∗(X)\hat{f}(X)\neq f^{*}(X).

by a multiplicative Chernoff bound; here X(k+1)X_{(k+1)} denotes the (k+1)(k+1)st nearest neighbor of x0x_{0}.

Now assume x0∉∂p,γx_{0}\notin\partial_{p,\gamma}. To bound the probability of E3E_{3}, we consider the following sampling process for (X1,Y1),…,(Xn,Yn)(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n}) relative to x0x_{0}:

Pick Xk+1X_{k+1} from the marginal distribution of the (k+1)(k+1)st nearest neighbor of x0x_{0}.

Pick kk points X1,…,XkX_{1},\dotsc,X_{k} independently from μ\mu restricted to B⁡(x0,∥x0−Xk+1∥)\operatorname{B}(x_{0},\|x_{0}-X_{k+1}\|).

For each XiX_{i}, independently pick the label YiY_{i} from the corresponding conditional distribution with mean η(Xi)\eta(X_{i}).

The distance from x0x_{0} to its (k+1)(k+1)st nearest neighbor is determined in the first step of this process, from the choice of Xk+1X_{k+1}. The kk nearest neighbors of x0x_{0} are the points X1,…,XkX_{1},\dotsc,X_{k} picked in the second step; their corresponding labels are Y1,…,YkY_{1},\dotsc,Y_{k}.

Suppose without loss of generality that x0∈Xp,γ−x_{0}\in\mathcal{X}_{p,\gamma}^{-}. It suffices to prove that, conditional on the event r\mathchar58=∥x0−Xk+1∥≤rp(x0)r\mathrel{\mathop{\mathchar 58\relax}}=\|x_{0}-X_{k+1}\|\leq r_{p}(x_{0}),

Observe that by definition of ηˉx0,r\bar{\eta}_{x_{0},r} and the assumption x0∈Xp,γ−x_{0}\in\mathcal{X}_{p,\gamma}^{-},

for each i=1,…,ki=1,\dotsc,k. Define Zi\mathchar58=ϕ(∥x0−Xi∥/r)(Yi−1/2)Z_{i}\mathrel{\mathop{\mathchar 58\relax}}=\phi(\|x_{0}-X_{i}\|/r)(Y_{i}-1/2) for i=1,…,ki=1,\dotsc,k. Since Z1,…,ZkZ_{1},\dotsc,Z_{k} are iid, the following bound holds by Chebyshev’s inequality:

The conclusion follows now from the definition of κp\kappa_{p}. ∎

A.6 Proof of Corollary 4.6

First, we bound sup⁡x∈supp⁡(μ)rp(x)\sup_{x\in\operatorname{supp}(\mu)}r_{p}(x). Let Rp\mathchar58=(Z0p/c0)1/dR_{p}\mathrel{\mathop{\mathchar 58\relax}}=(Z_{0}p/c_{0})^{1/d}, and fix any x∈supp⁡(μ)x\in\operatorname{supp}(\mu). By assumption, Rp≤r0R_{p}\leq r_{0}, so λ(B⁡(x,Rp)∩supp⁡(μ))≥c0λ(B⁡(x,Rp))\lambda(\operatorname{B}(x,R_{p})\cap\operatorname{supp}(\mu))\geq c_{0}\lambda(\operatorname{B}(x,R_{p})). Consequently,

Next, we bound μ(∂p,γ)\mu(\partial_{p,\gamma}). Pick any x0∈supp⁡(μ)x_{0}\in\operatorname{supp}(\mu). If η(x0)≥1/2+γ+ARpα\eta(x_{0})\geq 1/2+\gamma+AR_{p}^{\alpha}, then by the smoothness condition,

Hence ηˉx0,r≥1/2+γ\bar{\eta}_{x_{0},r}\geq 1/2+\gamma for all r≤Rpr\leq R_{p}. Similarly, if η(x0)≤1/2−γ−ARpα\eta(x_{0})\leq 1/2-\gamma-AR_{p}^{\alpha}, then

which implies ηˉx0,r≤1/2−γ\bar{\eta}_{x_{0},r}\leq 1/2-\gamma for all r≤Rpr\leq R_{p}. Since rp(x0)≤Rpr_{p}(x_{0})\leq R_{p} for all x0∈supp⁡(μ)x_{0}\in\operatorname{supp}(\mu), we conclude that

The claim now follows by using the (B,β)(B,\beta)-margin condition.

Finally, we bound κp\kappa_{p}. Fix any x0∈supp⁡(μ)x_{0}\in\operatorname{supp}(\mu) and 0≤r≤rp(x0)0\leq r\leq r_{p}(x_{0}) (so r≤Rp≤r0r\leq R_{p}\leq r_{0}). Then

The proof of Corollary 4.6 follows by combining the bounds on μ(∂p,γ)\mu(\partial_{p,\gamma}) and κp\kappa_{p} from Lemma A.4 with Theorem 4.5.

A.7 Proof of Theorem 4.3

Fix x0∈supp⁡(μ)x_{0}\in\operatorname{supp}(\mu). As in the proof of Theorem 4.5, we sample (X1,Y1),…,(Xn,Yn)(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n}) as follows:

Pick Xk+1X_{k+1} from the marginal distribution of the (k+1)(k+1)st nearest neighbor of x0x_{0}.

Pick kk points X1,…,XkX_{1},\dotsc,X_{k} independently from μ\mu restricted to B⁡(x0,∥x0−Xk+1∥)\operatorname{B}(x_{0},\|x_{0}-X_{k+1}\|).

For each XiX_{i}, independently pick the label YiY_{i} from the corresponding conditional distribution with mean η(Xi)\eta(X_{i}).

The distance from x0x_{0} to its (k+1)(k+1)st nearest neighbor is determined in the first step of this process, from the choice of Xk+1X_{k+1}. The kk nearest neighbors of x0x_{0} are the points X1,…,XkX_{1},\dotsc,X_{k} picked in the second step; their corresponding labels are Y1,…,YkY_{1},\dotsc,Y_{k}, so the regression estimate at x0x_{0} is η^(x0)=∑i=1kWiYi\hat{\eta}(x_{0})=\sum_{i=1}^{k}W_{i}Y_{i}, where

Define ϵi\mathchar58=Yi−η(Xi)\epsilon_{i}\mathrel{\mathop{\mathchar 58\relax}}=Y_{i}-\eta(X_{i}) and bi\mathchar58=η(Xi)−η(x0)b_{i}\mathrel{\mathop{\mathchar 58\relax}}=\eta(X_{i})-\eta(x_{0}) for i=1,…,ki=1,\dotsc,k. Then

where the last inequality follows from A.5 (below). Combining (6) and (7) concludes the proof. ∎

Fix X=x0∈supp⁡(μ)X=x_{0}\in\operatorname{supp}(\mu), and consider the sampling process for X1,…,Xk+1X_{1},\dotsc,X_{k+1} in the proof of Theorem 4.3 to define

Appendix B Interpolating kernel regression and semi-supervised learning

In this section we make some informal observations and notes on kernel regression in Reproducing Kernel Hilbert Space and semi-supervised learning. Full and rigorous exploration of these theoretically rich and practically significant topics is well beyond the scope of this paper and, likely, requires new theoretical insights. Still, we feel that some comments and observations may be of interest and could help connect our treatment of interpolation to other related ideas.

We start by discussing a simple setting already mentioned in Section 6. Consider the space H\mathcal{H} of real-valued functions ff with the norm defined as

This space is a reproducing kernel Hilbert Space corresponding to the Laplace kernel e−κ∣x−z∣e^{-\kappa|x-z|}. We can now define the minimum norm interpolant as

It is well-known that the η^\hat{\eta} can be written as a linear combination of kernel functions:

The coefficient αi\alpha_{i} can be obtained by solving a system of linear equations given by η^(xi)=yi\hat{\eta}(x_{i})=y_{i}.

It is however not necessary to solve this system to find the interpolating solution. Minimizing the norm directly, from the calculus of variations it follows that η^\hat{\eta} satisfies the following differential equation:

This equation should be solved in each interval (xi,xi+1)(x_{i},x_{i+1}) separately (here, assuming x1<⋯<xnx_{1}<\dotsb<x_{n}). The boundary conditions g(xi)=yig(x_{i})=y_{i} and g(xi+1)=yi+1g(x_{i+1})=y_{i+1} uniquely determine the solution of this second order ODE inside the interval.

Importantly, note that as κ→0\kappa\rightarrow 0, the solution tends to a linear interpolation between the sample points, since the differential equation becomes d ⁣⁡2η^d ⁣⁡x2=0\dfrac{\operatorname{d\!}{{}^{2}}\hat{\eta}}{\operatorname{d\!}{x^{2}}}=0, i.e., η^(x)=a+bx\hat{\eta}(x)=a+bx in each interval with the line passing through the samples at the ends of the interval. This observation connects RKHS interpolation in one dimension to simplicial interpolation analyzed in some detail in this paper. Higher dimensional RKHS kernel interpolation is significantly harder to analyze but some insight may be gained by considering a special case below.

B.2 Connections to semi-supervised learning

We will now discuss a discrete version of (9) on a graph and its connection to semi-supervised learning. We do not attempt any theoretical analyses of these methods here. Let G=(V,E)G=(V,E) be a (potentially weighted) graph with vertices xi,i=1,…,nx_{i},i=1,\ldots,n. Let WW be its adjacency matrix and LL the corresponding graph Laplacian. We can now consider the (finite-dimensional) space of functions ff defined on the vertices of the graph GG. The following definition of the norm is the discrete analogue of (8):

This norm defines a finite dimensional RKHS on the vertices of the graph GG. In the semi-supervised setting, where some of the vertices, x1,…,xkx_{1},\ldots,x_{k} have labels y1,…,yky_{1},\ldots,y_{k}, the interpolation problem becomes almost the same as before

Considering η^\hat{\eta} as a vector, we see that the analogue of the differential equation in (9) is the system of linear equations

The set of linear equations determining the minimum norm interpolating solution on the unlabelled points ii can be recast into a somewhat more intuitive form:

Here Ni{\cal N}_{i} is the set of neighbors of ii (i.e., nodes connected to ii by edges of the graph) and zi=∑i∈Niwijz_{i}=\sum_{i\in{\cal N}_{i}}w_{ij} is the weighted degree of the iith vertex.

The classifier for semi-supervised classification can be obtained by thresholding η^(xi)\hat{\eta}(x_{i}). This provides a graph-based interpolated semi-supervised learning algorithm similar to label propagation or interpolated graph regularization . Indeed, when κ→0\kappa\to 0, this scheme becomes label propagation. Interestingly, and consistently with the main story of this paper, it has been observed empirically in various works including the references above that interpolated semi-supervised learning typically provides optimal or near-optimal results compared to regularization.

If the graph corresponded to a (unweighted) hypercubic lattice in dd-dimensions, then the degree of each vertex is zi=2dz_{i}=2d. Thus, the interpolating solution has the property that at each unlabeled vertex, the inferred label value is proportional to the average of the assigned labels in the neighboring vertices. This is reminiscent of the interpolated nearest neighbor algorithms discussed in this paper.

While the solution of these equations generally depends on the structure of the neighborhood graph, there is a particularly simple case for which a closed form solution is easily obtained. This corresponds to the fully connected (unweighted) graph. The fully connected graph can be viewed as a local model for high-dimensional data. Similarly, it is used in the physics literature to mimic an infinite dimensional lattice.

It is easy to see with the above assumptions that the semi-supervised learning algorithm described above recovers the Bayes classifier when k→∞k\rightarrow\infty. Since each unlabeled vertex is equivalent, the solution ηiU=ηU\eta_{i}^{U}=\eta^{U} does not depend on ii. Thus, the minimum norm interpolating solution is constant on all the unlabeled points and is given by

The value of the interpolating regression function in this example is a constant and is independent of the number of unlabeled points. The plug-in classifier output is given at every site by f^=sign⁡(η^U)=sign⁡(n+−n−)\hat{f}=\operatorname{sign}(\hat{\eta}^{U})=\operatorname{sign}(n_{+}-n_{-}). Notice, as in d=1d=1, the classifier output does not depend on κ\kappa.

If p>12p>\frac{1}{2} and kk is large, then f^\hat{f} is therefore +1+1 with high probability, and for k→∞k\rightarrow\infty one recovers the Bayes classifier. Using Hoeffding’s inequality for the Binomially distributed n+n_{+}, the excess risk is exponentially small: