Local Kernel Ridge Regression for Scalable, Interpolating, Continuous Regression
Authors: Mingxuan Han
Chenglong Ye
Jeff Phillips
Venue: Accepted by TMLR ()

Abstract: We study a localized version of kernel ridge regression that can continuously, smoothly interpolate the underlying function values which are highly non-linear with observed data points. This new method can deal with the data of which (a) local density is highly uneven and (b) the function values change dramatically in certain small but unknown regions. By introducing a new rank-based interpolation scheme, the interpolated values provided by our local method continuously vary with query points. Our method is scalable by avoiding the full matrix inverse, compared with traditional kernel ridge regression.

Published in Transactions on Machine Learning Research (09/2022)
Local Kernel Ridge Regression for Scalable, Interpolating, Continuous Regression
Mingxuan Han School of Computing University of Utah
u1209601 @utah. edu
Chenglong Ye
chenglong, yee @ruky. edu
Dr. Bing Zhang Department of Statistics University of Kentucky
Jeff M. Phillips
jeffp@es utah edu
School of Computing University of Utah
Reviewed on OpenReview: https: // openreview. net/ forum? id- EDAK6F8yMM
Abstract
We study localized version of kernel ridge regression that can contimuously, smoothly interpolate the underlying function values which are highly non-linear with observed data points. This new method can deal with the data of which (a) local densily is highly uneven and (b) the function values change dramatically in certain small but unknown regions. By introducing new rank-based interpolation scheme, which can be interpreted as variable bandwidth Nadaraya- Watson Kernel Regression, the interpolated values provided by our local method can be proven to continuously vary with query points. Our method scalable by avoiding the full matrix inverse, compared with traditional kernel ridge regression.
1 Introduction
For decades the idca interpolating data has been scorned machine learning, since it was prone to overfit, and then not generalize to new data. However in the last decade, models with significantly more parameters than data points have become the norm, and these models achieve the best generalization when they have little or no error on the training data. This suggests we should revisit models that can be designed to explicitly interpolate the data.
In this paper we study (mostly low-dimensional) regression problems, but with somewhat non-standard assumptions. First, we assume the noise is very small so the model should more-or-less interpolate the training data. Second, we assume individual data points may be important for the model a local sense. The data may be strategically placed either due to generation from physical simulation, or because it was chosen adaptively fill areas that lacked coverage. Third, we assume that there is a smooth or continuous underlying model behind the data, e-g., driven by physics or some natural phenomenon. However, the underlying model itself can be complicated with many local features at various scales, which makes the exact recovery from random sampling challenging. Hence the estimate should be smooth, and at least continuous. Fourth, the data size may be reasonably large so (a) a full matrix inverse may be impractical in building the model, and (b) the model evaluation should be much faster than linear in the data size.
One approach towards dealing with these challenges is to keep the estimated model "local" that for any query q E d the estimated model only depends on nearby data points. Evaluation time should only depend these local points or simple parameters. Starting from this possible vantage, propose localized kernel

---

Published in Transactions on Machine Learning Research (09/2022)
ridge regression (local KRR) which adjusts the bandwidth/scalle locally, is very scalable, and moreover is continuous while also nearly interpolating the data.
One motivation for this modeling and new algorithm is the observation (Belkin et al.
2018b 2019a that many
modern large-scale machine learning models achieve 0 or nearly-O error on training
data. This has led to the
Hastie et al. 2020 Belkin let al. 2018b; Liang & Rakhlin 2020). A second motivation arises in modeling physics simulations (Pope
study of simpler models where interpolation can be achieved (Belkin et al. 2019b a
2013 Van Oijen & De Goey 2000) where data is generated according to well-behaved and constrained
pattern, and with very little error compared to the traditional complexity of the model and sampling of the data. Such data may be generated by simulating non-linear PDE systems with uneven local density as more data points are generated and more detail is required where the PDE evolves slowly, and fewer data points and less detail where it evolves rapidly. Models with a global scale parameter or traditional error measures (e.g., RMSE) can miss important features. Both settings have a scale where the cost. of model learning and evaluation non-trivial.
To deal with challenges posed in such datasets, our local KRR builds (nearly) interpolating regression models centered around set of model points M chosen to adapt to data density. Using fast k-nearest neighbor (KNN) search, we determine which local data X near each my E M to be included in building a local model. At any query point 9 we use a weighted average of nearby pre-built models, as determined by their model points, We introduce a new based interpolation scheme to ensure that this weighted average continuously varies with the query points, and hence the full model is continuous. As elaborated on in Nadaraya- Watson kernel regression estimator
Section 4.1 this can be interpreted as a over the local models. Our analysis guarantees the continuity of such estimator. This avoids either full KNN search on X at query time, avoids truncating continuous kernel (causing discontinuity), and avoids building a simplicial complex to allow for consistent barycentric interpolation. Local KRR is a locally-defined, nearly-iuterpolating model where the most expensive part of evaluation is a KNN query (for small k) on a reduced set of model points M C X. It is non-parametric and locally-optimized, yet efficient and scalable. We empirically demonstrate interpolation ability, generalization accuracy, and elficiency of this new local KRR model by comparing with global models, discontinuous local models, and non-learned interpolating local models. For completencss, we also identify some standard data assumptions under which can guarantee our model still achieves convergence in the standard global lg sense.
2 Background and Related Work
Kernel Ridge Regression estimator. Let (X,w) C Rd x I be a set of n observed data points with explanatory variables in X and response variable in y. We will leverage positive definite kernel k R2 x Rd + R. and will by default use the un-normalized Gaussian kernel k(æ, q) = cxp(- 1I /ll2/b?) for a bandwidth parameter >>0 0. This gives rise to a kernel matrix K E Raxa where K.j = k(a;, xr)) for Is, By E X. The kernel regression estimator at point q E RRd has the form: p(4) = Σ 0,,((zi, 4), where a is known as dual coefficients of KRR, e.g. in Murphy| (2012). The kernel ridge regression (KRR) estimator augments the standard least squares formulation with ridge term to regularize towards simpler model with a robust solution. Specifically, given ridge parameter its goal is to minimize Ei-(u(a y;)² + na? Ka, for which the optimal a can be solved as = (K + +/,,)) y.
Extensive work has shown the generalization propertics of and risk bounds for KRR (c.f. Steinwart et al. ) [Eberts& Steinwart : Cui et al.] (2021): ). In an idea related to our proposed local KRR extension, Zhang et al. (2013) designed a scalable variant of KRR that splits the data into parts and takes their average. The main distinctions with our local KRR approach are that Zhang el.al.'s a model (a) splits data randomly, not defined by local structure, (b) enforces a global bandwidth parameter so cannot adjust to local variation in data density and model complexity, (c) does not aim to or achieve interpolation.
Other work, knn-svm Hable (001)) and local-svm Meister & Steinwarl (0016). also considered kernel regression models where the data is split to form many smaller models, and these works defined the splits into spatially contiguous regions- as does our method. However these (a) local-svm combines local models into a full one by selecting a nearest one (different from our weighting scheme) and does not ensure global continuity since the models do not match at the boundary bet nearest models. For knn-svm, the model is defined
2

---

Published in Transactions on Machine Learning Research (09/2022)
implicitly, and built at query time the k nearest points so the issue persists. (b) Neither aims for, analyzes, the (near-)interpolation properties. (c) local-svm defines local regions by volume, and hence do not. refine subject to local data density.
Interpolation Learning. Recent work (Belkin et al. 2018b Bartlett et al. 2020; Liang & Rakhlin 2020
has observed that, state-of-the-art methods for large complex models (e-g., deep learning, those exhibiting double descent (Belkin et al.) 2019ah) which have more parameters than data tend to perform best on
out-of-sample data. Belkin et al. (2018b) established a theoretical foundation of generalization error analysis of such phenomenon based on the local interpolation scheme by using an unusual kernel called a singular kernel Shepard (1968). Follow on work (Belkin et al.) 2019b) showed that Nadaraya- Watson kernel regression (NWKR, a form of kernel regression that does not "learn" the a parameters) with a singular kernel has a type of guarantee of statistical optimality, Liang & Rakhlin (020) explained how using a nonlinear kernel, in the higl-dimennsioal setting (e.g., d > n); and with the ridge parameter ") goes to zero, global KRR can achieve zero training error and generalize well on new data; see also related results of Liu et al. (2021); Elkhalil et al. 2020 Karoui (2010). However, we are unaware of comparable work about generalization for KRR with ridge parameter 1 0 small, and in low moderate dimensions (n > d).
2018b) identified another method to learn an interpolating model to (X, without using kernels. This method builds a simplicial complex on X and simply interpolates the values of each yi within each simplex using zy as a vertex (e.g., using barycentric coordinates). They show for data generated by smooth manifolds, this will converge with enough samples. However, the size of simplicial complex (combinatorially describing each simplex) unfortunately grows exponentially in the dimension d. While our work is inspired by this, it devises a new interpolation scheme that does not rely the simplicial complex it only uses (fast) nearest neighbor searching. We also show similar convergence analysis is amenable to our method.
Belkin et al.
Fast Nearest Neighbor Search. The main computation cost of our algorithm will be invoking KNN search and its relatives. This is an area of computer science that has seen enormous recent progress (Arya & Mount 1993 Andoni et al.] 2015 Li et al.) 2019 Ram & Sinha 2019 [Brnhharssoon 2021 Johnson et al., 2021), and is nearing maturity in theory and practice. There are several related operations such libraries can perform, which we will employ in our algorithms: p-NNx (q) returns the p nearest neighbors of q in X; p-NNx (4) returns only the pth nearest neighbor of q in X; and RANGEX (g. returns all points in X within a distance from q. For instance, on 100-NNx(qi) queries (c.g., on SIFTIM (Jégou et al. 2011 with
n = 100, d = 128), libraries can return beyond 1000 queries a second (Bernhardsson 2021). We employ a state-of-the-art method FAISS (Johnson et al. 2021
3 Local Kernel Ridge Regression
Consider again the dataset X C RR.. Instead of building one global KRR model over the X,, we propose data density-adaptive localized version of KRR.
3.1
Constructing Local KRR
To construct the local KRR model, we first need to determine a set of model points M C X. Then build a model μ4 around cach model point my E M.
Determine M. We employ the following method to choose the set M. It is iterative, and starts with arbitrary point mi x; E X as the center of the first region. It then sets the radius of mı that it contains exactly L points. We maintain for each point z; X by how many regions it has been covered. In cach iteration, we choose an arbitrary (this is done at random) point that has been covered fewer than / times. In particular, if there are some points which have never been covered, we choose among these first. If all points have been covered least (j < t) times but some only new center among those covered only j + times (not ones covered j times). This repeated until the coverage of each point is at least times.
times (again j+ < t), then we choose a
3

---

Published in Transactions on Machine Learning Research (09/2022)
Learning the local models. When building the local KRR, for each model point m, E M, determine (-NNx(m)) and by Im; x-NNx (m:)x(m))ll- The by becomes the bandwidth for the kernel k used in KRR,
following similar approaches for variable bandwidth selection (Loftsgaarden & Quesenberry 1965 Terrell
1992). Next we retrieve the ncarby points on which to build our model is an expansion factor; by default À = 3; a common choice of where truncate a kernel (Robert, 1995p. ensures the points used to build each model sufficiently cover the neighborhood beyond 1 bandwidth, and also ensures the models provide sufficiently overlapping coverage of X for guarantees evaluation. In practice, we keep À small so that the models are constant size and run efficiently.
& Scott 1992 Fan & Gijbels
X RANGEX
Here
Now we can construct our local KRR model pú, so it evaluates a query point E R as pi(g) Ex,ex, ayk(zz. where a (K; + m,1d=1(()) K2 is the kernel matrix K restricted to X,, and y E R is the restriction of y to points in X,. We summarize the training of the local KRR model in Algorithm |
Algorithm LOCAL KRR TRAINING
1: Determine M.
2: for my E M do
3: Get x-NNx (m:); set b; = m| x-,Nx(m,)ll
4: Retrieve Xj- = RANGEX (m,, Ab;)
Learn KRR model (+(+) on (") with 6;.
end for
3.2 Evaluating Local KRR
To evaluate at query point RRd instead of using the single evaluation function (z()), use a weighted average evaluation based on the k closest model points My C. M. This only needs a fast Mg = N-NN call over the model points and the rest pro-trained. Let t 2 k: 2 2+ as we will see in Section 4.1, this will allow us to provide results on continuity the local KRR model. We generally use 1= k d+3 to be conservative. The evaluation of a query point q then a weighted average over these i models p(g) Em.EM, Pot((q) with the proportion values p; E (0, 1] and Σ p: = explained next.
KNN Model Interpolation. Given query point g E Rd, let M, m,, m2, my C M the k-nearest neighbors q in M. Set r ||g myll- Adopt the convention that rj S Tj+1- Assign weights each of the points in My as wy = (rx 77)/(1z 71). Observe that this ensures wy 0 and w1 1. Finally, we assign each model μj proportion Pj wz/(( (E wz). We summarize how to evaluate p(g) using KNN model interpolation in Algorithm 2
Algorithm LOCAL KRR EVALUATION q
Get M, {m,, m2, ... ma) (-NNM(9)
Let ry = llq m,ll for j [1.. Set weight wy Iu-ii for [1...k]
Σ (g
Return the weighted average evaluation of g p(q)=
Σ Wi
4 Properties of Local KRR
The desired properties for our regression model μ are:
1. it adapts to the local density
is efficient evaluate near data X
3. it is continuous
4. it (nearly) interpolates data points
4

---

Published in Transactions on Machine Learning Research (09/2022)
The adaption to the local density follows by how adaptively choose model points M, set the bandwidth, and will be demonstrated empirically in Section P
The efficiency is due to each model being built on a constant number of points (determined by constant parameters and A). After selecting model points, each training requires a fast RANGEX call, and then constant time training. Similarly, the evaluation requires a fast k-NN call, and then evaluating the k pre-trained models, which each require a constant (depending on E and A) amount of operations. We demonstrate its empirical scalability in Section 5
We show the contimuity and the (near))inteepollatinn properties in the following two subsections.
4.1
Continuity of Local KRR
Some technical calculation-focused proofs are deferred to Appendix A
Interpretation as NWKR estimator. We first note that weighted interpolation of the local models can be interpreted varisble-bandwiith Nadaraya- Watson kernel regression (NWKR) estimator with triangle kernel K (g, m) = max 1 |m - qll/h} with bandwidth h = r& chosen as the distance to the kth nearest neighbor model. Under NWKR, let W = =m,, ((g, m.) and the estimator at Q is 1 Em,m (1)) With this triangle kernel, only the nearest & - models have non-zero weights, which are (
W
and (via multiplying by TT/(TR - rI)) are proportional to the weights considered in our interpolation w; = 7k-ri re-r While the continuity of fixed-bandwidth NWKR is known (Ferraty et al. 2011], the following analysis can be interpreted as analyzing the continuity for this variable bandwidth condition under the triangle kernel.
Continuity of KNN Interpolation. The first step towards formalizing the contimuity of local KRR is to understand the continuity of the KNN interpolation scheme. Consider continuously moving the position of a query point q. As any model point moves into the set of k nearest of q, its weight is initially zero. So there is boundary effect discontinuity in the contribution of the nearest k models when suddenly a diifferett set of points are the i ncarest. Also, when two points change their relative position in the sorted order, they have the same weight. So again there is no discontinuity in the weight as the relative order of model points change within that nearest k.
To formalize that KNN interpolation proportions (the Py values) are continuous in choice of 9, we only neod to assume that M is in general posilion (Todd 1969;; that is, no d + points are all equidistant from any query point. We prove Lipschitz bound with respect to the movement of q, and it will depend RR = H E r1). Recall, that rj is the distance from to the jth closest model point. For k 2 d + 2, and M in general position, then Rs >
quantity
The Lipschitz factor will be smaller when (r) and Rs is large. The same is true for simplicial interpolation based on the standard Delaunay complex Aurenhammer et al. (2013). It is common when dealing with that complex to assume general position, since it is not well-defined otherwise.
Let Rg = = - rz). Recall the proportion pj for model My in the weighted average is normalized by the total weight, and is written
-Ty (Ta - r1) Dj Ts ri Ek_ (Tk r:)
Th Ts
Now we need to consider moving the query point E Id by some small amount ò E RR To formalize this consider some other point g' = = ōu, where u E Rdd and ||u||| = 1. So the perturbation is a change in Rd, but we will only care about the magnitude 6.
Lemma 4.1. Consider moving the query q a distance at most 8 < Re, of model μj, denoted Pj. changes by at most 48/Rg.
4. The proportion of model μ comprised
Now make a stronger assumption on M, so it is evenly distributed. Its density may change, but should do so in some local way. For 7 2 0 we say M is (r, diisttrbutee if any point q E Rd satisfies R > y. If a point set is (0, d + 2)-distributed, it is degenerate in the sense that d + 2 points are equidistance to some
5

---

Published in Transactions on Machine Learning Research (09/2022)
point q, and so not in general position. For k 2 d can expect y > 0. We can now understand the continuity bound Lemma 4.1 in terms of 7-
Corollary 4.2. For model points M that are (7, distribbuted consider moving the query q a distance at most 8 S /4. The proportion of model μ comprised of model μj, denoted Pj, changes by at most 5 4
7
Proof. Recall that this implies Rs 2 if M is (7, k) -istrribbd Thus R4 YKs, and 6/Rs S 6/(r)))
Let us interpret this bound. As a point q moves distance 6, then the proportion that a model can change moves at most proportional to o; it has a Lipschitz parameter __ This parameter has two parts.
TET"
This first is rk: this has the same unit and scaling as the data. All Lipschitz bounds that convert between two things with different units (e.g., q and pý) need some such term make sense. In this case, however, we note that r is local. That is, if the local density of points changes, then the change in effect from q on Júj also changes. Note that using the more generic form in Lemma 4.1 had similar effect.
The second part is 7, which is a notion of how close a point set is to degeneracy, mcasured in a density agmostic way. This value has no unils, so it does not interfere with the properties discussed about ry- Now recall that as the point set becomes more degenerate, then 7 goes towards 0. That makes the stability worse; a smaller change in q can effect a larger change in Pş- So less degenerate distributions M have more stable relationship between the movement of the query, and the change in the functions used.
Local KRR model stability. Now based on the Lemma 4.1 and Lipschitz continuity each Hy we can show that the model is is Lipschitz continuous.
Theorem 4.3. Assume M is (r, k) distributed, for any M4 E M, for q. E Rd with llq - q'll < Rx/4 and (يا((9) - py(2)) S L- 8 and py(g) (-T,T), then Ie(g") p(g) S My(--1yT
R2 + L8.
From the Theorem we can see that our model μ is Lipschitz continuous bounded by constant depending on the Lipschitz constant of each μj and also depend on the range of each μíz And a bounded range of μyy in Meister &
[-T,7] is common for kernel ridge regression (Steinwart et al., 2009 Fberts & Steinwart, 2013
Steinwart 2016). The (k- - 1) shown in the bound is typically Odd) e.g. k = +.3.
Let's compare this to the simplicial-baryeentric barycentric interpolation (Belkin et al. 2018b). Both require general position, although ours can loosen this requirement by increasing k. Theirs relies on knowing the (d 1)- simplex (e.g, triangle in R²') that query point lies in to determine the proportions. Our KNN interpolation does not, and only needs the k-nearest neighbors. As the complexity of the simplicial complex grows exponential in d, ours should be more generally applicable.
4.2 Local KRR (Nearly) Interpolates (X,y)
We now show how KRR and our proposed local KRR model nearly interpolates a data set (X, y) C Rd + R. We say a function g Rd R w-interpolates (X, y) if for all (z,, y) (X. v) we have that (g(2,) ysl S v. While a kernel regression model can O-interpolate any data set, it is unstable, and may result in extreme model, c.g., with very large gradient. However, for ridge model to nearly interpolate, need to include some assumptions on the niceness of the data, otherwise, the ridge regularization will force the learned model to favor other simplicity factors over fitting the data.
Small noise, good fit assumption. We now describe a noise property for data (X,y) Rd R, for kernel ridge regression models that leads to bounds in interpolation. It involves parameters for noise v and goodness of fit A. We say data that satisfy this have an (v, A)-fit. This assumes a known reproducing kernel (e.g., a Gaussian) with a known bandwidth b.
Specifically, (v, A)-fit data (X, y) should be so that for each yi there of so ly, - sA1 S and so there exists a model fa (q) Ex,ex āzk(z,, 9) so that fa;);) = y and so the the norm aā S A. We think of this fa potential generating model.
6

---

Published in Transactions on Machine Learning Research (09/2022)
For any reproducing kernel k, and any data where Fa + 2, for all i i j and Y E [-T.T] with T bounded, there exists some interpolating solution. However, this solution is often unstable (involving k- and so a small ridge term 17 is added, and solve for ân = arg min, Ev - fa(z))) + +na" Ka instead. By duality, there exists some decreasing function of 7 we call h(n), so the ridge solution o, is the one that minimizes Ei(h- - [x())) restricted to a Ka S h(n). Hence, for any reproducing kernel K, and data set (X,y) it is (v, A)-fit for some parameters V and A. The ones with very small and small A are those where points with nearby 24 values also have similar y/ values.
Ka
As motivation for this model, consider data generated from some real smooth phenomenon with high precision. So if E [- we may expect v/T < 10-6 is the amount of precision in the results. Similarly, if the data X is well-spread and close explanatory values (ri) have similar response values (y:), the model will not have heavy reliance on individual data points, and hence the corresponding model coefficients āz will be
bounded. For instance, when k(z,,) x) = 1, and fa built on a fairly uniform average over about elements, then ā₃ = / and e"Kä is roughly T².
A Analysis of the small noise, good fit setting. Now we analyze the accuracy of the model f Sa where á = (K +nT)= is chosen to minimize the standard sum of squared errors on (X, y) (the observed data, not on (X, with a ridge penalty
a
Error, = >'i: - (a(=)) + "arK..
For (v, A)-fit data, the potentially generating model å satisfies that (aror, S > D² + + nv² + nA. Moreover (rror,(()) I (ā)) To understand how well this interpolates, meaning how much deviation can occur on a single data point, in the worst case this sum of squared errors is due to a single data point (x), We can summarize in the following lemma.
Lemma 4.4. Consider an A)-it data set of size n. Then for the learned â = (K+m/)- -'y yy with corresponding model fa, ensures that for each (2i, vc) it nearly interpolates so 1y; fa()) S Vnv² +
Local Interpolation. Next we adapt the (v, A)-fit assumption to a model that requires a local fit, with adaptive bandwidth. With high density in X, this allows to reduce the bandwidth and fit more variation in the Yi values, relative to i, change. In sparse regions, larger bandwidths allow for smoother functions that generalizes over a greater span. The next definition is based on the common bandwidth selection method via the distance to the kth nearest neighbor.
We say data set (X, y) and associated reproducing kernel k is locally (v, A, if for any point x; E X, we can define a bandwidth as b,- = ||,, - x-NNx (zx)]l, so that the data set X, = N-NNx(x; is (v, A)-ii.
This modeling almost directly applies to our proposed local KRR algorithm. Under the local (v, A, f--fit assumption, each model μy constructed at a model point m, E X nearly interpolates the data X, = {x E X ||æ - myll S Ab,} } it is built on, with error at most /n,v² + nA where n, = [Xy]- To apply this to the full local KRR model, we need a property on the relationship between M and X.
Consider two subsets X,M C Rd parameter t, > O and real value X >> 0. For any E X let Ma. C M be the lo points in M closest to z. And non-symmetricallyy given any point m E M let Xm.eA C X be the points in X that are within a distance Xb(E) from m, where b(E) = ||m- - E-NN (m)]l- We say X, M C Rd are (A, E, ke) -ba-anaace if for any point x E X, for each m E My.k we have that x is in Xm,,.e This captures that X and M have similar densities in a key way, and that their point sets do not vary in density too abruptly. The difference in l # k allows for different sizes in |X| and [M]. In principal the ratio l/k increases with X (recall M C X), but will also depend on the intrinsic dimensionality in some way. Hence, we just leave the two parameters separate. The factor À >>1 (even after adjusting for size) accounts for some strange boundary conditions and imbalance between the data sets. We define this as radius \b(k), and not the (Ak)th closest point since the use of Λb(k) corresponds to how we build a model around each m using a kernel with bandwidth defined as b(k).
7

---

Published in Transactions on Machine Learning Research (09/2022)
Theorem 4.5. If the data X and model points M are (X, E, -bbblance, and locally (v, k)-fit, then μ built using local KRR with ridge parameter b region expansion X, and bandwidth set by k-NNN (m;), then μ will (Vrmax 2² mA) interrotle (X, y), where Amax = m [Xyl-
Proof. The -.J,,aa-alannee assumption implies that at any point E X, then l nearest neighbors in M = {mm..... .,.. (which are used to create the model l at x) have in the set X, defined by each my. By the locally (v, A, ki-fit assumption, we know each local model Pj will (v (Vnju² -i)t-eeteelollte the data in X2 including 2, via Lemma 4.4
Finally, we recognize that j is weighted average of the models associated with M (that is μ = Σj=1 for Ej= 1). Hence, if cach is at most Vn,u² + n S V Rmmax 1² + nA far from interpolating 2, then the total interpolation error is at most Σ PVV Umax D² + n = Vnmax 2² + nA.
An important implication is that local KRR guarantees to be closer to interpolation than global KRR.
4.3 l-- D-Convergence Rate of Local KRR
In this subscction, we show that local KRR- in addition to being locally adaptive, near-interpolating, and continuous also can achieve standard learning bounds. Specifically we provide generalization error bound in an & sense under some standard and mild assumptions. We provide two bounds depending on the algorithms behavior as the number of data points goes to infinity; we either consider (a) the number of local models is fixed and the number of data points per model goes to infinity, or (b) the number of local models also goes to infinity but the number of data points per model is fixed. In either case, need to establish few slandard assumptions on the models built.
(A1): We assume we use clipped KRR for each model (Steinwart et al. 2009p. That is, for a local data set (X,, (1)) assume that Ti = minfy= E }. T- = max E y that T = T3 - Ti and that if we predict the learned local model (x) returns a value outside of some [TY, Tyl where T 2 T2, Ti S T1, and Tţ Ti < 2T, then we clip the nearest value in the range TH]; namely, Ti Tī. Note that this cannot increase the empirical error, and will not affect the near-nneerpolation property the continuity property of our model.
(A2): Let B, = fi E Rd || - m,Il < }b:} be the ball on which the ith local model is Iearned, and for which its convergence properties will be well-defined. Local KRR requires k local models (recall k = (d)) is a constant), to evaluate a query. Let B, = {i E Rd HHB I I B, # 0}| 2 k} be domain within at least k such balls, and is where the local models are defined. We assumee our data is drawn from a measure a whose support is contained in Bg. Note that the way we determine M (ensuring X C B, with 2 k) in Section 3 guarantees that the support of is in B, as the number of models |M] goes to infinite with n. And in the case where M is fixed assume this has been achieved.
(A3): Assume f" Rd + R is true generating model of our data (X, v) and that f" E L2( Lo (Rd), where Tpp are standard Lebesgue spaces in Rd. As a smoothness assumption we make one of two refinements. Either (A3.1) we assume f E (o)) where (a) is Besov-like space with smoothness parameter s; see Eberts & Steinwarl (0033) for standard definitions of function spaces. Or (A3.2) we assume y/4 = f + ec, and f" satisfies a (C1, a) Hölder condition s.t., (ff - ff(()) S C -|z-== - =' We also assume the conditional variance of 6; as Var(y; X = =x) satisfies (C), az) HHölerr condition.
Next we observe that our global prediction p is the weighted average p(x) = Σ/ Po1s) of l local predictions, where each pj E [0, 1] and Σ/= Pj 1. Hence if for true function fq, each l4 satisfies El|m - f²113 < Ry for some diminishing crror rate R,, over the domain B,, then we have
Bl| f 113 = E Σ p EI( E ps - 12 = E p3 Σ 11 1 < &
= -kR,,
so the global model μ also achieves that error rate, up to a constant factor k. Thus what remains is to show each local model achieves a certain error rate; it will then generalize to the global model. For the fixed
8

---

Published in Transactions on Machine Learning Research (09/2022)
number of models, we use Eberts & Steinwart 2013 's optimal KRR analysis, and as the number of models going to infinity we adapt Belkin et al. [0188a] analysis of a linear interpolation scheme.
Fixed number of local models. When we assume the number of local models |M] is fixed, each local
model observes at least n' points, as n goes to infinity, and hence n' goes to infinity. Under assumptions (A1), (A2), and (A3.1) we can apply our settings to Theorem 3.1 and Corollary 3.2 in Eberts & Steinwart 20133) which gives the near-optimal (up to E) Eg error rate for KRR.
Lemma 4.6. For constant > 0 and any E> O, a local model μ4 observing n' points (Xt,y E R X from o, 1 salisfying (A3. 1), and learning according to (A1), would have
EGIB, - 1' C - -
Applying the observation that local KRR is weighted average of local models, and so inherit their convergence bound (up to constant k), can state.
Corollary 4.7. Consider for data (X,y) E Rd x R from o, fr, under assumptions (A1), (A2), and (A8.1), as the number of data points n goes to infinity, with fized set of local models each observing at least n' data points. Then the model learned by local KRR satisfies, for some constant C > O, and any & 0,
Erlly 11 kC (nd) (n/))zzaate 24
Number of local models goes to infinity. When the number of points in each local model is fixed, but n goes to infinity, so the number of local models goes to infinity as well, then we adapt the approach of Belkin et al. (2018a) that considered a linear interpolation scheme. They assume the strong Hölder smoothness in (A3.2), and consider what happens as 6,, the radius of each local piece, goes to 0. Their approach shows that with (A3.2) as 8, is small enough, then a constant prediction for each region is sufficient, and we can use that T bounded in (A1) show each local model μ4 close to the constant function. Their Theorem 3.2 bounds F|jn f"ll using 4 terms. Similar to how they observe in Corollary 3.3 as
oc three of those
terms go to 0 in setting: the first since by (A2) we assume the support of a is in B3, and the second and third since by snoothness (A3.2) and bounded T in (A1) that as ó, = Ab; + 0, then our prediction μi is always within T of the constant function, which is good enough. Taking the weighted average of k local models, only the variance (y: - (xx;)) remains, as follows:
Lemma 4.8. Consider for data (X,y) E TRd x R from o, f" under assumptions (A1), (A2), and (A3.2).
For 8max = À max, bi, the model u learned by local KRR satisfies
2k Ensoolys f= KC maxl
2k:
C max]
d d+2
As the number of data points and number of local models goes to infinity, 6max À Aax; d; - 0, we have
2k
E 1111- S + d
En En No 1(y (x:)))
5 Experimental Results
We examine the performance of local KRR against alternative models, and evaluate by several metries. For the examples for d we have ground truth values on a fine grid (Xa.yG), in addition to input data (X,y).
Prediction Errors: For tasks in R² it is useful to plot the error e; p(z,) Yi over (z, v: E (Xa, Vc)
fine grid, with blue as positive error, red a negative error, and white near 0 error, in Figure 2
RMSE: The root mean square error, RMSE = / (x,ex)(((xx) 9. This is measured over a fine grid
(Xa. or test data.
Worst Case Error: The worst case error, denoted foo, shows how far the model from interpolating the true data (XG,yG): for max (n((;)) ysl-
9

---

Published in Transactions on Machine Learning Research (09/2022)
Global KRR 6 0.2 1) le 5 Local KRR. & 3
NWKR & = 0.01
Local SVM
Singulor NWKI b 0.2
Global KRR & = 0.1 le-
5 Local KRR. & 1n
NWKR & 0.05
RNN KNN SVM
Singular NWKR b 0.1
laly W
Figure 1: Fit models on synthetic data d
Relative Error: Relative error =
||--)(a) 18
is relevant for our d = 3 physics simulation example. We
show maximum over held out set.
Average Curvature: This caplures the smoothness of the model for = 1. The discrete curvature at a ve
point (xi, Wi can be defined C; ((1)++ (1))))aaa where 24, vf are discrete derivativos, and x!! v. are the discrete second derivatives. The Average Curvature avgC is the average on all grid points.
Tuning Hyper-parameters. For clarity we define several variables for local KRR, but some (like À and k) do not noticeably affect the result, other than in runtime, after their constraints are met. There are two modeling parameters that one could tune: the points per local region l to control how large the local models should be, and ridge parameter 1) to control how close to interpolating the data. Except in illustrative 1d examples, we choose these via grid search on a training set using RMSE. We consider l E {10, 20, 30, 40, 50} and r) E (1e-7, le-5, 1e-3, 1, 10). For global KRR (build one single global regression model using KRR), singular kernel method, and NWKR (with Gaussian) we select bandwidth parameter b from b {0.01, 0.05, 0.1, 1, 5), and global KRR also choose 77 from {1e-7, 1e-5, 1e-3, 1, 10). In the real data experiment of methane data, we also compare knn-svm (Hable) 2013) and local-svm (Meister & Steinwart 2016) with our methods. For knn-svm, we tune k & {10, 20, 30, 40, 50} and 1/, b in the same ranges as local KRR. For local-svm we tune, on each local model, radius r E {0.05, 0.1, 0.2, 0.3, 0.4) and 17 and b as elsewhere. For Interpooero compared in 2D simulation example, we tune number of neighbors in {10, 20, 30, 40, 50). We observe in all case that the selected values are typically near the median one, and do not largely affect the evaluation in the middle of these ranges.
5.1 1D Simulation
As an illustrative warm-up, we compare for d = our method (local KRR) with singular kernel NWKR (Belkin let al.) 2019b), NWKR with Gaussian kernel, global KRR with Gaussian kernel, Local-SVM (Meister & Steinwart 2016) and KNN-SVM (Hable, 2013). The first five points match those in (Belkin et al. 2019b), and we extend more show varying density. In Figure [ the r is the ridge parameter; the b the bandwidth; the is the exponent parameter in singular kernel (Belkin al.) 2019b). We calculate the avgC for all the fitted curves shown in the Figure [ and the average curvature values from top-left to bottom-right are {12.93, 11.85} for global KRR. (12.05, 11.01} for local KRR, (171.99, 12.19) for Caussian (17.79, 18. 18.18) for local SVM and KNN SVM, and (20.07, 133.14} for singular-NWis Our local KRR interpolates the data, and is the most smooth one with smallest. avgC = 11.01.
The singular-NWKR interpolates the data, but exhibits strange artifacts of either cusps the points if bandwidth too small, or blocky gaps if bandwidth too large. We attempted to tune parameters remove these, but could not; these artifacts even more prominent in the examples in their paper (Belkin et al.
10

---

Published in Transactions on Machine Learning Research (09/2022)
True Y
KNN SVM
Local SVM
-2
-2
-2
-1
5.0
-2.5
0.0
2.5
5.0
5500
-25
0.0
2.5
5.0
-5.0
-2.5
2.5
5.0
-5.0
-2.5
25
5.0
Local
Interpolator
Simplex
Global
0
0
-2
-2
-2
-2
Figure Synthetic data. First line, L-RR: input data, true model, then prediction errors on knn-svm, local-svm. Second line, L-RR prediction errors on local KRR, RBF-Interpolator, simplex, global KRR.
2019b The second case also has large average curvature. The GaussianNWKKR eit does not interpolate the data (bottom figure) or interpolates but is similarly blocky and has high curvature (top figure). The global KRR may not interpolate the data (top figure), but in this small example can be tuned to perform reasonable bly well (bottom figure). As the number of data points n grows, it will be harder to interpolate; see Lemma 4.4 We used the following paramcters: &t 3, = le-5, b= 0.1 for knn-svm and for local-svm radius 0.5, 17 le-5, b 0.1. One can observe strange behavior in between the data points: large drops discontinuities.
5.2 2D Simulation
As we discussed in Section local KRR can adapt the local density of the data. We design a 2D experiment where I1 and i2 are explanatory variables, and a response variable is generated with undulating pattern, shown in Figure 2 The input data changes in density as the response variable requires more resolution More specifically, we place fine grid with a data point (x1.22) E [-5, 5] x [-5,5] every 0.05, which is 201 by 201 grid. The training data is selected using i.i.d. Bernoulli trials with probability equal to
1 + 22 )))
x2) max log1ol1 + log +
2
2
2
And the response variable y is defined as:
V2 2 +
V2 U2, 22
V2
V2
21 2
x22
2
sin (=11)= + sin( _20(2
=2||=2)) =2)) + 5 coos(21) sin
In Figure pleft), the blue points are training data. (x1, 22s's chosen by Bernoulli trials, and the red points are model points. Model points are denser when the training data is dense (i.e. in the upper right and bottom left corner). The second plot shows changes much faster in the top-right and bottom-right regions (the blue color means positive values of y, the red color means negative values). The data density distribution is similar to the variation of the function value in y.
We compare methods on this 2D synthetic dataset, showing the prediction errors on grids (Figure 2). and RMSE and foc error in Table 1 This experiment also shows Belkin et al. (2019h)'s 's Simplex Method and python's 's
11

---

Published in Transactions on Machine Learning Research (09/2022)
knn-svm
local-svm
local KRR
global KRR
RBF- Interpooroooator simplex method
RMSE
0.098
0.225
0.021
0.248
0.063
0.024
8oc
3.17
9.67
2.24
3.76
2.44
2.73
Table 1: RMSE and &oo error on 2D synthetic data.
Interpotroolltto which does not have continuity guarantees. To tune prper-parameters we randomly select percent of the data from training data as validation, then perform grid-search over syper-parameter space. We use the best performed rrper-parameters on that validation set evaluate the whole 201 by 201 grid. Details better than global KRR, and much better when the input data is sparse (in the middle). While global KRR and local SVM demonstrates significant error in many regions, for the interpolation- prasddd ones (local KRR, We also note the
more thorough sensitivily analysis for local KRR is in Appendix [ Local KRR performs
RBF-Interpolator, and simplex) the error is mainly concentrated near the boundaries. simplex method, even in 2d where the complex size can be linear, it about 40 times slower than local KRR.
5.3 Real Data: Slovakian Precipitation
We next compare local KRR with global KRR on benchmark GIS dataset of Slovakia precipitation Neteler & Mitasova 2008). We choose up to 16000 training points randomly beyond that size takes too long for global KRR, and is already too expensive (more than 12 hours) for the Simplex Method - retaining 180,000 for testing. Figure (left) shows visually that local KRR has less extreme prediction errors than global KRR. We also plot the RMSE and &o error a function of training size in Figure 3 right), averagod over trials, showing 1 std range. Local KRR has the least error under RMSE for all sample sizes, and always near the least error and the least at large sample sizes compared to any other method under loo. Although local-svm and knn-svm are conceptually similar, they perform worse empirically.
Prediction error by Local KRR.
ovez Training Size
1-inf error over Training Size
750
-10cel-emm
knn-eve -avrr-auessan slobal MRB
-Knn-em
mxtte-Starstann
3label
650
Prediction error by Global KRR
600
550
500
20
Figure 3: L-R: Prediction error, RMSE, and too error for Slovakian Percipitation
5.4 3D Combustion Simulation
We next apply these techniques to a physics-driven combustion simulation. The explanatory variables X measure the mass fraction and temperature, and the response variables y capture the so-called source terms required to perform physics simulations. The data gathered from figh-fi simulation of methane using the gril2 mechanism (Frenklach et al.) 2021) in Spitfire (Hansen et al. 2020), and then its dimension reduced to 3 using linear transform (Sutherland & Parente 2009). This sort of example directly motivates the model we study: the data has variable density and requires locally different resolution, the data is of high accuracy so the model should basically interpolate the data, and such a model is to be used in simulating a PDE so is required to be smooth. Figure 4 shows the simulation lines and variable density. As such we use relative error show how the model fits the observed data in a way that adapts to the local scale. We still build single set of model parameters for each local model.
12

---

Published in Transactions on Machine Learning Research (09/2022)
The methane data set has 44,000 observations. We evaluate under 67/33 train-test split under five different random seeds then, report the average and the standard deviation of the maximum relative errors (max rel. er) from five trials. We compare our local KRR with 4 other methods including global KRR, NWKR, Simplex method, singular kernel method, knn-svm and local-svm. Although all methods can achieve high r-squared value greater than 0.95, in this combustion task, the maximum relative error is the more important metric to guarantee that PDE simulations will work. From the results in the Table 2 we observe local KRR data has 1 order of magnitude less maximun relative error than the second best method, simplex method (which takes about 10x as long local KRR in this setting). Compared to other methods, the improvement in error is 2 to 8 orders of magnitude. This shows how local KRR model is efficient, continuous, can adapt to data locally, and nearly interpolate.
max rel.er
(std dev)
local KRR
5.58
(1.10)
global KRR
2.01 e6
(4.70 e5)
NWKR
(Gaussian)
7.26 e8
(1.15 e8)
6.01
NWKR (Singular)
2784
(52.2)
Simplex
66.35
(37.81)
knn-svm 4845
(828)
local-svm 719
(69.0)
Table 2: Maximum relative error for combustion data, averaged over 5 trials, and standard de- viation.
Figure 4: Data points for combustion data, with vary- ing density, color-coded by relative error (mostly near 0) from local KRR.
6 Discussion
Main novelties. The main novelty of our paper is that it efficiently learns a locally-adapted regression function that is (a) contimuous and (b) nearly interpolates the training data. Most ML work does not strive for either of these goals, but they are natural, especially in modern over-parameterizet settings. For instance, similar divide-and-conquer kernel methods are not careful to blend local models together and either have discontinuities at model boundaries, do not nearly-interpolate points near these boundaries. Other methods that meet these goals (like singular NWKR or high-degree polynomial regression) do not adapt to the local data properties, or (like building simplicial complex) are very expensive.
Nearly interpolation property The nearly interpolation goal is to have a local model precisely fit the local data, with error tolerance for each interpolation points controlled by v. Indeed the dependence on mmax in our Theorem 4.5 shows that this is signéficantly easier for locally-learned model (where n small or a constant) than a globol model. For not every data set can KRR achieve small interpolation error (v) without large ridge parameter 17, and the data parameter A precisely characterizes this.
Time complexity and high dimensions. The time complexity is as follows where TN is the (practically efficient, even in high dimensions) time for nearest neighbor search on N data points. Choosing model points takes (n + +T)) time with a priority queue. Building models takes OKnn(T, + +(²) time if on average a point is in (kk) modcls. Making a query/prediction takes time ( M + kE). As = O(d), that is lincar in dimension, and we can make the number of points in a model O(), hence these are very efficient, even in high dimensions
13

---

Published in Transactions on Machine Learning Research (09/2022)
References
Alexandr Andoni, Piotr Indyk, Thijs Laarhoven, Nly Razenshteyn, and Ludwig Schmidt. Practical and optimal Ish for angular distance. In NeurIPS, 2015.
Sunil Arya and David M. Mount. Approximate nearest neighbor queries in fixed dimensions. In Proceedings of the Fourth Annual ACM-SIAM Symposium Discrete Algorithms, SODA '93, pp. 271 280, USA, 1993. Socicty for Industrial and Applicd Mathematics. ISBN 0898713137.
Franz Aurenhammer, Rolf Klein, and Der-Tsal Lee. Voronoi diagrams and Delavnay triangulations. World Scientific Publishing Company, 2013.
Peter L. Bartlett, Philip M. Long, Gábor Lugosi, and Alexander Tsigler. Benign overfitting in linear regression. Proceedings of the National Academy of Sciences, 107(43))30063 30070, 2020. ISSN 0027-8424.
doi:
URL https://ww..pin pnas. org/ wwg/conteen/11744//0000
Mikhail Belkin, Daniel Hsu, and Partha Mitra. Overfitting or perfect fitting? risk bounds for classification and regression rules that interpolate, 2018a.
Mikhail Belkin, Daniel Hsu, and Partha P. Mitra. Overfitting or perfect fitting? risk bounds for classification and regression rules that interpolate. Tn Proceedings of the 32nd International Conference on Newral Information Processing Systems, NIPS'18, 18, pp. 2306 2317, Red Hook, NY, USA, 2018b. Curran Associates Inc.
Mikhail Belkin, Daniel J. Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machine-learning practice and the classical bias variinneeeee trade-off. Proceedings of the National Academy of Sciences, 116: 15849 - 15854, 2019a.
Mikhail Belkin, Alexander Rakhlin, and Alexandre B. Tsybakov. Does data interpolation contradict statistical optimality? In Kamalika Chaudhuri and Masashi Sugiyama (eds.), Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pp. 1611-1619. PMLR, 16-18 Apr 2019b. URL https://proceedings mlr. p v89/belkin19a.1 html
Erik Bernhardsson. Ann benchmarks. httppr//gthhuu.comm//reeee///ame
2021. Accessed:
2021-10-12.
Hugo Cui, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová. Generalization error rales in kernel regression: The crossover from the noiseless noisy regime. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 10131 10143. Curran Associates, Inc., 2021. URL https::/proceeddingg/ neurips. cc/ cc/paper/2021/ Papeee-aapprpppop
Mona Eberts and Ingo Steinwart. Optimal regression rates for svms using gaussian kernels. Electronic Journal of Statistics, 7:1 42, 2013.
Khalil Elkhalil, Abla Kammoun, Xiangliang Zhang, Mohamed-Slim Alouini, and Tareq Al-Naffouri. Risk convergence of centered kernel ridge regression with large dimensional data. IEFE Trunsactions on Signal Processing, 68:1 1588, 2020. ISSN 1941-0476. doi: 2.1109/tsp.2020.29777 URL http:
//dx. d/dx..oi..org org. 11PP/TSP.P... 2020. 2975939
Jianqing Fan and Irene Gijbels. Variable Bandwidth and Local Linear Regression Smoothers. The Annals of Statistics, 20(4):2008 2036, 1992. doi: 10. 01//4//0//1163448800 URL http::::///oi..r.rg/./. 1000.2211/a00 1176348900
Frédéric Ferraty, Ali Laksaci, Amel Tadj, and Philippe Vieu. Kernel regression with functional response. Electronic. Journal of Statistics, 5:159 1711 2011.
14

---

Published in Transactions on Machine Learning Research (09/2022)
M. Frenklach, H. Wang. C.L. Yu, M. Goldenberg, C.T. Bowman, R.K. Hanson, D.F. Davidson, E.J. Chang, G.P. Smith, D.M. Golden, W.C. Gardiner, and V. Lissianski. htttp:/ccmbustion berkeley ber.eley.....grriiiieee.e mech/ releases. html, 2021.
Robert Hable. Universal consistency of localized versions of regularized kernel methods. J. Mach. Learn. Res., 14(1):153-186, jan 2013. ISSN 1532-4435.
Mike Hansen, Elizabeth Armstrong, James Sutherland, Josh McConnell, John Hewson, and Robert Knaus. www./sanddillabb//pptttt.... Spitfire. http://////.ttaooooo/////aaaaaaa
2020.
Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan Tibshirani. Surprises in high-dimensional
ridgeless least squares interpolation, 2020.
Jeff Johnson, Matthijs Douze, and Hervé Jégou. Billion-scale similarity search with gpus. IEEE Transactions on Big Data, 7:535-547, 2021.
Herve Jégou, Matthijs Douze, and Cordelia Schmid. Product quantization for nearest neighbor search. IEEE Transactions on Pullern Analysis and Machine Intelligence, 33(1):117-125 128, 2011. doi: 10.1109/TPAMI. 2010.57.
Noureddine El Karoui. The spectrum of kernel random matrices. The Annals of Statistics, 38(1):1 50, 2010. doi: 10. 10.1214/08-A0550 AOS648. URL ftttp:://doi.rr//// 1214/08-A06648
Wen Li, Ying Zhang, Yifang Sun, Wei Wang, Mingjie Li, Wenjie Zhang, and Xuemin Lin. Approximate nearest. neighbor scarch on high dimensional data experimenssssss analyses, and improvement. IEEE Transactions Knowledge and Data Engineering, 32(S):1475-144) 1488, 2019.
Tengyuan Liang and Alexander Rakhlin. Just interpolate: Kernel "ridgeless" regression can generalize. The Anmals of Statistics, 48(3), Jun 2020. ISSN 0090-5364. doi: 10.1214/19-a0s1849 A0S1884
URL http:
/dx.doo..rr/11
Fanghui Liu, Zhenyu Liao, and Johan A. K. Suykens. Kernel regression in high dimensions: Refined analysis beyond double descent, 2021.
D. O. Loftsgaarden and C. P. Quesenberry. A Nonparametrie Estimate of a Multivariate Density Function. The Annols of Mathematical Statistics, 36(3):1049 - 1051, 1965. doi: 10. 00.1214/aams/11777000 URL
https://ddi.orz org/ 10. 1214 http///mmmmm////000////000
Mona Meister and Ingo Steinwart. Optimal learning rates for localized svms. Journal of Machine Learning Research, 17119)):1-44, 2016. URL http://jmlr.... orgg /v/4ers/117/14-022 httpp////lllloooggpppppss/////// html
Kevin P. Murphy. Machine Leorning: A Probabilistic Perspective. The MIT Press, 2012. ISBN 0262018020.
Markus Neteler and Helena Mitasova. Open source GIS: a GRASS GIS approach, volume 773. 2008. ISBN 978-0-387-35767-6, doi: -88107777778-0337776657748
Stephen B Pope. Small scales, many species and the manifold challenges of turbulent combustion. Proceedings of the Combustion Institute, 34(1):1-31, 2013.
Parikshit Ram and Kaushik Sinha. Revisiting kd-tree for nearest neighbor search. Proceedings of the 25th acm sigkdd international conference on knowledge discovery data mining, pp. 1388-1388, 2019.
Christian P. Robert. Simulation of truncated normal variables. Statistics and Computing, 122):111-15, Jun 1995. ISSN 1573-1375. doi: 10.1007 10.1007/bb0011433922
Donald Shepard. A two-dimensional interpolation function for irregularly-ppaced data. In Proceedings of the 1968 23nd ACM National Conference, ACM '68, pp. 517-524, New York, NY, USA, 1968. Association for Computing Machinery. ISBN 9781450374866. doi: 00.1145/80116..1000 URL https: https::/doll https://ddi..rg////
15

---

Published in Transactions on Machine Learning Research (09/2022)
Ingo Steinwart, Don R Hush, Clint Scovel, et al. Optimal rates for regularized least squares regression. In COLT, pp. 79-93, 2009.
James C Sutherland and Alessandro Parente. Combustion modeling using principal component analysis. Proceedings of the Combustion Institute, 32(1):1563-115) 2009.
George R. Terrell and David W. Scott. Variable Kernel Density Estimation. The Anmals of Statistics, 20(3): 1236 1265, 1992. doi: 00.1214/a05/117634477777 URL https://doi...... http:///oooor//////// org/10..//... 1tt../////////7774//55900000
J. A. Todd. Geometry and symmetry. The Mathematical Gazette, 53(355):333-33 334, 1969. doi: 10.2307/36 00007733149999 4997.
JA Van Oijen and LPH De Goey. Modelling of premixed laminar flames using flamelet-generated manifolds. Combustion Science and Technology, 161(1):113 137, 2000.
Yuchen Zhang, John Duchi, and Martin Wainwright. Divide and conquer kernel ridge regression. In Shai Shalev-Shwartz and Ingo Steinwart (eds.), Proceedings of the 26th Annual Conference on Learning Theory, volume 30 of Proceedings of Machine Learnóng Research, pp. 592 617, Princeton, NJ, USA, 12 14 Jun 2013. PMLR.
16

---

Published in Transactions on Machine Learning Research (09/2022)
A Proofs for Continuity
We provide the proof that the interpolation proportions Lipschitz continuous.
Recall, that ry is the distance from lo the jth closest model point. Recall Rr Ei= (ri ri), and let R = - (*k rz). And recall the jth weight = (r& - rr)/ - rı). And the sum of all weights is W = سن= wy The proportion Pj for model 2j in the weighted average is normalized by the total weight, and is written
Wj
Pk rj
(rk 11)-11
Tk Tj. Rx
Pj
Tk
11
10 >x=1 (Tk rr;)
Now we need consider moving the query point q E Idd by some small amount á E R. To formalize this consider some other point q' = q + 8u, where Red and ||u|| = 1. So the perturbation is change in Rd but: will only care about the magnitude 6.
Consider moving the query point q distance at most S 4. The proportion of model comprised of model μj, denoted pj, changes by at most 45 Rx.
Lemma A.1 (Lemma 4.1.
Proof. Under the ò-perturbation, each jth distance ry is updated to a value r in the range [ry - 6, + 6). So (rk - r5) changes by al most 26, except for (rk - rR) which is always 0. Note this is true of the jth distance even the sorted order of the models change. is also true of the term effecting the contribution of the model that used to be ranked j, even its ranking changes during this smaller perturbation of g.
Given the model μj is Pj proportion of model j before the change, let p/ be the proportion after the change. Then we can analyze
Tk - Tj
Tk ri +
Pj p/
> (Tk - -r;)
Ex=1 (rk - r; + 26)
Tj
Tk [3 + 28
R3 + (k 1)28"
We first get an upper bound
Tk T3
Tk r4 28
Pj p
Rb.
Ra + (k 1)28
Tkk rū R4 (k: 11k-194 11)2
Tk T3 - 28 Rs.
R R4 + (k: - 1)28
R₁ + (k 1)28 R
Rx(rk r3)) + (rk (k
1)26 (rx TT) R 28
Ry(RR + (k - 1)26)
=2 (rk (k)(( 1) + Rs.
Ri + (k 1)28)
17

---

Published in Transactions on Machine Learning Research (09/2022)
Now using that (rx- - -r) s } E (rk - -r:) - } Rs can show
(Tk: - r (k - 1) + Rs
Pj - 1/3 S 82
R& R2 + (k- 1)26)
(rk - rK)(k - 1) + Pix
< 82
Rs (Rx)
(Tk rK 1)
20
+ 26
R
Rr
k-1 Rs
S 26
+ 25
R
R
1
S 20
+ 26
iRk
RR:
(k-
1)/j + 1
k
20
<
20
20/ Ra.
Rs.
Rx
Similarly the lower bound, using that 8 S R implies 2k0 v R2/2, it can be shown as
p' D3 S TT, + 26
Tk -rj
D3
Rs 2(k 116
Fkk - -r; + 25 Rg
Th -rj R4 2(k 1)6
R& 2(k 1)6 Rs
Rx
Rs 2(k 1)6
Ra(Fk r3)) + R5
(rk- - R2 + 2(T r;) 1)8
RB(1i, 2(k 1)8)
26 (rk r,)((k 1) + R6
R& - 2(k: 1)8)
k-1 R4 + Rs
S 28
R4 2(k - 1)8)
k-1 + 1
28
Rs 2(k - 1)8
k
26
Rs 2kō
k:
< 26
46 / Rk-
R4. R2/2
18

---

Published in Transactions on Machine Learning Research (09/2022)
And then need to provide the Cipschitz result applied to the full local KRR model.
Theorem A.2 (Theorem 4.3 Assume for any P4 E M, for q, Rd with 11g d'll ò R2/4 and
(1((7) - (5(2)) & L and pg(4) E (-T,T), then (m(d) - p(2) - 46(k-1)T + Lò.
Ra
Proof. We start with upper bound,
M(d') u())= Σ ( p((((()))
Σ (0)(()(0)) + +8) -pp
Σ (o) - p))p((g) + L8
100 P))11ug()) + L8
46 RI:
(1) + L6
48(k 1)T
+ LS.
Rs
Then similarly for the lower bound we have:
mld') m(a)
(15(((())) 1.5) (3)5(())
-104, - p;)1lu;(01 - L8
45(k 1)f
Ld.
R₁
B Syper-paramettr Sensitivity Analysis of 2D Simulations
We provide a sensitivity analysis 10-fold cross-validation on training data for the 2D simulation in Section 5.2 In that 10-fold validatii,,,,nn, each time we pick 10 percent of training data as the validation set, build the model based on the rest of the training data, and record the RMSE the 10 percent predictions. We perform this procedure 10 times for random scenario and report the average RMSE
of our local KRR method in Table 5 is based on a
We observe that there is very little variation with the choice of (default is k = + 2 and continuity requires k 2 d + 1) and (default À = 3, following literature). However, & is somewhat sensitive. This the main scale parameter that controls how many data points are required to measure a distinct patch of the model. If is too small, the modeling in that patch may be too noisy, and if is too large, may miss numerous small features. This problem dependent, parameter balancing the noise the size of features and may not be surprising it benefits from some tuning.
19

---

Published in Transactions on Machine Learning Research (09/2022)
8
RMSE
k
RMSE
10
0.055
0.035
20
0.031
RMSE
0.034
30
0.056
2
0.031
0.031
40
0.087
0.031
0.032
50
134
0.032
Table 3: RMSE with fixed X 4,
Table 4: RMSE with fixed À = 4,
Table 5: RMSE with fixed = 4,
20, tuned
k = 4, tuned
= 20, tuned
20


Reviews:
[{"review_id": "8iA5IvOSuk", "reviewer": "Reviewer_Xqyg", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Novelty \n\nThe novelty of this paper is not convincing in the current submission. The proposed algorithm can be viewed as a variant of a combination of kNN regression and kernel regression, and I think such an idea has already been widely used. For example, I can easily find the R package (https://cran.r-project.org/web/packages/KernelKnn/KernelKnn.pdf). However, there is no discussion about the difference between the proposal and such related approaches in this paper. Although I think the contribution of the theoretical part of this paper is original and worth publishing, discussion of such related work is necessary.\n\n## Quality \n\n1. The proposal uses three levels of local search, selection of M, l-NN, and RANGE, and all of them are based on different search strategies. However, it is not clearly demonstrated why such three steps are required in the proposal. For example, how about just finding k-NN for each data point and applying Algorithm 2? I imagine that they may be necessary to obtain some theoretical properties or empirical effectiveness, while it is hard to see it in the current submission. Please clearly explain the validity of the proposed strategy.\n1. Theoretical analysis is strong: it is thorough and convincing. \n1. Empirical analysis is not thorough. In particular, only two simple synthetic datasets and one real-world datasets are used, which is not convincing. Since the proposal is a simple approach, it is interesting to see how it works on various types of datasets. I strongly recommend adding experiments on a number of real-world datasets (from small to large, low-dim to high-dim).\n\n## Presentation\n\n1. This paper is overall clearly written and easy to follow.\n1. Presentation of Section 4 can be improved. Although the authors carefully explain    all the necessary materials, it is a bit difficult to understand the motivation of each theoretical result. In particular, it is not clear which part is a straightforward consequence from the kNN property itself, and which part is non-trivial due to the localness of the proposal. For example, is the first half of Section 4.2 (\"Small noise, good fit assumption\" and \"Analysis of the small noise, good fit setting.\") a direct consequence from the standard kNN? If so, they can be removed. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes local kernel ridge regression, where the key idea is to perform kernel ridge regression to each local region determined based on k-NNs.\nAlthough the proposed algorithm itself is simple, the authors thoroughly study its theoretical properties in terms of continuity, interpolation, and convergence.\nIn addition, the proposal is empirically examined on synthetic and real-world datasets compared to relevant regression methods.\n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Please refer to the Novelty, Quality 1, 3, and Presentation 2 in the above comment.", "broader_impact_concerns": "I do not have any concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "UeHt63PnBN", "reviewer": "Reviewer_j594", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strong aspects\nOne strong aspect of the paper would be the theoretical part showing that Local KRR has some good properties.\nThe continuity of Local KRR will be of particular importance in practice because some neighbor-based local models, such as nearest neighbor methods, are not continuous.\n\nAnother strong aspect of the paper would be the adaptivity and computational scalability of Local KRR.\nDifferent from global KRR and NWKR, Local KRR does not require accessing the entire training set for each individual KRR both in the training and prediction phases: each KRR handles only the points in small neighborhoods.\nThis enables Local KRR to interpolate functions with varying smoothness.\nStandard global KRR and NWKR have a difficulty of handling such functions because the bandwidth of the kernel function determines the global smoothness of the fitting function.\nThe computational scalability is also of importance in practice.\nNWKR requires computing distances between the input $x$ and all the training points which takes $O(N)$ time, while global KRR requires $O(N^3)$ time for computing the inverse of the Gram matrix.\nIn Local KRR, only the neighbor search takes time dependent on $N$ (but only $o(N)$ time), and the training and predictions of each KRR takes time dependent only on the size of the local neighborhood, which is typically far smaller than $N$.\n\n\n## Weak aspects\nThis paper has two major weakness, as follows.\n\n### Weakness 1. \nIn the 2D simulation experiment, the proposed Local KRR performs worse than `RBFInterpolator` in `scipy.interpolate`.\nThis raises a question on the practical utility of Local KRR.\n\n[Detail]\n\nI used the `neighbor=10` and the default values for the other options in `RBFInterpolator`.\nWith the option `neighbor=10`, `RBFInterpolator` fits KRR for each $x$ using the 10 neighbors of the input $x$, similar to the proposed Local KRR.\nIt attained RMSE around 0.09, far smaller than the reported 0.747 of Local KRR.\nOf course, in theory, `RBFInterpolator` with `neighbor` option is not continuous at the boundary where the membership of the neighbor changes.\nHowever, given the far smaller RMSE, we will not be able to ignore `RBFInterpolator` just because it is not continuous.\nIn practice, accurate prediction may be more worthy than the continuity.\n\n```\n### This code produces RMSE = 0.09020604653847476 ###\n\nimport numpy as np\nfrom scipy import interpolate\n\ndef prob(x):\n    p = np.log10(1 + np.log10(1 + 0.5 * (0.5 * x[:, 0] + 0.5 * x[:, 1])**4))\n    return np.maximum(0.05, p)\n\ndef response(x):\n    z1 = x[:, 0] / np.sqrt(2) + x[:, 1] / np.sqrt(2)\n    z2 = - x[:, 0] / np.sqrt(2) + x[:, 1] / np.sqrt(2)\n    return np.sin(2 * z1 * np.abs(z1)) + np.sin(0.5 * z2 * np.abs(z2)) + 5 * np.cos(z1) * np.sin(z2)\n\n# data generation\nnp.random.seed(0)\nt = np.linspace(-5, 5, 201)\nX1, X2 = np.meshgrid(t, t)\nx_mesh = np.stack([X1.flatten(), X2.flatten()], axis=1)\ny_mesh = response(x_mesh)\np = prob(x_mesh)\nq = np.random.rand(p.size)\nx, y = x_mesh[p >= q], y_mesh[p >= q]\n\n# training and prediction\nrbf = interpolate.RBFInterpolator(x, y, neighbors=10)\nz_mesh = rbf(x_mesh)\n\n# RMSE\nprint(np.sqrt(np.mean((y_mesh - z_mesh)**2)))\n```\n\n### Weakness 2. \nLocal KRR is actually a interpolation of multiple KRRs using the Nadaraya-Watson estimator.\nThis fact is completely overlooked in the paper.\n\n[Detail]\n\nRecall that Local KRR makes a prediction for the input $x$ as\n$$\ny = \\frac{\\sum_{i=1}^k w_i \\mu_i(x)}{\\sum_{i=1}^k w_i} = \\sum_{i=1}^k p_i \\mu_i(x), \\quad \\mathrm{where} \\; p_i = \\frac{w_i}{\\sum_{i=1}^k w_i}, \\quad w_i = \\frac{r_k - r_i}{r_k - r_1}, \\quad r_i = i\\mathrm{-th} \\min_{m \\in M}\\\\|x - m\\\\| .\n$$\nThe choice of the weight $w_i$ is identical to the triangular kernel with the bandwidth equals to $r_k$.\nThe triangular kernel is given as\n$$\nK(\\\\|m - x\\\\|) = \\max\\\\{0, 1 - \\\\|m - x\\\\|\\\\} .\n$$\nThe scaled version with the bandwidth $h > 0$ is given as\n$$\n\\frac{1}{h} K\\left(\\frac{\\\\|m - x\\\\|}{h}\\right) = \\frac{1}{h^2}\\max\\\\{0, h - \\\\|m - x\\\\|\\\\} .\n$$\nUsing this triangular kernel, the interpolation of multiple KRRs using the Nadaraya-Watson estimator can be then expressed as\n$$\ny = \\frac{\\sum_{m \\in M}\\max\\\\{0, h - \\\\|m - x\\\\| \\\\} \\mu_m(x)}{\\sum_{m \\in M}\\max\\\\{0, h - \\\\|m - x\\\\|\\\\}}.\n$$\nSetting $h = r_k$ yields Local KRR because $\\max\\\\{0, r_k - \\\\|m - x\\\\|\\\\}$ is nonnegative only for $\\\\|m - x\\\\| \\ge r_k$, i.e, the first $k$ nearest neighbor of $x$.\n\nThis finding opens up several questions, as follows.\n\nFirst, what is the advantage of using the triangular kernel?\nOther kernels, such as Epanechnikov kernel and Quartic kernel, can also induce similar estimators using only the $k$ nearest neighbors.\nThat is, the proposed weight $w_i$ is not a unique choice for the purpose considered in the paper.\n\nSecond, some of the theoretical findings can be explained by the Nadaraya-Watson estimator.\nFor example, the continuity of the fitted function (one of the major contributions) would be trivial from the continuity of the Nadaraya-Watson estimator.\nThe interpolation property as well as the convergence rate reported in the paper can be explained as well.\nFor example, [Ref1] provides an asymptotic convergence of the Nadaraya-Watson estimator for interpolating multiple functions.\n\n[Ref1] Kernel regression with functional response, Electric Journal of Statistics, 2011.\n\n### Minor Weakness.\nI found that Theorem 4.3 can be improved by omitting the factor $L \\delta$.\n$$\n| \\mu(q') - \\mu(q) | = \\left| \\sum_{i=1}^{k-1} (p'_i \\mu'_i(q') - p_i \\mu_i(q)) \\right| \n$$\n\n$$\n\\le \\left| \\sum_{i=1}^{k-1} (p'_i T - p_i T) \\right|\n$$\n\n$$\n\\le T \\sum_{i=1}^{k-1} \\left| p'_i - p_i \\right| \n$$\n\n$$\n\\le \\frac{4 \\delta (k - 1)T}{\\bar{R}_k} .\n$$\nFor this bound to hold, only the Lipschitz continuity of $p$ and the boundedness of $\\mu$ will be sufficient.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "### Research Target\nThis paper considers a supervised regression with a special focus on the data interpolation.\n\n### Overview of the method\nThis paper proposes a method called Local KRR, as the smoothed weighted average of multiple KRRs.\nLet $D = \\\\{x_n, y_n\\\\}\\_{n=1}^N$ be a training set.\nLocal KRR proceeds as follows.\nFirst, we select a subset of input points $M \\\\subset \\\\{x_n\\\\}\\_{n=1}^N$ so that all the inputs reside in the close neighborhoods of $M$.\nSecond, for each $m_j \\in M$, we fit KRR $\\mu_j$ using the data points in $D$ within the close neighborhood of $M$ so that $\\\\|m_j - x_n\\\\| \\le \\tau_m$ for some $\\tau_m > 0$ (set as $3 \\times$ the distance to the $\\ell$-th nearest neighbor in the paper).\nThird, at the prediction phase for the input $x$, the predictions made by each KRR $\\mu_j(x)$ is averaged using a distance-based weight as follows.\nLet $m_i \\in M$ be the $i$-th nearest neighbor of the input $x$, and denote its distance to $x$ by $r_i = \\\\|m_i - x\\\\|$.\nThe weight for KRR $\\mu_i$, corresponding to $m_i$, is chosen as $w_i = \\frac{r_k - r_i}{r_k - r_1}$, where $k$ is the number of neighboring KRRs used for prediction.\nThe resulting predictions is then given as $\\frac{\\sum_{i=1}^k w_i \\mu_i(x)}{\\sum_{i=1}^k w_i}$.\n\n### Some properties of Local KRR\nIn the paper, the authors reported that Local KRR has the following four properties.\n1. it adapts to the local density\n2. it is efficient to evaluate near data X\n3. it is continuous\n4. it (nearly) interpolates data points\n\nThe first property follows from the fact that each KRR is trained using data points in a small neighborhood.\nSuch a local model can adapt to the smoothness of the underlying function within the small neighborhood.\nSome local models get very smooth if the underlying function is almost constant in the small neighborhood, while some other local models may become sharper if the underlying function changes dynamically in the small neighborhood.\nNote that, it would be difficult to fit a function with both smooth and sharp local behaviors using a single model such as KRR and Nadaraya-Watson Kernel Regression (NWKR) with a kernel function with a common bandwidth.\n\nThe second property follows from the fact that the prediction of Local KRR is made by the weighted average of only $k$ KRRs.\nBecause each KRR is trained using a small number of points, the prediction made by each KRR is computationally fast, in contrast to global KRR that uses all the training points for prediction.\n\nThe third and fourth properties follows from the fact that the weight $w_i$ is continuous / Lipschitz continuous in the neighborhood of any input $x$.\nThe authors provided rigorous proofs for these properties.\n\n### Experimental results\nThe authors performed experiments on 1D, 2D simulation data, and 2D and 3D real world data (Slovakian Precipitation, Combustion Simulation).\nThe results show that Local KRR performs better than some standard interpolation methods such as NWKR, KRR, and some of the neighbor-based local models.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "### Request 1.\nAs I pointed out in Weakness 1, the performance of Local KRR seems to be inferior to `RBFInterpolator` of `scipy.interpolate`.\nThis raises a question of the effectiveness of Local KRR.\nThe selection of the baseline methods as well as their hyperparameter choices may require reconsideratoin.\n\n### Request 2.\nAs I pointed out in Weakness 2, Local KRR can be interpreted as the interpolation of multiple KRRs using the Nadaraya-Watson estimator.\nBecause there are rich literatures on the Nadaraya-Watson estimator, exhaustive discussions will be essential, in particular, on the relationship with the prior studies in the literatures, and on the novelty of Local KRR over the existing literatures.\nThe current paper completely misses this relationship.", "broader_impact_concerns": "There is no ethical concern.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "PF94pEkmTj", "reviewer": "Reviewer_DPQ8", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n1. **Motivation**: The authors claimed the importance of studying interpolation from the perspective of recent advances in overparametrization and physical modeling. The introduction is well-written and motivating.\n2. **Simplicity**: The proposed local model is simple in terms of practical implementation, which requires only choosing the center points of the local models fitted independently. In addition, simplicity leads us to understand the model behavior in a transparent way. Indeed, the analysis provided in Section 4 can be followed rather effortlessly.\n3. **Continuity of the model**: Surprisingly, many localized kernel models do not exhibit a good continuity property and curvature. In contrast, the proposed model has nice Lipschitz continuity (with respect to a query point), which is practically important as well, particularly in physical modeling.\n\n## Weaknesses\n\n1. **Dependency on data dimension**: Although the authors stated that the focus of this paper is not extremely high-dimensional data, I cannot avoid mentioning that the proposed model requires computational time linearly proportional to the data dimension due to $k = \\Theta(d)$ points chosen in the query phase. Perhaps, it is interesting to see the empirical behavior of the proposed model with $k \\ll O(d)$, regardless of whether the chosen points $M$ are in general position or not.\n2. **Many hyperparameters**: The proposed model increases the number of hyperparameters compared with the global kernel ridge regression. Specifically, the number of coverages $t$ (used when we choose the center points $M$), the number of nearest neighbors $\\ell$ (used when we adaptively choose the kernel bandwidth), the expanding factor $\\lambda$ (used when we form the base set on which a single local model is applied), and the number of nearest neighbors $k$ (used when we query a test point) are newly introduced. Whereas the authors claimed that most of the hyperparameters do not influence the empirical performance so much, I would recommend supporting this point by showing the insensitivity with some experimental results.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In recent machine learning studies, the interpolation regime has gained more attention after seeing the well-known double descent phenomenon.\nThe interpolation is also interesting from the perspective of modeling physical dynamics that may not admit a nice closed-form formula.\nDespite the classical machine learning paradigm paying more attention to the generalization perspective, this paper seeks a better interpolator with a localized version of kernel ridge regression.\nThe proposed model integrates multiple local kernel ridge regressors that change the kernel bandwidth adaptively to the local data density.\nThis approach is effective compared with the global model, particularly when the number of total data grows.\nThe local kernel ridge regression also has nice properties such as continuity (with respect to a query point) and guarantee on the (pointwise) interpolation error.\nLastly, the authors provide extensive experiments to see that the local kernel ridge regression indeed has better curvature and interpolation error globally on synthetic and real datasets.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "## Major requests\n\n- **Hyperparameter sensitivity**: As mentioned in Weaknesses 2, having experiments to test hyperparameter sensitivity would be nice.\n- **On general position (1)**: Some points could be made clearer. First, you can state its definition in Section 3.2 (when you first mention general position) or earlier; the definition is firstly stated in Section 4.1 in the current version, making the validation of the statement in Section 3.2 slightly time-consuming.\n- **On general position (2)**: Second, its definition stated in Section 4.1 could be corrected: \"no $d+2$ points are all equidistant from any query point\" -> \"given $d+1$ points, no point are equidistant from any other $d$ points\". If I understand correctly, the latter corresponds to one of the standard definition of general position (no $d+1$ points are contained in any $d$-dimensional hyperplane). Let's consider the case $d=2$, then three ($d+1$) points contained in a straight line are not in general position (let's say this is *[condition A]*). Indeed, we can choose three points contained in a straight line such that the \"middle\" point is equidistant from the other two points. If we only preclude the case where four ($d+2$) points are not equidistant from each other, [condition A] is not necessarily satisfied. Hence, your definition does not seem to recover the standard one. Please correct me if I misunderstand.\n- **On general position (3)**: Lastly, validation of the fact that $t \\ge k \\ge d+2$ implies $M_q$ is unique in Section 3.2 requires a little bit of heavy work (at least for me). The inequality $t \\ge k$ was okay, but especially $k \\ge d+2$ was a bit heavy (partly because the definition could be slightly incorrect as mentioned above). This part could be explained a little bit more.\n- **Some explanation in Section 4.1**: The first paragraph of \"Continuity of KNN Interpolation\" was ambiguous. I would like to understand more what you refer to by \"model point\" and \"the set of $k$ nearest neighbors\" (neighbors of which point?), \"boundary effect\", and \"jumping effect\". Explaining what they refer to must make this paragraph easier to follow.\n\n## Minor requests\n\n- `\\citet` and `\\citep` can be used more appropriately. When you do not use a citation as a noun, `\\citep` is a more appropriate choice. For example, citations in the fourth paragraph of Section 1, the first paragraph of \"Learning the local models\" in Section 3.1, the second paragraph of \"Continuity of KNN Interpolation\" in Section 4.1, the first paragraph of Section 5.1, and throughout Section 5.2 to Section 5.4 are more natural to be written with `\\citep`.\n- The sixth line of \"Kernel Ridge Regression estimator\" in Section 2: \"dual coefficient\" -> \"dual coefficients\"\n- In \"Interpolation Learning\" of Section 2, what the existing work has done is usually written in the past tense. For example, \"Belkin et al. (2018b) establish a theoretical foundation\" -> \"established\", \"Follow on work (Belkin et al., 2019b) shows that\" -> \"Follow-up work (...) showed that\", \"Liang & Rakhlin (2020) explains how using\" -> \"explained\" in the first paragraph.\n- In the last line of the first paragraph of \"Kernel Ridge Regression estimation\" in Section 2, do you mean that there does not exist a prior work about generalization for KRR with *moderately large ridge parameter* (namely $\\eta$ is moderately large than zero), instead of $\\eta = 0$ or small?\n- The sixth line of \"Fast Nearest Neighbor Search\" of Section 2: \"$p\\_{th}$\" -> \"$p$th\", to make it consistent with other parts of the paper.\n- The second line of \"Determine $M$\" in Section 3.1: $c\\_1$ -> $m\\_1$? Also in the same line, do you consider the $\\ell\\_2$ hypersphere when you cover $\\ell$ points around the chosen point $m\\_1$? (please make it explicit)\n- In the third paragraph of \"Continuity of KNN Interpolation\": \"The Lipschitz factor will be larger\" -> \"will be smaller\"? I'm looking at the statement of Theorem 4.3.\n- In the third paragraph of \"Continuity of KNN Interpolation\", do you need to mention \"the standard Delaunay complex\"? Since many readers do not have background knowledge of it, this explanation may be distracting. If you feel necessary, I would recommend you add more explanation.\n- Right after the proof of Corollary 4.2: \"Lets\" -> \"Let us\"\n- In Theorem 4.3: Please make it explicit that the theorem assumes $M$ is $(\\gamma,k)$-distributed.\n- In the second paragraph of \"Local Interpolation\" in Section 4.2: $b\\_i = \\ell\\mathrm{NN}(x\\_i)$ -> $b\\_i = \\\\|x\\_i - \\ell\\mathrm{NN}(x\\_i)\\\\|$\n- In the ninth line of the fourth paragraph of \"Local Interpolation\" in Section 4.2, what do you refer to by \"this\" in \"We define this as a radius...\"?\n- You do not seem to have defined \"global KRR\", which becomes a little bit confusing in the latter part of the paper because there are several confusing terminologies such as \"global KRR\", \"local KRR\", and \"full local KRR\".\n- In Section 4.3, after \"Fixed number of local models\", the notation of expectations should be corrected to $\\\\mathbb{E}$ occasionally. In Lemma 4.6, the real value set should be denoted by $\\\\mathbb{R}$.\n- In Figure 1, I prefer to see the name \"singular NWKR\" in the captions of the right-most column.", "broader_impact_concerns": "This section does not apply to the current submission.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a new local kernel ridge regression method and analyzes its theoretical properties. The proposed method is simple yet effective, and the theoretical analysis of the adaptive KRR method is a good contribution. Moreover, the authors successfully address some of the concerns raised by the reviewers.

However, as reviewer j594 pointed out, the adaptive kernel width is relatively standard in practice; the main contribution is its theoretical analysis. Moreover, the hyperparameter sensitivity could be addressed more in the camera-ready version.  Therefore, we ask the authors to address the following issues in the camera-ready manuscript. 

Main change:
1. The adaptive kernel width learning is a standard technique in practice, as pointed out by reviewer j594; claiming to propose an adaptive kernel method in this paper is too strong. As I wrote above, the key contribution would be the theoretical analysis of the adaptive kernel width method. So, please rewrite the introduction to emphasize the theoretical analysis of a Nadaraya-Watson-type method (not proposing adaptive KRR itself).

2. More careful evaluation of the hyperparameter selection raised by Reviewer j594. "The proposed method has a performance comparable to RBFInterpolator of scipy.interpolate, if we carefully choose several hyperparameters". In the revised paper, please clearly describe the hyper-parameter tuning procedure (90 percent of data is treated as training data, and 10 percent of data as test data?). Also, please run cross-validation or hold-out validation for the illustrative case and report the performance.

3. Related to Question 2. If the performance is comparable to scipy.interpolate in the illustrative example, please provide another example that the proposed method works, and discuss the reason carefully. If the performance of the proposed method is better than the global one in Question 2, the authors do not need to address this request.

  
Minor change:
1. Please make the font size of Figures 1 and 3 (left panels) bigger.
2. The font of Figure 2 is not consistent. Please use the same font for all figures.

Decision: Accept (poster)