Supervised Kernel Thinning

Authors: albert  gong
kyuseong  choi
raaz  dwivedi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:38:07

Abstract: The kernel thinning algorithm of Dwivedi & Mackey (2024) provides a better-than-i.i.d. compression of a generic set of points. By generating high-fidelity coresets of size significantly smaller than the input points, KT is known to speed up unsupervised tasks like Monte Carlo integration, uncertainty quantification, and non-parametric hypothesis testing, with minimal loss in statistical accuracy. In this work, we generalize the KT algorithm to speed up supervised learning problems involving kernel methods. 
Specifically, we combine two classical algorithms---Nadaraya-Watson (NW) regression or kernel smoothing, and kernel ridge regression (KRR)---with KT to provide a quadratic speed-up in both training and inference times. We show how distribution compression with KT in each setting reduces to constructing an appropriate kernel, and introduce the Kernel-Thinned NW and Kernel-Thinned KRR estimators. We prove that KT-based regression estimators enjoy significantly superior computational efficiency over the full-data estimators and improved statistical efficiency over i.i.d. subsampling of the training data. En route, we also provide a novel multiplicative error guarantee for compressing with KT.  We validate our design choices with both simulations and real data experiments.

Main paper: 
 introduction: 
 In supervised learning, the goal of coreset methods is to find a representative set of points on which to perform model training and inference. On the other hand, coreset methods in unsupervised learning have the goal of finding a representative set of points, which can then be utilized for a broad class of downstream tasks-from integration [10, 9] to non-parametric hypothesis testing [8] . This work aims to bridge these two research threads.Leveraging recent advancements from compression in the unsupervised setting, we tackle the problem of non-parametric regression (formally defined in Sec. 2). Given a dataset of n i.i.d. samples, (x i , y i ) n i=1 , we want to learn a function f such that f (x i ) ≈ y i . The set of allowable functions is determined by the kernel function, which is a powerful building block for capturing complex, non-linear relationships. Due to its powerful performance in practice and closed-form analysis, non-parametric regression methods based on kernels (a.k.a "kernel methods") have become a popular choice for a wide range of supervised learning tasks [13, 19, 23] .There are two popular approaches to non-parametric kernel regression. First, perhaps a more classical approach, is kernel smoothing, also referred to as Nadaraya-Watson (NW) regression. The NW estimator at a point x is effectively a smoothing of labels y i such that x i is close to x. These weights are computed using the kernel function (see Sec. 2 for formal definitions). Importantly, the NW estimator takes Θ(n) pre-processing time (to simply store the data) and Θ(n) inference time for each test point x (n kernel evaluations and n simple operations).Another popular approach is kernel ridge regression (KRR), which solves a non-parametric least squares subject to the regression function lying in the reproducing kernel Hilbert space (RKHS) of a specified reproducing kernel function. Remarkably, KRR admits a closed-form solution via inverting the associated kernel matrix, and takes O(n 3 ) training time and Θ(n) inference time for each test point x.Our goal is to overcome the computational bottlenecks of kernel methods, while retaining their favorable statistical properties. Previous attempts at using coreset methods include the work of Boutsidis et al. [4] , Zheng and Phillips [30] , Phillips [17] , which depend on a projection type compression, having similar spirit to the celebrated Johnson-Lindenstrauss lemma, a metric preserving projection result. So accuracy and running depend unfavorably on the desired statistical error rate. Kpotufe [14] propose an algorithm to reduce the query time of the NW estimator to O(log n), but the algorithm requires super-linear preprocessing time.Other lines of work exploit the structure of kernels more directly, especially in the KRR literature. A slew of techniques from numerical analysis have been developed, including work on Nyström subsampling by El Alaoui and Mahoney [11] , Avron et al. [1] , Díaz et al. [7] . Camoriano et al. [5] and Rudi et al. [21] combine early stopping with Nyström subsampling. Though more distant from our approach, we also note the approach of Rahimi and Recht [20] using random features, Zhang et al. [29] using Divide-and-Conquer, and Tu et al. [27] using block coordinate descent.Our contributions. In this work, we show how coreset methods can be used to speed up both training and inference in non-parametric regression for a large class of function classes/kernels. At the heart of these algorithms is a general procedure called kernel thinning [10, 9] , which provides a worst-case bound on integration error (suited for problems in the original context of unsupervised learning and MCMC simulations). In Sec. 3, we introduce a meta-algorithm that recovers our two thinned non-parametric regression methods each based on NW and KRR. We introduce the kernel-thinned Nadaraya-Watson estimator (KT-NW) and the kernel-thinned kernel ridge regression estimator (KT-KRR).We show that KT-NW requires O(n log 3 n) time during training and O( √ n) time at inference, while achieving a mean square error (MSE) rate of n -β β+d (Thm. 1)-a strict improvement over uniform subsampling of the original input points. We show that KT-KRR requires O(n 3/2 ) time during training and and O( √ n) time during inference, while achieving an near-minimax optimal rate of m log n nwhen the kernel has finite dimension (Thm. 2). We show how our KT-KRR guarantees can also be extended to the infinite-dimension setting (Thm. 3). In Sec. 5, we apply our proposed methods to both simulated and real-world data. In line with our theory, KT-NW and KT-KRR outperform standard thinning baselines in terms of accuracy while retaining favorable runtimes. 
problem setup: 
 We now formally describe the non-parametric regression problem. Let x 1 , . . . , x n be i.i.d. samples from the data distribution P over the domain X ⊂ R d and w 1 , . . . , w n be i.i.d. samples from N (0, 1). Then define the response variables y 1 , . . . , y n by the follow data generating process:EQUATIONwhere f ⋆ : X → Y ⊂ R is the regression function and v i ≜ σw i for some noise level σ > 0. Our task is to build an estimate for f ⋆ given the n observed points, denoted by S in ≜ ((x 1 , y 1 ), . . . , (x n , y n )).Nadaraya-Watson (NW) estimator. A classical approach to estimate the function f ⋆ is kernel smoothing, where one estimates the function value at a point z using a weighted average of the observed outcomes. The weight for outcome y i depends on how close x i is to the point z; let κ : R d → R denote this weighting function such that the weight for x i is proportional to κ(∥x i -z∥/h) for some bandwidth parameter h > 0. Let k : R d × R → R denote a shift-invariant kernel defined asEQUATIONThen this smoothing estimator, also known as Nadaraya-Watson (NW) estimator, can be expressed asEQUATIONwhenever the denominator in the above display is non-zero. In the case the denominator in ( 3) is zero, we can make a default choice, which for simplicity here we choose as zero. We refer to the estimator (3) as FULL-NW estimator hereafter. One can easily note that FULL-NW requires O(n) storage for the input points and O(n) kernel queries for inference at each point.Kernel ridge regression (KRR) estimator. Another popular approach to estimate f ⋆ is that of non-parametric (regularized) least squares. The solution in this approach, often called as the kernel ridge regression (KRR), is obtained by solving a least squares objective where the fitted function is posited to lie in the RKHS H of a reproducing kernel k, and a regularization term is added to the objective to avoid overfitting. 1 Overall, the KRR estimate is the solution to the following regularized least-squares objective, where λ > 0 denotes a regularization hyperparameter:EQUATION)Like NW, an advantage of KRR is the existence of a closed-form solutionEQUATIONEQUATIONNotably, the estimate f full,λ , which we refer to as the FULL-KRR estimator, can also be seen as yet another instance of weighted average of the observed outcomes. Notably, NW estimator imposes that the weights across the points sum to 1 (and are also non-negative whenever k is), KRR allows for generic weights that need not be positive (even when k is) and need not sum to 1. We note that naïvely solving f full,λ requires O(n 2 ) kernel evaluations to compute the kernel matrix, O(n 3 ) to compute a matrix inverse, and O(n) kernel queries for inference at each point. One of our primary goals in this work is to tackle this high computational cost of FULL-KRR.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: In Sec. 6, we discuss limitations as well as future work to address these limitations.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community.Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
speeding up non-parametric regression: 
 We begin with a general approach to speed up regression by thinning the input datasets. While computationally superior, a generic approach suffers from a loss of statistical accuracy motivating the need for a strategic thinning approach. To that end, we briefly review kernel thinning and finally introduced our supervised kernel thinning approach.Our generic approach comprises two main steps. First, we compress the input data by choosing a coreset S out ⊂ S in of size n out ≜ ∥S out ∥. Second, we apply our off-the-shelf non-parametric regression methods from Sec. 2 to the compressed data. By setting n out ≪ n, we can obtain notable speed-ups over the FULL versions of NW and KRR.Before we introduce the thinned versions of NW and KRR, let us define the following notation. Given an input sequence S in and output sequence S out , define the empirical probability measuresP in ≜ 1 n(x,y)∈Sin δ (x,y) and Q out ≜ 1 n out (x,y)∈Sout δ (x,y) .Thinned NW estimator. The thinned NW estimator is the analog of Full-NW except that S in is replaced by S out in (3) so that the thinned-NW estimator is given byEQUATIONwhenever the denominator in the display is not zero; and 0 otherwise. When compared to the FULL-NW estimator, we can easily deduce the computational advantage of this estimator: more efficient O(n out ) storage as well as the faster O(n out ) computation for inference at each point.Thinned KRR estimator. Similarly, we can define the thinned KRR estimator asEQUATIONα ′ ≜ (K ′ + n out λ ′ I nout ) -1    y ′ 1 . . . y ′ nout    ∈ R nout and K ′ ≜ [k(x ′ i , x ′ j )] nout i,j=1 ∈ R nout×noutgiven some regularization parameter λ ′ > 0. ). Our goal is to provide good computational benefits without trading off statistical error. Moreover, we may be able to do better by leveraging the underlying geometry of the input points and summarize of the input distribution more succinctly than i.i.d. sampling.A subroutine central to our approach is kernel thinning (KT) from Dwivedi and Mackey [10, Alg. 1]. We use a variant called KT-COMPRESS++ from Shetty et al. [22, Ex. 6 ] (see full details in App. A), which provides similar approximation quality as the original KT algorithm of Dwivedi and Mackey [10, Alg. 1], while reducing the runtime from O(n2 ) to O(n log 3 n). 2 Given an input kernel k ALG and input points S in , KT-COMPRESS++ outputs a coreset S KT ⊂ S in with size n out ≜ √ n ≪ n. In this work, we leverage two guarantees of KT-COMPRESS++. Informally, S KT satisfies (with high probability):EQUATIONEQUATIONwhere C 1 , C 2 > 0 are constants that depend on the properties of the input kernel k ALG and the chosen failure probability of KT-COMPRESS++, R in characterizes the radius of {x i } n i=1 , and N kALG (B 2 (R in ), 1/n out ) denotes the kernel covering number of H(k ALG ) over the ball B 2 (R in ) ⊂ R d at a specified tolerance (see Sec. 4.2 for formal definitions).At its highest level, KT provides good approximation of function averages. The bound (10) (formally stated in Lem. 1) controls the worst-case point-wise error, and is near-minimax optimal by Phillips and Tai [18, Thm. 3.1] . In the sequel, we leverage this type of result to derive generalization bounds for the kernel smoothing problem. The bound (11) (formally stated in Lem. 2) controls the integration error of functions in H(k ALG ) and is near-minimax optimal by Tolstikhin et al. [25, Thm. 1, 6] . In the sequel, we leverage this type of result to derive generalization bounds for the KRR problem.We show how the approximation results from kernel thinning can be extended to the regression setting. We construct two meta-kernels, the Nadaraya-Watson meta-kernel k NW and the ridge-regression meta-kernel k RR , which take in a base kernel k (defined over X only) and return a new kernel (defined over X × Y). When running KT, we set this new kernel as k ALG .A tempting choice of kernel for KT-NW is the kernel k itself. That is, we can thin the input points using the kernelEQUATION)This choice is sub-optimal since it ignores any information in the response variable y. For our supervised learning set-up, perhaps another intuitive choice would be to use KT withEQUATIONwhere ⊕ denotes vector concatenation (sox 1 ⊕ y 1 , x 2 ⊕ y 2 ∈ R d+1). While this helps improve performance, there remains a better option as we illustrate next.In fact, a simple but critical observation immediately reveals a superior choice of the kernel to be used in KT for NW estimator. We can directly observe that the NW estimator is a ratio of the averages of two functions:f numer (x, y)(•) ≜ k(x, •)⟨y, 1⟩ R and f denom (x, y)(•) ≜ k(x, •),over the empirical distribution P in (7) . Recall that KT provides a good approximation of sample means of functions in an RKHS, so it suffices to specify a "correct" choice of the RKHS (or equivalently the "correct" choice of the reproducing kernel). We can verify that f denom lies in the RKHS associated with kernel k(x 1 , x 2 ) and f numer lies in the RKHS associated with kernel k(x 1 , x 2 ) • y 1 y 2 . This motivates our definition for the Nadaraya-Watson kernel:EQUATIONsince then we do have f denom , f numer ∈ H(k NW ). Intuitively, thinning with k RR should simultaneously provide good approximation of averages of f denom and f numer over P in (see the formal argument in Sec. 4.1). When S out = KT-COMPRESS++(S in , k NW , δ), we call the resulting solution to (8) the kernel-thinned Nadaraya-Watson (KT-NW) estimator, denoted by f KT .As we show in Fig. 1 (a), this theoretically principled choice does provide practical benefits in MSE performance across sample sizes. Input sample size n2 3 2 1 2 1Mean Squared Errork(x1, x2) k((x1 y1), (x2 y2)) kNW((x1, y1), (x2, y2)) 2 10 2 14Input sample size n While with NW estimator, the closed-form expression was a ratio of averages, the KRR estimate (5) can not be expressed as an easy function of averages. However, notice that L Sin in ( 4) is an average of the function ℓ f : X × Y → R defined asℓ f (x, y) ≜ f 2 (x) -2f (x)y + y 2 for f ∈ H(k).Thus, there may be hope of deriving a KT-powered KRR estimator by thinning L Sin with the appropriate kernel. Assuming f ∈ H(k), we can verify that f 2 lies in the RKHS associated with kernel k 2 (x 1 , x 2 ) and that -2f (x)y lies in the RKHS associated with kernel k(x 1 , x 2 ) • y 1 y 2 . We now define the ridge regression kernel byEQUATIONand we can verify that f 2 (x) -2f (x)y lies in the RKHS H(k RR ). 3 When S out ≜ KT-COMPRESS++(S in , k RR , δ), we call the resulting solution to (9) the kernel-thinned KRR (KT-KRR) estimator with regularization parameter λ ′ > 0, denoted f KT,λ ′ . We note that the kernel k RR also appears in [12, Lem. 4] , except our subsequent analysis comes with generalization bounds for the KT-KRR estimator. Like for NW, in Fig. 1 (b) we do a comparison for KRR-MSE across many kernel choices and conclude that the choice ( 15) is indeed a superior choice compared to the base kernel k and the concatenated kernel (13) .• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
main results: 
 We derive generalization bounds of our two proposed estimators. In particular, we bound the mean squared error (MSE) defined by ∥f -f ⋆ ∥ 2 2 = E X (f (X) -f ⋆ (X)) 2 .Our first assumption is that of a well-behaved density on the covariate space. This assumption mainly simplifies our analysis of Nadaraya-Watson and kernel ridge regression, but can in principle be relaxed. Assumption 1 (Compact support). Suppose X ⊂ B 2 (R in ) ⊂ R d . Thus, the points x 1 , . . . , x n are drawn from a distribution with density p that satisfies 0 < p min ≤ p(x) ≤ p max for all x ∈ X .For the analysis of the NW estimator, we define function complexity in terms of Hölder smoothness following prior work [26] . Definition 1. For L > 0 and β ∈ (0, 2], a function f :X → R is (β, L)-Hölder if for all x 1 , x 2 ∈ X , |f (x 1 ) -f (x 2 )| ≤ L∥x 1 -x 2 ∥ β if β ∈ (0, 1] and |f (x 1 ) -f (x 2 ) -⟨∇f (x 1 ), x 2 -x 1 ⟩| ≤ L∥x 1 -x 2 ∥ β if β ∈ (1, 2],where f is assumed to be continuously differentiable for β ∈ (1, 2] but not for β ∈ (0, 1].Our next assumption is that on the kernel: We require k to be reproducing kernel to allow for valid analysis for the KT-NW estimator. While typically the NW estimator does not require a reproducing kernel several popular choices in practice, like Gaussian kernel, Laplace kernel, Matérn, Wendland, Sinc, and B-spline kernels, do satisfy this assumption. Assumption 2 (Shift-invariant kernel). k is a shift-invariant (2), reproducing kernel k(x 1 , x 2 ) = κ(∥x 1 -x 2 ∥/h) such that h > 0 and κ : R → R is bounded, L κ -Lipschitz, square-integrable, and decreasing with rate satisfyingEQUATION)Note that this assumption encomposses reproducing kernels with sub-Gaussian, sub-exponential, and poly tails. We now present our main result for the KT-NW estimator.Theorem 1 (KT-NW). Suppose Assum. 1 and 2 hold. Suppose eitherf ⋆ ∈ Σ(β, L f ) with β ∈ (0, 1] or f ⋆ ∈ Σ(β, L f ) for β ∈ (1, 2] and p ∈ Σ(β -1, L p ), L p > 0 and 2β > d.Then for any fixed δ ∈ (0, 1], the KT-NW estimator (8) EQUATIONwith probability at least 1 -δ, for some positive constant C that does not depend on n.See App. B for the proof. Remark 1. When κ from Assum. 2 has compact support, the condition 2β < d is not necessary.Tsybakov and Tsybakov [26] , Belkin et al. [3] show that FULL-NW achieves a rate of O(n -2β 2β+d ), which is minimax optimal for the (β, L)-Hölder function class. Compared to the ST-NW rate of n -β 2β+d , KT-NW achieves strictly better rates for all β > 0 and d > 0, while retaining ST-NW's fast query time of O( √ n). Note that our method KT-NW has a training time of O(n log 3 n), which is not much more than simply storing the input points.We present our main result for KT-KRR using finite-rank kernels. This class of RKHS includes linear functions, polynomial function classes . Theorem 2 (KT-KRR for finite-dimensional RKHS). Assume f ⋆ ∈ H(k), Assum. 1 is satisfied, and k has rank m ∈ N. Let f KT,λ ′ denote the KT-KRR estimator with regularization parameterλ ′ = O( m log nout n∧n 2 out). Then with probability at least 1 -2δ -2e- ∥f ⋆ ∥ 2 k c 1 (∥f ⋆ ∥ 2 k +σ 2 ), the following holds:EQUATIONfor some constant C that does not depend on n or n out .See App. C for the proof. Under the same assumptions, Wainwright [28, Ex. 13.19] showed that the Full-KRR estimator f full,λ achieves the minimax optimal rate of O(m/n) in O(n 3 ) runtime. When n out = √ n log c n, the KT-KRR error rates from Thm. 2 match this minimax rate in O(n 3/2 ) time, a (near) quadratic improvement over the Full-KRR. On the other hand, standard thinning-KRR with similar-sized output achieves a quadratically poor MSE of order m √ n . Our method and theory also extend to the setting of infinte-dimensional kernels. To formalize this, we first introduce the notion of kernel covering number.k : Z × Z → R with B k ≜ {f ∈ H : ∥f ∥ H ≤ 1}, a set A ⊂ Z and ϵ > 0, the covering number N k (A, ϵ) is the minimum cardinality of all sets C ⊂ B k satisfying B k ⊂ h∈C {g ∈ B k : sup x∈A |h(x) -g(x)| ≤ ϵ}.We consider two general classes of kernels. Assumption 3. For some C d > 0, all r > 0 and ϵ ∈ (0, 1), andB 2 (r) = x ∈ R d : ∥x∥ 2 ≤ r , a kernel k is LOGGROWTH(α, β) when log N k (B 2 (r), ϵ) ≤ C d log(e/ϵ) α (r + 1) β with α, β > 0 and POLYGROWTH(α, β) when log N k (B 2 (r), ϵ) ≤ C d (1/ϵ) α (r + 1) β with α < 2.We highlight that the definitions above cover several popular kernels: LOGGROWTH kernels include finite-rank kernels and analytic kernels, like Gaussian, inverse multiquadratic (IMQ), and sinc [9, Prop. 2], while POLYGROWTH kernels includes finitely-many continuously differentiable kernels, like Matérn and B-spline [9, Prop. 3] . For clarity, here we present our guarantee for LOGGROWTH kernels and defer the other case to App. E. Theorem 3 (KT-KRR guarantee for infinite-dimensional RKHS). Suppose Assum. 1 is satisfied and k is LOGGROWTH(α, β) (Assum. 3). Then f KT,λ ′ with λ ′ = O(1/n out ) satisfies the following bound with probability at least 1 -2δ -2eEQUATION)for some constant C that does not depend on n or n out .See App. E for the proof. When n out = √ n, ST-KRR achieves an excess risk rate of n -1/2 log α n for k satisfying LOGGROWTH(α, β). While KT-KRR does not achieve a strictly better excess risk rate bound over ST-KRR, we see that in practice, KT-KRR still obtains an empirical advantage. Obtaining a sharper error rate for the infinite-dimensional kernel setting is an exciting venue for future work. 
experimental results: 
 We now present experiments on simulated and real-world data. On real-world data, we compare our KT-KRR estimator with several state-of-the-art KRR methods, including Nyström subsampling-based methods and KRR pre-conditioning methods. All our experiments were run on a machine with 8 CPU cores and 100 GB RAM. Our code can be found at https://github.com/ag2435/npr. We begin with some simulation experiments. For simplicity, letEQUATIONand follow (1) to generate {y i } n i=1 (see Fig. 2 ). We let the input sample size n vary between 2 8 , 2 10 , 2 12 , 2 14 and always set the output coreset size to be n out = √ n. For NW, we use the Wendland(0) kernel defined byEQUATIONFor KRR, we use the Gaussian kernel defined byk(x 1 , x 2 ) ≜ exp(- ∥x1-x2∥ 2 2 2h 2) for h > 0.We select the bandwidth h and regularization parameter λ ′ (for KRR) using grid search. Specifically, we use a held-out validation set of size 10 4 and run each parameter configuration 100 times to estimate the validation MSE since KT-KRR and ST-KRR are random.Ablation study. In Fig. 1 , we compare thinning with our proposed meta-kernel k ALG = k NW to thinning with the baseline meta-kernels ( 12) and (13) . For our particular regression function (20) , thinning with (12) outperforms thinning with (13) . We hypothesize that the latter kernel is not robust to the scaling of the response variables. By inspecting (21) , we see that ∥(x 1 ⊕ y 1 ) -(x 2 ⊕ y 2 )∥ 2 is heavily determined by the y i values when they are large compared to the values of x i -as is the case on the right side of Fig. 2 (when X > 0). Since P is a uniform distribution, thinning with (12) evenly subsamples points along the input domain X , even though accuractely learning the left side of Fig. 2 (when X < 0) is not needed for effective prediction since it is primarily noise. Validating our theory from Thm. 1, the best performance is obtained when thinning with k NW (14) , which avoids evenly subsampling points along the input domain and correctly exploits the dependence between X and Y .In Fig. 1 , we perform a similar ablation for KRR. Again we observe that thinning with k ALG ((x 1 , y 1 ), (x 2 , y 2 )) = k(x 1 , x 2 ) outperforms thinning with k ALG ((x 1 , y 1 ), (x 2 , y 2 )) = k(x 1 ⊕ y 1 , x 2 ⊕ y 2 ), while thinning with k ALG = k RR achieves the best performance.Comparison with FULL, ST, RPCHOLESKY. In Fig. 3 (a), we compare the MSE of KT-NW to FULL-NW, ST-NW (a.k.a "Subsample"), and RPCHOLESKY-NW across four values of n. This last method uses the pivot points from RPCHOLESKY as the output coreset S out . At all n we evaluated, KT-NW achieves lower MSE than ST-NW and RPCHOLESKY-NW. FULL-NW achieves the lowest (b) Kernel ridge regression estimator with Gaussian kernel (22) . In Fig. 3 (b), we compare the MSE of KT-KRR to FULL-KRR, ST-KRR (a.k.a "Subsample"), and the RPCHOLESKY-KRR method from Chen et al. [6, Sec. 4.2.2] , which uses RPCHOLESKY to select landmark points for the restricted KRR problem. We observe that KT-KRR achieves lower MSE than ST-KRR, but higher MSE than RPCHOLESKY-KRR and FULL-KRR. In Fig. 3 (b), we also observe that KT-KRR is orders of magnitude faster than FULL-KRR across the board, with runtime comparable to ST-KRR and RPCHOLESKY-KRR in both training and testing. We hypothesize that RPCHOLESKY-while it provides a good low-rank approximation of the kernel matrix-is not designed to preserve averages.We now move on to our experiments on real-world data using two popular datasets: the California Housing regression dataset from Pace and Barry [16] (https://scikit-learn.org/1. 5/datasets/real_world.html#california-housing-dataset; BSD-3-Clause license) and the SUSY binary classification dataset from Baldi et al. [2] (https://archive.ics.uci.edu/ dataset/279/susy; CC-BY-4.0 license).California Housing dataset (d = 8, N = 2 × 10 4 ). Tab. 1(a) compares the test MSE, train times, and test times. We normalize the input features by subtracting the mean and dividing by the standard deviation and use a 80-20 train-test split. For all methods, we use the Gaussian kernel (22) with bandwidth h = 10. We use λ = λ ′ = 10 -3 for FULL-KRR, ST-KRR, and KT-KRR and λ = 10 (b) SUSY dataset.Table 1 : Accuracy and runtime comparison on real-world data. Each cell represents mean ± standard error across 100 trials.SUSY dataset (d = 18, N = 5 × 10 6 ). Tab. 1(b) compares our proposed method KT-KRR (with h = 10, λ ′ = 10 -1 ) to several large-scale kernel methods, namely RPCholesky preconditioning [7] , FALKON [21] , and Conjugate Gradient (all with h = 10, λ = 10 -3 ) in terms of test classification error and training times. For the baseline methods, we use the Matlab implementation provided by Díaz et al. [7] (https://github.com/ eepperly/Robust-randomized-preconditioning-for-kernel-ridge-regression). In our experiment, we use 4 × 10 6 points for training and the remaining 10 6 points for testing. As is common practice for classification tasks, we use the Laplace kernel defined by k(x 1 , x 2 ) ≜ exp(-∥x 1 -x 2 ∥ 2 /h). All parameters are chosen with cross-validation.We observe that KT-KRR achieves test MSE between ST-KRR and RPCHOLESKY preconditioning with training time almost half that of RPCHOLESKY preconditioning. Notably, our Cython implementation of KT-COMPRESS++ thinned the four million training samples in only 1.7 seconds on a single CPU core-with further speed-ups to be gained from parallelizing on a GPU in the future. 
conclusions: 
 In this work, we introduce a meta-algorithm for speeding up two estimators from non-parametric regression, namely the Nadaraya-Watson and Kernel Ridge Regression estimators. Our method inherits the favorable computational efficiency of the underlying Kernel Thinning algorithm and stands to benefit from further advancements in unsupervised learning compression methods.The KT guarantees provided in this work apply only when f ⋆ ∈ H(k) for some base kernel k. In practice, choosing a good kernel k is indeed a challenge common to all prior work. Our framework is friendly to recent developments in kernel selection to handle this problem: Dwivedi and Mackey [9, Cor. 1] provide integration-error guarantees for KT when f ⋆ / ∈ H(k). Moreover, there are recent results on finding the best kernel (e.g., for hypothesis testing [8, Sec. 4.2] ). Radhakrishnan et al. [19] introduce the Recursive Feature Machine, which uses a parameterized kernel kM (x 1 , x 2 ) ≜ exp(-(x 1 -x 2 ) ⊤ M (x 1 -x 2 )/(2h 2)), and propose an efficient method to learn the matrix parameter M via the average gradient outer product estimator. An exciting future direction would be to combine these parameterized (or "learned") kernels with our proposed KT methods for non-parametric regression. 
acknowledgements: 
 AG is supported with funding from the NewYork-Presbyterian Hospital. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We do not release models or data as part of this paper.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We include URL and licenses for baseline code and datasets used in Sec. 5.2.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: This section details the KT-COMPRESS++ algorithm of Shetty et al. [22, Ex. 6] . In a nutshell, KT-COMPRESS (Alg. 2) takes as input a point sequence of size n, a compression level g, a (reproducing) kernel functions k ALG , and a failure probability δ. KT-COMPRESS++ first runs the KT-COMPRESS(g) algorithm of Shetty et al. [22, Ex. 4] to produce an intermediate coreset of size 2 g √ n. Next, the KT algorithm is run on the intermediate coreset to produce a final output of size √ n.KT-COMPRESS proceeds by calling the recursive procedure COMPRESS, which uses KT with kernels k ALG as an intermediate halving algorithm. The KT algorithm itself consists of two subroutines:(1) KT-SPLIT (Alg. 3a), which splits a given input point sequence into two equal halves with small approximation error in the k ALG reproducing kernel Hilbert space and (2) KT-SWAP (Alg. 3b), which selects the best approximation amongst the KT-SPLIT coresets and a baseline coreset (that simply selects every other point in the sequence) and then iteratively refines the selected coreset by swapping out each element in turn for the non-coreset point that most improves MMD kALG error. As in Shetty et al. [22, Rem. 3] , we symmetrize the output of KT by returning either the KT coreset or its complement with equal probability.Following Shetty et al. [22, Ex. 6] , we always default to g = ⌈log 2 log n + 3.1⌉ so that KT-COMPRESS++ has an overall runtime of O(n log 3 n). For the sake of simplicity, we drop any dependence on g in the main paper. Algorithm 1: KT-COMPRESS++ -Identify coreset of size √ n Input: point sequence Sin of size n, compression level g, kernel kALG, failure probability δ SC ← KT-COMPRESS(g, Sin, δ) // coreset of size 2 g √ n SC++ ← KT(SC, kALG, g g+2 g (βn+1) δ) // coreset of size √ n return SC++EQUATION)Dwivedi and Mackey [10, Thm. 1, Rmk. 4] show thatP(E KT,δ ) ≥ 1 -δ.We restate [10, Thm. 4 ] in our notation: Lemma 1 (L ∞ guarantee for KT-COMPRESS++). Let Z ⊂ R d and consider a reproducing kernelk ALG : Z × Z → R. Assume n/n out ∈ 2 N . Let S KT ≜ KT-COMPRESS++(S in , g, k ALG , δ) andAlgorithm 3a: KT-SPLIT -Divide points into candidate coresets of size ⌊n/2⌋Input: kernel ksplit, point sequence Sin = (xi) n i=1 , failure probability δ S (1) , S (2) ← {} // Initialize empty coresets: S (1) , S (2) have size i after round i σ ← 0 // Initialize swapping parameter for i = 1, 2, . . . , ⌊n/2⌋ do // Consider two points at a time (1) .append(x); S (2) .append(x ′ ) end return (S (1) , S (2) ), candidate coresets of size ⌊n/2⌋(x, x ′ ) ← (x2i-1, x2i) // Compute swapping threshold ai ai, σ ←get_swap_params(σ, b, δ n ) with b 2 = ksplit(x, x)+ksplit(x ′ , x ′ )-2ksplit(x, x ′ ) // Assign one point to each coreset after probabilistic swapping θ ← 2i-2 j=1 (ksplit(xj, x) -ksplit(xj, x ′ )) -2 z∈S (1) (ksplit(z, x) -ksplit(z, x ′ )) (x, x ′ ) ← (x ′ , x) with probability min(1, 1 2 (1 -θ a i )+) Sfunction get_swap_params(σ, b, δ): ai ← max(bσ √ 2 log(2/δ), b 2 ) σ 2 ← σ 2 +b 2 (1+(b 2 -2ai)σ 2 /a 2 i )+ return (ai, σ)Input: kernel kALG, point sequence Sin = (xi) n i=1 , candidate coresets (S (1) , S (2) ) S (0) ← baseline_coreset(Sin, size = ⌊n/2⌋) // Compare to baseline (e.g., standard thinning)SKT ← S (ℓ ⋆ ) for ℓ ⋆ ← argmin ℓ∈{0,1,2} MMD k ALG (Sin, S (ℓ) ) // Select best coreset// Swap out each point in SKT for best alternative in Sin while ensuring no point is repeated in SKT for i = 1, . . . , ⌊n/2⌋ doSKT[i] ← argmin z∈{S KT [i]}∪(S in \S KT ) MMD k ALG (Sin, SKT with SKT[i] = z) end return SKT, refined coreset of size ⌊n/2⌋ define P in ≜ 1 n z∈Sin δ z and Q out ≜ 1 nout z∈SKT δ z .Then on event E KT,δ , the following bound holds:EQUATIONEQUATION)log 1 δ + d log 1 + L k ALG ∥kALG∥ ∞ (R kALG,n + R) , L kALG ≜ sup z1,z2,z3∈Z |kALG(z1,z2)-kALG(z1,z3)| ∥z2-z3∥2, andEQUATIONfor some universal positive constant c.Proof. ∥h∥ k ALG ≤1 |(P in -Q out )h| ≤ inf ϵ∈(0,1) Sin⊂A 2ϵ + 2∥kALG∥ 1/2 ∞,in nout W kALG (n, n out , δ, A, ϵ)whereEQUATION)for some universal postive constant c.Proof. The claim follows from replacing k in [9, Thm. 2] with k ALG and replacing the sub-Gaussian constant of KT with that of KT-COMPRESS++ in [22, Ex. 5 ].B Proof of Thm. 1: KT-NWOur primary goal is to bound E Sin [( f KT (x 0 ) -f ⋆ (x 0 )) 2 ] for a fixed x 0 ∈ X . Once we have this bound, bounding ∥ f KT -f ⋆ ∥ 2 2is as straightforward as integrating over x 0 ∈ X . Consider the following decomposition:EQUATION)Define the random variablesη i ≜ 1 ∥Xi-x0∥ h ≤ 1 for i = 1, 2, . . . , n.Also define the eventEQUATIONSince X i are i.i.d. samples from P, it follows that η i are i.i.d. Bernoulli random variables with parameterEQUATIONwhere c 0 > 0 depends only on d and κ (see Assum. 2). Denote the denominator terms in f and f KT byEQUATIONrespectively, and the numerator terms in f and f KT byEQUATIONWe now consider two cases depending on the event E.Case I: Suppose event E c is satisfied. It follows from (33) that p(x 0 ) = 0, in which case f (x 0 ) = 0. Since S out ⊂ S in , it necessarily follows that p KT (x 0 ) = 0 and f KT (x 0 ) = 0. Thus, we can bound ( 29) and ( 30) byE Sin f KT (x 0 ) -f (x 0 ) 2 I[E c ] = 0 and E Sin f (x 0 ) -f ⋆ (x 0 ) 2 I[E c ] = E Sin (0 -f ⋆ (x 0 )) 2 I[E c ] ≤ (f ⋆ ) 2 (x 0 )P(E c ) ≤ (f ⋆ ) 2 (x 0 )(1 -p) n ≤ (f ⋆ ) 2 (x 0 ) exp -Cnh dfor some positive constant C that does not depend on n. Note that these are low-order terms compared to the rest of the calculations, so we may ignore them in the final bound.Case II: Otherwise, we may assume event E is satisfied. Let us first bound (29) . On event E KT,δ (23) , we claim thatEQUATIONWe defer the proof to App. B.1.Letting X ≜ (X 1 , . . . , X n ) and Y ≜ (Y 1 , . . . , Y n ) denote the x and y components of S in , respectively, we can further decompose ( 30) byE Sin f (x 0 ) -f ⋆ (x 0 ) 2 I[E] = E X E Y |X f (x 0 ) -E Y |X f (x 0 ) 2 I[E] + E X E Y |X f (x 0 ) -f ⋆ (x 0 ) 2 I[E] ,where the first RHS term corresponds to the variance and the second RHS term corresponds to the bias. We claim thatEQUATIONEQUATIONfor some constant C > 0 that does not depend on either n or h. We defer the proofs to App. B.2 and B.3.Combining (35) to (37), we haveEQUATION)Note that h d ≤ 1, so the Cd log 2 n nh 2dterm dominates the C nh d term. Thus, the optimal choice of bandwidth h comes from balancingEQUATIONFinally, we must verify our growth rate assumption on p in (35) is satisfied. Since β > 0, we havep (32) ≥ c 0 p min h d (39) = c ′ 0 n -d 2β+2d =⇒ lim n→∞ p √ d n log n = ∞.Plugging (39) into (38) yields the advertised bound (17) .We first provide a generic result for approximating the numerator and denominator terms defined in (33) and (34). Lemma 3 (Simultaneous L ∞ bound using KT-COMPRESS++ with k NW ). Suppose k satisfies Assum. 2. Given S in , the following bounds hold on the event E KT,δ :EQUATIONwhere c a , c p > 0 are constants that do not depend on d or n.See App. B.1.1 for the proof. In the sequel, we will simply treat the log(1/δ) term as a constant, meaning the log n terms dominate in the expressions.With this lemma in hand, let us prove the claim (35). Define the following events:A ≜ { p KT (x 0 ) = 0} B ≜ { p KT (x 0 ) ̸ = 0} C ≜ p(x 0 ) ≥ p 2 .On event E, consider the following decomposition:EQUATIONEQUATION)Bounding (42). Note that almost surely, we have| f (x 0 )| ≤ Y max and | f KT (x 0 )| ≤ Y max .Thus, we haveE Sin f KT (x 0 ) -f (x 0 ) 2 I[C c ] ≤ 4Y 2 max P n p(x 0 ) < np 2 (i) = 4Y 2 max P n i=1 η i -np < np 2 -np (ii) ≤ c 0 exp -c 1 nh d ,where (i) follows from subtracting np from both sides of the probability statement and (ii) follows from concentration of Bernoulli random variables (see App. B.2).Bounding (43). Note that on event E KT,δ ∩ C, we haveEQUATION)where step (i) follows from applying (40) and substituting p ≥ p 2 . Hence the events E KT,δ and A ∩ C are mutually exclusive with probability 1, thereby yieldingE Sin f KT (x 0 ) -f (x 0 ) 2 I[E KT,δ ∩ A ∩ C] = 0.Bounding (44). On the event B ∩ C, we have f KT (x 0 ) = AKT(x0) pKT(x0) and f (x 0 ) = A(x0) p(x0) , which yieldsf KT (x 0 ) -f (x 0 ) = AKT pKT -A(x0) p(x0) = AKT(x0)• p(x0)-A(x0)• pKT(x0) p(x0)• pKT(x0) = ( AKT(x0)-A(x0))• p(x0)+ A(x0)•( p(x0)-pKT(x0)) p(x0)• pKT(x0) ≤ | AKT(x0)-A(x0)|• p(x0)+ A(x0)•| p(x0)-pKT(x0)| p(x0)• pKT(x0)We can invoke (40) and (41) to bound | p(x 0 ) -p KT (x 0 )| and | A KT (x 0 ) -A(x 0 )| respectively. Thus, we haveE Sin f KT (x 0 ) -f (x 0 ) 2 I[E KT,δ ∩ B ∩ C] ≤ ca √ d n log n• p(x0)+cp √ d n log n• A(x0) p(x0)• pKT(x0) 2 ≤ 2d•log 2 n n ca pKT(x0) 2 + A(x0) p(x0) 2 cp pKT(x0) 2 (i) ≤ 2d•log 2 n n c 2 a +Y 2 max c 2 p pKT(x0) 2 (ii) ≤ Cd log 2 n nh 2d ,for some positive constant C that does not depend on n, where step (i) uses the fact that A(x0) p(x0) ≤ Y max and step (ii) uses the lower bound on p KT (x 0 ) from (45). Combining (42) to (44), we haveE Sin f KT (x 0 ) -f (x 0 ) 2 I[E KT,δ ] ≤ c 0 exp -c 1 nh d + Cd log n nh 2d .Note that the second term dominates so that we may drop the first term with slight change to the value of the constant C in the bound (35).B.1.1 Proof of Lem. 3: Simultaneous L ∞ bound using KT-COMPRESS++ with k NW We first decompose k NW as k NW ((x 1 , y 1 ), (x 2 , y 2 )) = k 1 ((x 1 , y 1 ), (x 2 , y 2 )) + k 2 ((x 1 , y 1 ), (x 2 , y 2 )),whereEQUATIONEQUATION)and note thatEQUATIONThis fact will be useful later for proving simultaneous L ∞ approximation guarantees for A and p.Given that k satisfies Assum. 2, we want to show that k NW defined by ( 14) satisfies the Lipschitz and tail decay properties, so that we may apply Lem. 1. Note thatEQUATIONWe claim that kernel k NW satisfiesEQUATIONBy [10, Rmk. 8], we haveEQUATIONwhere κ † is defined by (16) . Applying (49) and (50), we haveEQUATIONwhere the last inequality follows from Assum. 2 for some universal positive constant c ′ .Since Assum. 1 is satisfied, R is constant. Applying (53) to M kNW (n, n out , d, δ, R) as defined by (25) , we have the boundM kNW (n, n out , d, δ, R) ≤ c ′′ log nout δ log 8 δ + 5 √ d log nfor some positive constant c ′′ . Substituting this into (24), we have∥(P in -Q out )k NW ∥ ∞ ≤ c 1 ∥k∥ ∞ (1+Y 2 max ) nout log nout δ log 8 δ + 5 √ d log n ≤ c 2 ∥k∥ ∞ (1+Y 2 max ) nout √ d( √ log n + log(1/δ)) 2 ,for some positive constants c 1 , c 2 . By definition,∥(P in -Q out )k NW ∥ ∞ = sup z∈Z ⟨(P in -Q out )k NW , k NW (•, z)⟩ kNW .Define k 1 and k 2 by ( 46) and (47), respectively, and note that k1 (•, z), k 2 (•, z) ∈ H(k NW ) for all z ∈ Z.We want to show thatsup z∈Z ⟨(P in -Q out )k NW , k 1 (•, z)⟩ kNW ≤ c 2 ∥k∥ ∞ (1+Y 2 max ) nout √ d( √ log n + log(1/δ)) 2 andsup z∈Z ⟨(P in -Q out )k NW , k 2 (•, z)⟩ kNW ≤ c 2 ∥k∥ ∞ (1+Y 2 max ) nout √ d( √ log n + log(1/δ)) 2 ,which would imply (40) and (41) (after simplifying all terms besides n, d, and δ).The first inequality follows from replacing all occurrences of the test function k NW (•, (x, y)) in the proof of Lem. 1 with the function k 1 (•, x) and noting that ⟨k NW (•, (x i , y i )), k 1 (•, (x, y))⟩ kNW = ⟨k 1 (•, x i ), k 1 (•, x)⟩ k1 from the fact that H(k NW ) = H(k 1 ) ⊕ H(k 2 ) (48).The second inequality follows from replacing all occurrences of the test function k NW (•, (x, y)) in the proof of Lem. 1 with the function k 2 (•, x) and noting that ⟨k NW (•, (x i , y i )), k 2 (•, (x, y))⟩ kNW = ⟨k 2 (•, (x i , y i )), k 2 (•, (x, y))⟩ k2 , again from the fact that H(k NW ) = H(k 1 ) ⊕ H(k 2 ) (48).Proof of claim (50). We leverage the fact that the Lipschitz constants defined by (26) satisfies the following additivity property. Letting Z = S in , we haveL kNW = sup z1,z2,z3∈Z |kNW(z1,z2)-kNW(z1,z3)| ∥z2-z3∥2 ≤ sup z1,z2,z3∈Z |k1(z1,z2)-k1(z1,z3)| ∥z2-z3∥2 + sup z1,z2,z3∈Z |k2(z1,z2)-k2(z1,z3)| ∥z2-z3∥2 = L k1 + L k2 .We proceed to bound L k1 and L k2 separately. Note thatL k1 = L k .Applying the definition (26) to L k2 , we haveL k2 = sup z1=(x1,y1) z2=(x2,y2) z3=(x3,y3) |k(x1,x2)y1y2-k(x1,x3)y1y3| √ ∥x2-x3∥ 2 +∥y2-y3∥ 2 = sup z1=(x1,y1) z2=(x2,y2) z3=(x3,y3) |y1|•|k(x1,x2)y2-k(x1,x2)y3+k(x1,x2)y3-k(x1,x3)y3| √ ∥x2-x3∥ 2 +∥y2-y3∥ 2 = sup z1=(x1,y1) z2=(x2,y2) z3=(x3,y3) |y1|•|k(x1,x2)(y2-y3)+(k(x1,x2)-k(x1,x3))y3| √ ∥x2-x3∥ 2 +∥y2-y3∥ 2 ≤ sup z1=(x1,y1) z2=(x2,y2) z3=(x3,y3) |y1|•|k(x1,x2)(y2-y3)| √ ∥x2-x3∥ 2 +∥y2-y3∥ 2 + sup z1=(x1,y1) z2=(x2,y2) z3=(x3,y3) |y1|•|(k(x1,x2)-k(x1,x3))y3| √ ∥x2-x3∥ 2 +∥y2-y3∥ 2 ≤ Y max ∥k∥ ∞ + L k Y 2 max .Putting together the pieces yields the claimed bound.Proof of claim (51). We aim to show that R kNW,n is not much larger than R k,n . Note that k NW can be rewritten ask NW ((x 1 , y 1 ), (x 2 , y 2 )) = (1 + y 1 y 2 )k(x 1 , x 2 ).We define the setsEQUATIONEQUATIONnoting thatR k,n = inf Γ and R kNW,n = inf Γ ⋆by definition (27) .Suppose r ∈ Γ. Then for any (x 1 , y 1 ), (x 2 , y 2 ) such that∥x 1 -x 2 ∥ 2 + ∥y 1 -y 2 ∥ 2 ≥ r 2 + 4Y 2 max , it must follow that ∥x 1 -x 2 ∥ 2 ≥ r 2 (since ∥y 1 -y 2 ∥ 2 ≤ 4Y 2max by triangle inequality). Since r satisfies (54), it must follow that|k(x 1 , x 2 ) • (y 1 y 2 + 1)| ≤ |k(x 1 , x 2 )|(Y 2 max + 1) ≤ ∥k∥ ∞ n (Y 2 max + 1) (49) = ∥kNW∥ ∞ n , meaning r 2 + 4Y 2 max ∈ Γ ⋆, where recall Γ ⋆ is defined by (55). Thus, we haveR kNW,n ≤ R k,n + 4Y 2 max ≤ R k,n + 2Y max as desired.Suppose event E (31) is satisfied. Define the shorthand for the variance:EQUATIONwhere we have used the fact thatE[Y | X = •] = f ⋆ (•)by assumption (1) .Note that on event E, we haveσ 2 (x 0 ; X) (56) = E Y |X n i=1 Yik(Xi,x0) n i=1 k(Xi,x0) - n i=1 f ⋆ (Xi)k(Xi,x0) n i=1 k(Xi,x0) 2 = E Y |X n i=1 vik(Xi,x0) n i=1 k(Xi,x0) 2 = Var[v 1 ] • n i=1 k(Xi,x0) 2 ( n i=1 k(Xi,x0)) 2 ,where recall v 1 , . . . , v n are i.i.d. random variables with Var[v i ] = σ 2 by (1). Taking the expectation w.r.t. X 1 , . . . , X n and leveraging symmetry, we haveEQUATION)σ 2 X can be bounded byσ 2 X ≤ E X k 2 (X1,x0) ( n i=1 k(Xi,x0)) 2 I n i=1 η i ≤ np 2 + 2 np 2 E X k 2 (X 1 , x 0 ) (i) ≤ P n i=1 η i ≤ np 2 + 2 np 2 X k 2 (x 1 , x 0 )p(dx 1 ),where step (i) follows from the fact thatk 2 (X1,x0) ( n i=1 k(Xi,x0))2 ≤ 1. Using Bernstein's inequality [28, Prop. 2.14], the first term can be bounded byP n i=1 η i ≤ np 2 = P n i=1 η i -np ≤ -np 2 ≤ exp - (np) 2 2(np(1-p)+np/3) ≤ exp -c 1 nh dfor some universal positive constant c. Applying the fact that p is bounded by Assum. 1 and κ is square-integrable by Assum. 2, we can bound the second term by2 np 2 X k 2 (x 1 , x 0 )p(dx 1 ) (i) ≤ 2 np 2 p max h d R d κ 2 (u)du (ii) ≤ 4 (np) 2 c 1 h d ≤ c2 n 2 h d ,for some positive constant c 2 that does not depend on either h or n. Substituting these expressions into (57) yields a bound on E X σ 2 (x 0 ; X) I[E] as desired.Define the following shorthand for the bias:b(x 0 ; X) ≜ E Y |X f (x 0 ) -f ⋆ (x 0 ).We state a more detailed version of the claim.Lemma 4 (Bias of Nadaraya-Watson). Suppose Assum. 1 and 2 are satisfied and the event E (31) holds.If f ⋆ ∈ Σ(β, L f ) for β ∈ (0, 1], L f > 0, then the following statements hold true for any x 0 ∈ X :(I.a) If k is compactly supported, then b 2 (x 0 ; X) ≤ L 2 f h 2β . (I.b) If k has tail decay satisfying (16), then b 2 (x 0 ; X) ≤ c • L 2 f (h 2β ∨ h 2β-d ) for some positive constant c. Now suppose f ⋆ ∈ Σ(β, L f ) for β ∈ (1, 2], L f > 0 and the density p of the marginal distribution of X satisfies p ∈ Σ(β -1, L p ), L p > 0. Let σ X be defined by (57). Then the following statements hold for any x 0 ∈ X :(II.a) If k is compactly supported, then b 2 (x 0 ; X) ≤ L f + ∥∇f (x 0 )∥L p p -1 min h 2β + σ 2 X . (II.b) If k has tail decay satisfying (16), then b 2 (x 0 ; X) ≤ c • L f + ∥∇f (x 0 )∥L p p -1 min (h 2β ∨ h 2β-d ) + σ 2X for some positive constant c. EQUATIONwhere step (i) follows from our assumption that f ⋆ ∈ Σ(β, L f ) and step (ii) follows from our assumption that k is compactly supported, so k(x, x 0 ) = 0 whenever ∥x -x 0 ∥ > h.Proof of claim (I.b). Consider the following decomposition of the bias:b(x 0 ; X) ≤ n i=1 k(Xi,x0)I[∥Xi-x0∥≤h](f ⋆ (Xi)-f ⋆ (x0)) n k=1 k(X k ,x0)I[∥Xi-x0∥≤h] + n i=1 k(Xi,x0)I[∥Xi-x0∥>h](f ⋆ (Xi)-f ⋆ (x0)) n k=1 k(X k ,x0).Note that the first RHS term can be bounded using (58). To bound the second RHS term, we introduce the following high-probability eventE ≜ {|I| = O(nh d/2 )}, where I ≜ {i ∈ [n] : ∥X i -x 0 ∥ = O(h 1/2 )}.On this event, we may apply the Lipschitz property of f ⋆ again to obtainn i=1 k(Xi,x0)I[∥Xi-x0∥>h](f (Xi)-f (x0)) n k=1 k(X k ,x0) ≤ L f n i=1 k(Xi,x0)∥Xi-x0∥ β k∈I k(X k ,x0) (i) ≤ L f n i=1 κ( ∥X i -x 0 ∥ h )∥Xi-x0∥ β k∈I k(X k ,x0) ≤ L f n max u≥0 κ(u)(uh) β k∈I k(X k ,x0), where step (i) follows from the fact that k is a shift-invariant kernel by Assum. 2. Note that when κ has tail decay satisfying (16) , we havemax x≥0 κ(u)(uh) β = c • h βfor some positive constant c. This implies on the event E, the following bound holds b(x 0 ; X) ≤ L f h β + L f • c • nh β nh d/2 = c • L f (h β ∨ h β-d/2 ). B.3.2 Bias when β ∈ (1, 2]Proof of claim (II.a). For notational simplicity, let x 0 = 0. We further assume the density p of X satisfies (i) p ∈ [p ℓ , p u ] for some 0 < p ℓ < p u < ∞ and p ∈ Σ(β -1, L p ). Define B ≜ ∥∇f (0)∥ 2 and denote θ ≜ (L p , L f , p ℓ , p u , B) to be the collection of parameters. Constants c > 0 in this proof may differ from line to line, but depends only up to parameters in θ.By definition,b 2 (0) ≜ E n i,j=1 G i G j I[E 0 ]whereG i = (f (Xi)-f (0))k(Xi) n k=1 k(X k ) . Further define A(X i , X j ) ≜ k(Xi)k(Xj ) ( n k=1 k(X k )) 2 I[E 0 ] ≥ 0 for i ̸ = j,and denote E i,j to be the conditional expectation of law X i , X j , while fixing all other randomnessX k , k ̸ = i, k ̸ = j. Under independence of X i , i ∈ [n],it is safe to say that E i,j is the expectation of marginal law of X i , X j , while fixing other randomness as constants.DefineR(x) = f (x) -f (0) -⟨∇f (0), x⟩so that by simple add-substract algebra, we haveE i,j [G i G j I[E 0 ]] = ⟨∇f (0), x i ⟩⟨∇f (0), x j ⟩A(x i , x j )p(x i )p(x j )dx i dx j + 2 ⟨∇f (0), x i ⟩R(x j )A(x i , x j )p(x i )p(x j )dx i dx j + R(x i )R(x j )A(x i , x j )p(x i )p(x j )dx i dx j .Observe that A(x i , x j ) are even functions for both arguments, and x → ⟨∇f (0), x⟩ is an odd function.So we can see ⟨∇f (0), x i ⟩A(x i , x j )p(0)dx i = 0.E i,j [G i G j I[E 0 ]] = ⟨∇f (0), x i ⟩⟨∇f (0), x j ⟩A(x i , x j )(p(x i ) -p(0))(p(x j ) -p(0))dx i dx j + 2 ⟨∇f (0), x i ⟩R(x j )A(x i , x j )(p(x i ) -p(0))p(x j )dx i dx j + R(x i )R(x j )A(x i , x j )p(x i )p(x j )dx i dx j .From the assumptions, we see |R(x i )| ≤ L f ∥x i ∥ β so that we haveE i,j [G i G j I[E 0 ]] ≤ c ∥x i ∥ β ∥x j ∥ β A(x i , x j )dx i dx j + c ∥x i ∥ β ∥x j ∥ β A(x i , x j )p(x i )dx i dx j + c ∥x i ∥ β ∥x j ∥ β A(x i , x j )p(x i )p(x j )dx i dx j ≤ c ∥x i ∥ β ∥x j ∥ β A(x i , x j )dx i dx j .So overall we haveE[G i G j I[E 0 ]] ≤ cE ∥Zi∥ β ∥Zj ∥ β k(Zi)k(Zj ) ( n k=1 k(Z k )) 2 I[E 0 ] ,where Z i 's are independent and uniform measures on the domain -such integration is established by pulling out all the p u , p ℓ for all p's of X i 's. Sob 2 (0) ≤ cE n i=1 |Zi| β k(Zi) n k=1 k(Z k ) 2 I[E 0 ]. Now let's induce the same reasoning as done in β ∈ (0, 1] case, so thatb 2 (0) ≤ O(h 2β ) ∨ O(h 2β-d ).Proof of claim (II.b). The proof follows by similar logic as (II.a) combined with the truncation argument of (I.b).We rely on the localized Gaussian/Rademacher analysis of KRR from prior work [28] . Define the Gaussian critical radius ε n > 0 to be the smallest positive solution to the inequalityEQUATIONB H (3) is the ∥•∥ k -ball of radius 3 and w i i.i.d.∼ N (0, 1).Assumption 4. Assume that ∥f ⋆ ∥ k ∈ B k (R) and f KT,λ ′ ∈ B k (c † R), for some constant c † > 0 .Note that for any g ∈ B k (c † R), we have∥g∥ ∞ ≤ sup x∈X ⟨g, k(•, x)⟩ k (i) ≤ sup x∈X ∥g∥ k ∥k(•, x)∥ k ≤ ∥g∥ k ∥k∥ ∞ ≤ c † R ∥k∥ ∞ ≜ B,where step (i) follows from Cauchy-Schwarz. Thus, the function class B k (c † R) is B-uniformly bounded. Now define the Rademacher critical radius δ n > 0 to be the smallest positive solution to the inequalityEQUATIONand ν i = ±1 each with probability 1/2.Finally, we use the following shorthand to control the KT approximation error term, EQUATIONLet f KT,λ ′ denote the KT-KRR estimator with regularization parameterλ ′ ≥ 2ξ 2 n where ξ n ≜ ε n ∨ δ n ∨ 4 √ C m (∥f ⋆ ∥ k + 1)η n,k .Then with probability at least 1 -2δ -2e- nδ 2 n c 1 (b 2 +σ 2 ) , we have ∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ c ξ 2 n + λ ′ ∥f ⋆ ∥ 2 k + cδ 2 n .(64) where recall δ is the success probability of KT-COMPRESS++ (23) .See App. D for the proof. We set λ ′ = 2ξ 2 n , so that (64) becomes∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ 3cξ 2 n ∥f ⋆ ∥ 2 k + cδ 2 n .(65) It remains to bound the quantities ε n (59), δ n (60), and η n,k (61). We claim thatε n ≤ c 0 σ R m n(66)δ n ≤ c 1 b m n (67) η n,k ≤ c 2 √ m•log nout•log(1/δ) nout.(68) for some universal positive constants c 0 , c 1 , c 2 . Now setR = ∥f ⋆ ∥ k δ = e -1/R 4 . Thus, we have ξ n ≤ c ′ ( σ ∥f ⋆ ∥ k ∨ b ∨ 4 √ Cm ∥f ⋆ ∥ k ) √ m √ n∧nout .for some universal positive constant c ′ . Substituting this into (65) leads to the advertised bound (18) .Proof of claim (66). For finite rank kernels, μj = 0 for j > m. Thus, we have2 n n j=1 min{ε 2 , μj } = 2 n √ mε 2 .From the critical radius condition (59), we want2 n √ mε 2 ≤ R 4σ ε 2 , so we may set ε n ≃ σ R m n .Proof of claim (67). By similar logic as above, we have 2n n j=1 min{δ 2 , µ j } = 2 n √ mδ 2 .From the critical radius condition (60), we want2 n √ mδ 2 ≤ 1 b δ 2 , so we may set δ n ≃ b √ 2 m n .Proof of claim (68). Consider the linear operator T : H → R m that maps a function to the coefficients in the vector space spanned by{ϕ i } m i=1 . Note that ∥T ∥ = ∥T f ∥ ∞ ∥f ∥ k ≤ ∥k∥ ∞ Since the image of T has dimension m, we have rank(T ) ≤ m. Moreover, ∥k∥ ∞ ≤ µ 1 • R 2in . Now we can invoke [24, Eq. 14] with ϵ = a/n out to obtainN k (B d 2 (R in ), a/n out ) ≤ N (T, a/n out ) ≤ (1 + µ 1 R 2 in n out /a) m .Taking the log on both sides and substituting this bound into (61), we haveη n,k = a 2 nout (2 + W k (n, n out , δ, R in , a nout )) ≤ a 2 nout (2 + log nout log(n/nout) δ • log 1 δ + log N k (B d 2 (R in ), a nout ) ) ≤ a 2 nout (2 + log nout log(n/nout) δ • log 1 δ + m log 1 + 2∥T ∥nout a ) ≤ c √ m•log nout•log(1/δ) noutfor some positive constant c that doesn't depend on m, n out , δ.We rescale our observation model (1) by ∥f ⋆ ∥ k , so that the noise variance is (σ/∥f ⋆ ∥ k ) 2 and our new regression function satisfies ∥f ⋆ ∥ k = 1. Our final prediction error should then be multiplied by ∥f ⋆ ∥ 2 k to recover a result for the original problem. For simplicity, denote σ = σ/∥f ⋆ ∥ k . For notational convenience, define an eventE = ∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ c(ξ 2 n + λ ′ ), and our goal is to show that E occurs with high-probability in terms of P, the probability regarding all the randomness. For that end, we introduce several events that are used throughout,EQUATIONwhere δ n is defined in (60) and H is the RKHS generated by k hence star-shaped. Further, we introduce two technical events A KT (u), B KT defined in (82) and ( 95) respectively, which are proven to occur with small probability, and define a shorthandE good ≜ A c KT (ξ n ) ∩ B c KT ∩ E conc ∩ E KT,δ .Equipped with these shorthands, observe the following inequality,EQUATIONwhere the second inequality is because E c lower ⊆ E c lower ∩ E due to the assumption λ ′ ≥ 2ξ 2 n ≥ 2δ 2 n . If we are able to show the set inclusion {E good ∩ E lower } ⊆ {E ∩ E lower } and that P(E c good ) is small, we are able refine (70) to the followingP(E) ≥ P(E good ∩ E lower ) + P(E c lower ) ≥ 1 -P(E c good ) -P(E c lower ) + P(E lower ) = 1 -P(E c good ), where the last quantity 1 -P(E c good ) would be large. To complete this proof strategy, we claim the set inclusion {E good ∩ E lower } ⊆ {E ∩ E lower } (71) to hold and prove it in App. D.1 and further claimEQUATIONwhich verify in App. D.2.Putting the pieces together, claims (71) and (72) collectively impliesP(E) ≥ 1 -c ′′ δ + e -c ′ nδ 2 n /(B 2 H ∧ σ 2 )as desired.There are several intermediary steps we take to show the set inclusion of interest (71). We introduce the shorthand∆ KT ≜ f KT,λ ′ -f ⋆ .By invoking Propositions and basic inequalities to come, we successively show the following chain of set inclusionsE good ∩ E lower ⊆ E good ∩ E lower ∩ {∥ ∆ KT ∥ 2 nout ≤ c(ξ 2 n + λ ′ )} (73) ⊆ E good ∩ E lower ∩ {∥ ∆ KT ∥ 2 n ≤ c(ξ 2 n + λ ′ )} (74) ⊆ E good ∩ E lower ∩ E (75) ⊆ E lower ∩ E.(76) Note that step (76) is achieved trivially by dropping E good . Further note that (74) is the crucial intermediary step after which we may apply uniform concentration across n independent samples. Proof of (74) leverages on the Proposition to come (Prop. 1) that allows ∥ • ∥ n and ∥ • ∥ nout to be exchangeable for finite rank kernels.Recovering step (73) Since f KT,λ ′ and f ⋆ are optimal and feasible, respectively for the central optimization problem of interestEQUATIONWith some algebra , may refine (77) toEQUATION)where ∆ KT = f KT,λ ′ -f ⋆ .Suppose that ∥ ∆ KT ∥ nout < ξ n , then we trivially recover (73) by adding λ ′ > 0. Thus, we assume that ∥ ∆ KT ∥ nout ≥ ξ n .Under the assumption ∥ ∆ KT ∥ nout ≥ ξ n , which is without loss of generality, we utilize the basic inequality (78) and control its stochastic component1 nout nout i=1 v ′ i ∆ KT (x ′ i ) ,with a careful case work to follow, which is technical by nature.Case where ∥ f KT,λ ′ ∥ k ≤ 2: Under such case, we introduce a technical eventA KT (u) ≜ ∃g ∈ F \ B 2 (δ n ) ∩ {∥g∥ nout ≥ u} such that 1 nout nout i=1 v ′ i g(x ′ i ) ≥ 3∥g∥ nout u , for any star-shaped function class F ⊂ H. Since ∥f ⋆ ∥ k = 1, triangle inequality implies ∥ ∆ KT ∥ k ≤ ∥ f KT,λ ′ ∥ k + ∥f ⋆ ∥ k ≤ 3.Moreover, on the event E lower (⊆ E good ), we have ∥ ∆ KT ∥ 2 > δ n . Thus, we may apply ∆ KT to the event A c KT (ξ n ) with F = B H (3) (i.e., the H-ball of radius 3) to attainEQUATIONUpper bounding the stochastic component of the basic inequality (78) by ( 79) and dropping the -∥ f KT,λ ′ ∥ 2 k term in (78), we have1 2 ∥ ∆ KT ∥ 2 nout ≤ c 0 ξ n ∥ ∆ KT ∥ nout + λ ′ .As a last step under the case ∥ f KT,λ ′ ∥ k ≤ 2, apply the quadratic formula (specifically, if a, b ≥ 0 andx 2 -ax -b ≤ 0, then x ≤ a 2 + b) to obtain ∥ ∆ KT ∥ 2 nout ≤ 4c 2 0 ξ 2 n + 2λ ′ .Case where ∥ f KT,λ ′ ∥ k > 2: Under such case, by assumption we have∥ f KT,λ ′ ∥ k > 2 > 1 ≥ ∥f ⋆ ∥ k .Thus, we may derive the following∥f ⋆ ∥ k -∥ f KT,λ ′ ∥ k < 0 and ∥f ⋆ ∥ k + ∥ f KT,λ ′ ∥ k > 1,which further implies the following inequality∥f ⋆ ∥ 2 k -∥ f KT,λ ′ ∥ 2 k = {∥f ⋆ ∥ k -∥ f KT,λ ′ ∥ k }{∥f ⋆ ∥ k + ∥ f KT,λ ′ ∥ k } ≤ ∥f ⋆ ∥ k -∥ f KT,λ ′ ∥ k . (80) Further writing f KT,λ ′ = f ⋆ + ∆ KT and noting that ∥ ∆ KT ∥ k -∥f ⋆ ∥ k ≤ ∥ f KT,λ ′ ∥ k holds throughtriangle inequality, we may further refine (80) as∥f ⋆ ∥ k -∥ f KT,λ ′ ∥ k ≤ 2∥f ⋆ ∥ k -∥ ∆ KT ∥ k ≤ 2 -∥ ∆ KT ∥ k ,so that the basic inequality in (78) reduces toEQUATION)We again introduce a technical event that controls the stochastic component of (81), which isB KT ≜ ∃g ∈ F \ B 2 (δ n ) ∩ {∥g∥ k ≥ 1} : (82) 1 nout nout i=1 v ′ i g(x ′ i ) > 4ξ n ∥g∥ nout + 2ξ 2 n ∥g∥ k + 1 4 ∥g∥ 2 nout ,for a star-shaped function class F ⊂ H.By triangle inequality, we have∥ ∆ KT ∥ k ≥ ∥ f KT,λ ′ ∥ k -∥f ⋆ ∥ k > 1,and on event E lower (⊂ E good ), we have ∥ ∆ KT ∥ 2 > δ n . Thus, we may apply g = ∆ KT to the event B c KT , and the resulting refined basic inequality is1 2 ∥ ∆ KT ∥ 2 nout ≤ 4ξ n ∥ ∆ KT ∥ nout + (2ξ 2 n -λ ′ )∥ ∆ KT ∥ k + 2λ ′ ≤ 4ξ n ∥ ∆ KT ∥ nout + 2λ ′ on the event B c KT ∩ E lowerwhere the second inequality is due to the assumption that λ ′ ≥ 2ξ 2 n . We apply the quadratic formula (specifically, if a, b ≥ 0 andx 2 -ax -b ≤ 0, then x ≤ a 2 + b) to obtain ∥ ∆ KT ∥ 2 nout ≤ 4c 2 0 ξ 2 n + 2λ ′ .Putting the pieces together, we have shown∥ ∆ KT ∥ 2 nout ≤ c(ξ 2 n + λ ′ ) on the event A c KT (ξ n ) ∩ B c KT ∩ E lower ,which is sufficient to recover (73).Recovering step (74) We now upgrade events{∥ ∆ KT ∥ 2 nout ≤ c(ξ 2 n + λ ′ )} =⇒ {∥ ∆ KT ∥ 2 n ≤ c ′ (ξ 2 n + λ ′ )}by exploiting the events E conc ∩ E KT,δ (subset of E good ) that were otherwise not used when recovering (73). For this end, the following result is a crucial ingredient, which shows that ∥•∥ nout and ∥•∥ n are essentially exchangeable with high-probability, Proposition 1 (Multiplicative guarantee for KT-COMPRESS++ with k RR ). Let C m ≜ 1/µ m and suppose δ n satisfies (63). Then on event E KT,δ ∩ E conc , where E KT,δ and E conc are defined in (23) and (69) respectively, we have(1 -4C m • η n,k )∥g∥ n ≤ ∥g∥ nout ≤ (1 + 4C m • η n,k )∥g∥ n (83)uniformly over all g ∈ H such that ∥g∥ 2 > δ n .See App. D.3 for the proof.An immediate consequence of Prop. 1 is that{∥ ∆ KT ∥ 2 nout ≤ c(ξ 2 n + λ ′ )} =⇒ {∥ ∆ KT ∥ 2 n ≤ c ′ (ξ 2 n + λ ′ )}on the event E KT,δ ∩ E conc ∩ E lower ,, which is sufficient to recover (74).Recovering step (75) Our last step is to show{∥ ∆ KT ∥ 2 n ≤ c ′ (ξ 2 n + λ ′ )} =⇒ {∥ ∆ KT ∥ 2 2 ≤ c ′′ (ξ 2 n + λ ′ )}.Such result can be immediately shown on the event E conc by observing that ∆ KT ∈ H, by our assumption that f ⋆ ∈ H and by the definitionf KT,λ ′ ∈ argmin f ∈H(k) L nout (f ) + λ ′ ∥f ∥ 2 k .It suffices to show the appropriate bounds for the following four probability termsP(A KT (ξ n )), P(B KT ), P(E c conc ), P(E c KT,δ ). Fix the shorthand B H ≜ ∥k∥ 2 ∞ R 2 < ∞.We know from [10] that P(E c KT,δ |S in ) ≤ δ and then we may apply [28, Thm. 14.1] to obtain a high probability statement,EQUATION)Now we present two Lemmas that bound the P(• | S in ) probability of events A KT (ξ n ) and B KT ,Lemma 5 (Controlling bad event when ∥ f KT,λ ′ ∥ k ≤ 2). Suppose u ≥ ξ n .Then for some constant c > 0,EQUATION)where σ = σ/∥f ⋆ ∥ k .See App. D.4 for the proof. Note that by plugging in ξ n into (85) results in a probability that depends on S in (as ξ n depends on S in ). By invoking the definition of ξ n , we may further refine the probability bound of A KT (ξ n ) byP(A KT (ξ n ) | S in ) ≤ δ + e -cnδ 2 n /B 2 H + e -cnδ 2 n / σ 2 (86) Lemma 6 (Controlling bad event when ∥ f KT,λ ′ ∥ k > 2). For some constants c, c ′ > 0,EQUATION)where σ = σ/∥f ⋆ ∥ k .See App. D.5 for the proof. It is notable that ξ n in the probability bound of (87) contains a term ε n defined in (59) that is a function of S in . Invoking the definition of ξ n , we observe the probability upper bound (87) can be refined toEQUATIONwhich does not depend on S in .Putting the pieces together, we have the following probability bound for some constants c, c ′ > 0,P(E c good | S in ) ≤ P(A KT (ξ n ) | S in ) + P(B KT | S in ) + P(E c conc | S in ) + P(E c KT,δ | S in ) (84)(86)(88) ≤ c δ + e -c ′ nδ 2 n /(B 2 H ∧ σ 2 ) thereby implying P(E c good ) ≤ c ′′ δ + e -c ′ nδ 2 n /(B 2 H ∧ σ 2 )for some constant c ′′ .Fix g ∈ H. Denote ⟨g, h⟩ = g(x)h(x)dx as the inner product in the L 2 sense. By Mercer's theorem [28, Cor. 12.26] , the k-norm of g has a basis expansionEQUATION)The assumption ∥g∥ 2 ≥ δ n implies that on the event E conc (69), we have1 2 δ n ≤ ∥g∥ 2 -1 2 δ n ≤ ∥g∥ n (90)Moreover, g must be a non-zero function. Note that g 2 ∈ H(k RR ) (see App. F.3). Thus, we may apply Lem. 12 to f 1 = f 2 = g and a = 1, b = 0 to obtainEQUATIONThe LHS can be expanded as∥g∥ 2 n -∥g∥ 2 nout = ∥g∥ n -∥g∥ nout • ∥g∥ n + ∥g∥ nout >0 by (90).Thus, we may rearrange (91) and combine with (89) to obtainEQUATIONOn event E conc , we haveEQUATIONThus, we haveEQUATION)Using (94) to refine (92), we have on event E KT,δ ∩ E conc :∥g∥ n -∥g∥ nout ≤ 4C m ∥g∥ n • η n,k .With some algebra, this implies with probability at least1 -δ -exp(-c ′ nδ 2 n /B 2 F ): (1 -4C m • η n,k )∥g∥ n ≤ ∥g∥ nout ≤ (1 + 4C m • η n,k )∥g∥ nuniformly over all non-zero g ∈ H such that ∥g∥ 2 > δ n .∥ f KT,λ ′ ∥ k ≤ 2Recall E KT,δ and E conc defined by ( 23) and (69). Also recall that E KT,δ ∩ E conc combined with the assumption ∥g∥ 2 ≥ δ n invokes the event (83). Our aim is to showEQUATIONso that we have a probability boundP(A KT (u)) ≤ P(E c KT,δ ) + P(E c conc ) + P(Z n (2u) ≥ 2u 2). The first RHS term can be bounded by δ (see (23) ). The second RHS term can bounded by (84). The third term can be bounded byP(Z n (2u) ≥ 2u 2 ) = P Z n (u) ≥ u 2 /2 + u 2 /2 (i) ≤ P Z n (u) ≥ uε n /2 + u 2 /2 (ii) ≤ e -nu 2 8 σ 2 ,where (i) follows from our assumption that u ≥ ε n and (ii) follows from applying generic concentration bounds on Z n (u) (see [28, Thm. 2.26, Eq. 13 .66]). Putting together the pieces yields our desired probability bound (85).Proof of claim (95). Consider the event A KT (u) ∩ E KT,δ ∩ E conc . The norm equivalence established on the event E KT,δ ∩ E conc in Prop. 1 is an important ingredient throughout.Let g ∈ H be the function that satisfies three conditions: ∥g∥ 2 ≥ δ n , ∥g∥ nout ≥ u, and1 nout nout i=1 v ′ i g(x ′ i ) ≥ 3∥g∥ nout u.Define the normalized function g = u • g/∥g∥ nout so that it satisfies ∥ g∥ nout = u and alsoEQUATION)By triangle inequality, the LHS of (96) can be further upper bounded byEQUATION)Recall the chosen g satisies ∥g∥ nout ≥ u. Observe thatv i g(x i ) (1) = (y i -f ⋆ (x i ))g(x i ) = -f ⋆ (x i )g(x i ) + y i g(x i ),so we may apply Lem. 12 with f 1 = f ⋆ , f 2 = g and a = -1, b = 1. Thus, on the event E KT,δ ∩ E conc , we haveEQUATIONThus, we may rearrange (97) and combine with ( 96) and ( 98) to obtain1 n n i=1 v i g(x i ) ≥ 3u 2 -u ∥g∥n out ∥g∥ k (∥f ⋆ ∥ k + 1) • η n,kNote that∥g∥ k ∥g∥n out = ∥g∥ k ∥g∥2 • ∥g∥2 ∥g∥n • ∥g∥n ∥g∥n out .We tackle each term in turn. First, ≤ 2 on event E conc ∩ E KT,δ . Taken together,EQUATION) As u ≥ ξ n ≥ 4 √ C m (∥f ⋆ ∥ k + 1)η n,kby assumption, we have therefore found g with norm ∥ g∥ nout = u satisfying1 n n i=1 v i g(x i ) ≥ 3u 2 -u 2 = 2u 2 .We may further show that∥ g∥ n = u ∥g∥n out ∥g∥ n ≤ u ∥g∥n ∥g∥n out ≤ u • 2 on event E conc ∩ E KT,δ ,where the last inequality follows from the fact that ∥g∥ 2 ≥ δ n and by applying (103). So we observe2u 2 ≤ 1 n n i=1 v i g(x i ) ≤ sup ∥ g∥n≤2u 1 n n i=1 v i g(x i ) = Z n (2u) D.5 Proof of Lem. 6: Controlling bad event when ∥ f KT,λ ′ ∥ k > 2Our aim is to show for any g ∈ ∂H with ∥g∥ k ≥ 1,1 nout nout i=1 v ′ i g(x ′ i ) ≤ 2ξ n ∥g∥ nout + 2ξ 2 n ∥g∥ k + 1 16 ∥g∥ 2 nout with high probability.Note that it is sufficient to prove our aim for g ∈ ∂H with ∥g∥ k = 1-by proving only for g with ∥g∥ k = 1, then for any h ∈ ∂H with ∥h∥ k ≥ 1, we may plug g = h/∥h∥ k into1 nout nout i=1 v ′ i g(x ′ i ) ≤ 2ξ n ∥g∥ nout + 2ξ 2 n + 1 16 ∥g∥ 2 nout (100)to recover the aim of interest. So without loss of generality, we show (100) for all g such that g ∈ ∂H and ∥g∥ k = 1.Let B KT denote the event where (100) is violated, i.e. there exists g ∈ ∂H with ∥g∥ k = 1 so thatEQUATIONWe prove the following set inclusion,EQUATIONwhere we know the RHS event of (102) has probability bounded by ce -nξ 2 n /(c ′ σ 2 ) which is proven in [28, Lem. 13.23] . So the set inclusion (102) implies a bound over the event B KT ,P(B KT ) ≤ P(E c KT,δ ) + P(E c conc ) + ce -nξ 2 n /(c ′ σ 2 ) ,where P(E c KT,δ ) ≤ δ by (23) and P(E c conc ) by ( 84). Choose g so that ∥g∥ k = 1 and (101) holds. Condition ∥g∥ k = 1 as well as the condition (89) resulting from a finite rank kernel k impliesδ n ≤ 1 ≤ ∥g∥ k ≤ √ C m ∥g∥ 2 . Invoke Prop. 1 for the choice of g that satisfies ∥g∥ 2 ≥ δ n / √ C m ≥ δ n ,so that on the event E KT,δ ∩ E conc , we have the following norm equivalence,EQUATIONThen we have the following chain of inequalities, which holds on event E conc ∩ E KT,δEQUATIONEQUATIONEQUATIONwhere step (i) follows from triangle inequality and step (ii) follows from our assumption (101) to bound the first term and our approximation guarantee (98) to bound the second term. By definition of ξ n , we have3 2 ξ n -2 √ C m (∥f ⋆ ∥ k + 1) • η n,k ≥ 2ξ n .Using this to refine (104), we have1 n n i=1 v i g(x i ) ≥ 2ξ n ∥g∥ n + 2ξ 2 n + 1 16 ∥g∥ 2 n ,which directly implies the inclusion (102) as desired.E Proof of Thm. 3: KT-KRR guarantee for infinite-dimensional RKHSWe state a more detailed version of the theorem:Theorem 5 (KT-KRR guarantee for infinite-dimensional RKHS, detailed). Assume f ⋆ ∈ H(k) and Assum. 1 is satisfied. If k is LOGGROWTH(α, β), then for some constant c (depending on d, α, β),EQUATIONwith probability at least 1 -2δ -2e- nδ 2 n c 1 (∥f ⋆ ∥ 2 k +σ 2 ) .If k is POLYGROWTH(α, β) with α ∈ (0, 2), then for some constant c (depending on d, α, β), f KT,λ ′ with λ = O(nEQUATION)with probability at least 1 -2δ -2e- nδ 2 n c 1 (∥f ⋆ ∥ 2 k +σ 2 ) .We state a generic result for infinite-dimensional RKHS that only depends on the Rademacher and Gaussian critical radii as well as the KT approximation term, all introduced in App. C. Theorem 6 (KT-KRR). Let f ⋆ ∈ H(k) and Assum. 1 is satisfied. Let δ n , ε n denote the solutions to (59), (60), respectively. Denote f KT,λ ′ with regularization parameter λ ′ ≥ 2η n,k , where η n,k is defined by (61). Then with probability at least 1 -2δ -2e- nδ 2 n c(∥f ⋆ ∥ 2 k +σ 2 ) -c 1 e -c2 n∥f ⋆ ∥ 2 k ε 2 n σ 2 , we have ∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ U full + U KT , where U full ≜ c ε 2 n + λ ′ [∥f ⋆ ∥ k + 1] 2 + cδ 2 n and U KT ≜ c • η n,k [∥f ⋆ ∥ k + 1] 2 .See App. F for the proof. The term U full follows from the excess risk bound of FULL-KRR f full,λ . The term U KT follows from our KT approximation. Clearly, the best rates are achieved when we choose λ = 2η n,k .The strategy for each setting is as follows:1. Bound the Gaussian critical radius (60) using [28, Cor. 13.18] , which reduces to finding ε > 0 satisfying the inequalityEQUATIONand μ1 ≥ μ2 ≥ . . . ≥ μn ≥ 0 are the eigenvalues of the normalized kernel matrix K/n, where K is defined by (6). 2. Bound the Rademacher critical radius (59) using [28, Cor. 14.5] , which reduces to solving the inequalityEQUATIONwhere (µ j ) ∞ j=1 are the eigenvalues of the k according to Mercer's theorem [28, Thm. 12.20] and b is the uniform bound on the function class. 3. Bound η n,k (61) using the covering number bound N (B d 2 (R in ), ϵ) from Assum. 3.In the sequel, we make use of the following notation. LetEQUATION= 1 + R in and L k (r) ≜ C d log 2 r β according to [15, Eq. 6] , where C d is the constant that appears in Assum. 3.We begin by solving (107). Lemma 7 (Critical Gaussian radius for POLYGROWTH kernels). Suppose Assum. 1 is satisfied and k is POLYGROWTH with α < 2 as defined by Assum. 3. Then the Gaussian critical radius satisfiesEQUATIONProof. [15, Cor. B.1] implies thatμj ≤ 4 L k (Rn) j-1 2 α for all j > L k (R n ) + 1Let k be the smallest integer such thatk > L k (R n ) + 1 and 4 L k (Rn) k-1 2 α ≤ ε 2 .By Assum. 1, R n is a constant, so the first inequality is easily satisfied for large enough nEQUATIONThen2 √ n n j=1 min{ε 2 , μj } ≤ 2 √ n kε 2 + n j=k+1 4 L k (Rn) j-1 2 α (i) ≤ 2 √ n kε 2 + 4L k (Rn) 2/α 2/α-1 k 1-2/α (110) ≤ 2 √ n 2 -α L k (R n )ε 2-α + 4•2 2-α L k (Rn) 2/α-1 ε 2-α ,where step (i) follows from the approximationn-1 j=k 4 L k (Rn) j 2 α ≤ 4L k (R n ) 2/α ∞ k t -2/α dt = 4L k (R n ) 2/α 1 2/α-1 k 1-2 α .To solve (107), it suffices to solve2c √ n 2 -α L k (R n )(1 + 16 2/α-1 )ε 2-α ≤ βε 2 =⇒ ε ≥ 2c β 2 2+α 2 -α L k (R n )(1 + 32α 2-α ) 1 2+α • n -1 2+α .Since ε n is the smallest such solution to (107) by definition, we have (109) as desired.We proceed to solve (108). Lemma 8. Suppose Assum. 1 is satisfied and k is POLYGROWTH with α < 2 as defined by Assum. 3. Then the Rademacher critical radius satisfiesε 2 n ≃ b 4 2+α 2 -α L k (R n )(1 + 32α 2-α ) 2 2+α n -2 2+α .Proof. Thus, we can solve the following inequality2 n n j=1 min{δ 2 , μj } ≤ 1 b δ 2 ,Following the same logic as in the proof of Lem. 7 but with β = 1/b yields the desired bound.Finally, it remains to bound (61). We haveη n,k = a 2 nout (2 + W k (n, n out , δ, R in , a nout )) ≤ a 2 nout (2 + log nout log(n/nout) δ • log 1 δ + log N k (B d 2 (R in ), a nout ) ) ≤ a 2 nout (2 + log nout log(n/nout) δ • log 1 δ + C d nout a α (R in + 1) β ) ≤ a 2 nout 2 + log nout log(n/nout) δ • log 1 δ + C d (Rin+1) β a α n α 2 out ≤ n α 2 -1 out • a 2 2 + log nout log(n/nout) δ • C d (Rin+1) β a αfor some universal positive constant c.In summary, there exists positive constants c 0 , c 1 , c 2 such thatε 2 n ≤ c 0 σ ∥f ⋆ ∥ k 4 2+α n -2 2+α δ 2 n ≤ c 1 b 4 2+α n -2 2+α η n,k ≤ c 2 a 2 n -2-α 2 out log n out Setting λ ′ = c 2 a 2 n -2-α 2 out log n out , we have ∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ c ε 2 n + λ ′ + η n,k • [∥f ⋆ ∥ k + 1] 2 + c ′ δ 2 n ≤ c∥f ⋆ ∥ 2 2+α k n -2 2+α + [∥f ⋆ ∥ k + 1] 2 n -2-α 2 out log n out + c ′ b 4 2+α n -2 2+α .We begin by solving (107). Lemma 9 (Critical Gaussian radius for LOGGROWTH kernels). Under Assum. 1 and LOGGROWTH version of Assum. 3, Gaussian critical radius satisfiesε 2 n ≃ σ 2 ∥f ⋆ ∥ 2 k log(2e• ∥f ⋆ ∥ k 4σ √ n) α n • L k (R n )C ′′ α for some constant C ′′α that only depends on α. where we ignore log-log factors.Proof. [15, Cor. B.1] implies thatμj ≤ 4 exp 2 -2 j-1 L k (Rn) 1 α for all j > L k (R n ) + 1Let k be the smallest integer such thatk > L k (R n ) + 1 and 4 exp 2 -2 j-1 L k (Rn) 1 α ≤ ε 2 .By Assum. 1, R n is a constant, so the first inequality is easily satisfied for large enoughEQUATIONConsider the following approximation: . Note that we can simplify the exponent byEQUATION-mL k (R n ) 1/α = -(L k (R n )/2) -1/α L k (R n ) 1/α = -2 1/α . Note that k ′ = k -1 ≥ L(R n ) = 2m -α. Thus, we can absorb the b α-1 m -α term in (113) into k and obtain the following boundn-1 ℓ=k 4 exp 2 -2 ℓ L k (Rn) 1 α ≤ C ′ α k ′ 2e ε -2 1/α, where C ′ α depends only on α. Plugging this bound into (112), we have2 √ n n j=1 min{ε 2 , μj } ≤ 2 √ n kε 2 + C ′ α k ε 2e 2 1/α (111) ≤ 2c √ n L k (R n ) log 2e ε α (ε 2 + C ′ α ε 2e 2 1/α ) ≤ 2c √ n L k (R n ) log 2e ε α C ′′ α ε 2 1/(1∨α)for some constant C ′′ α that only depends on α and universal positive constant c. To solve (107), it suffices to solve2c √ n L k (R n ) log 2e ε α C ′′ α ε 2 1/(1∨α) ≤ βε 2 ,which is implied by the looser bound1 β 2 • 4c 2 n • L k (R n )C ′′ α ≤ ε 2 log 2e ε -α .The solution to (107) (up to log-log factors) isε ≃ log(2e•β √ n) α/2 √ n 4c 2 β 2 • L k (R n )C ′′ α .We proceed to solve (108). Lemma 10 (Critical Gaussian radius for LOGGROWTH kernels). Under Assum. 1 and LOGGROWTH version of Assum. 3, the Rademacher critical radius satisfiesδ 2 n ≃ b 2 log( 2e b • √ n) α n • L k (R n )C ′′ α .Proof. Following the same logic as in the proof of Lem. 9 but with β = 1/b yields the desired bound.Finally, it remains to bound (61). We haveη n,k = a 2 nout (2 + W k (n, n out , δ, R in , a nout )) ≤ a 2 nout (2 + log nout log(n/nout) δ • log 1 δ + log N k (B d 2 (R in ), a nout ) ) ≤ a 2 nout (2 + log nout log(n/nout) δ • log 1 δ + C d log( enout a ) α (R in + 1) β )for some universal positive constant c.In summary, there exists universal positive constants c 0 , c 1 , c 2 such thatε 2 n ≤ c 0 σ 2 ∥f ⋆ ∥ 2 k log(2e• ∥f ⋆ ∥ k 4σ √ n) α n δ 2 n ≤ c 1 b 2 log( 2e b • √ n) α/2 n η n,k ≤ c 2 a nout log( enout a ) α/2 R β/2in .Setting λ ′ = 2c 2 a nout log( enout a ) α/2 , we have∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ c ε 2 n + λ ′ + η n,k • [∥f ⋆ ∥ k + 1] 2 + c ′ δ 2 n ≤ c log(2e• ∥f ⋆ ∥ k 4σ √ n) α n + [∥f ⋆ ∥ k + 1] 2 c nout log( enout a ) α/2 + c ′ b 2 log( 2e b • √ n) α/2 n .F Proof of Thm. 6: KT-KRR Our first goal is to bound the in-sample prediction error. We relate∥ f KT,λ ′ -f ⋆ ∥ 2 n to ∥ f full,λ -f ⋆ ∥ 2 n, where the latter quantity has well known properties from standard analyses of the KRR estimator (refer to [28] ). Note that regularization parameter λ ′ of KT based estimator f KT,λ ′ is independently chosen from the regularization parameter λ of the estimator based on original samples f full,λ . For f full,λ , we choose the regularization parameterEQUATIONwhich is known to yield optimal L 2 error rates.Define the main event of interest,E ≜ {∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ c(ε 2 n + δ 2 n + λ ′ + η n,k )[∥f ⋆ ∥ k + 1] 2 }.Our goal is to show E occurs with high probability. For that end, we introduce several additional events that are used throughout this proof.For some constant c >, define the event of an appealing in-sample prediction error of f KT,λ ′ ,E KT,n (t) ≜ ∥ f KT,λ ′ -f ⋆ ∥ 2 n ≤ c t 2 + λ ′ + η n,k • (∥f ⋆ ∥ k + 1) 2 for t ≥ ε n .where η n,k is defined in (61). Recall E KT,δ is the event where KT-COMPRESS++ succeeds as defined by (23) .Further as f ⋆ and f KT,λ ′ are both in {f ∈ H : ∥f ∥ k ≤ R}, we may deduce that all the functions under consideration satisfies ∥f ∥ ∞ ≤ ∥k∥ ∞ ∥f ∥ k ≤ ∥k∥ ∞ R where ∥k∥ ∞ < ∞. Accordingly, we define a uniform concentration event,EQUATIONEvent ( 115) is analogous to the event E conc previously defined in (69) when dealing with finite rank kernels.We first show thatEQUATIONNotice that almost surely we haveEQUATIONNext invoking the event E KT,n (ε n ∨ δ n ) along with (117), we have∥ f KT,λ ′ -f ⋆ ∥ 2 2 ≤ 2c[(ε n ∨ δ n ) 2 + λ ′ + η n,k ] • (∥f ⋆ ∥ k + 1) 2 + δ 2 n ≤ c(ε 2 n + δ 2 n + λ ′ + η n,k )[∥f ⋆ ∥ k + 1] 2 . which recovers the event of E.The remaining task is to show E is of high-probability, which amounts to showing events E KT,n (t) and E ′ conc are of high-probability by reflecting on (116). From [28, Thm. 14 .1], we may immediately deriveP(E ′ conc ) ≥ 1 -c 1 e -c2 nδ 2 n ∥k∥ 2 ∞ R 2for some constants c 1 , c 2 > 0.We further claim thatP(E KT,n (t) | S in ) ≥ 1 -δ -e -nt 2 c 0 σ 2 -c 1 e -c2 n∥f ⋆ ∥ 2 k t 2 σ 2 (118)for some constants c 0 , c 1 , c 2 > 0. Proof of claim (118) is deferred to App. F.1. Plugging in t = ε n ∨δ n into (118), and invoking inequality ε n ∨ δ n ≥ δ n so as to decouple the dependence on S in , we haveP(E KT,n (ε n ∨ δ n ) | S in ) ≥ 1 -δ -e - nδ 2 n c 0 σ 2 -c 1 e -c2 ∥f ⋆ ∥ 2 k nδ 2 n σ 2which further impliesP(E KT,n (ε n ∨ δ n )) ≥ 1 -δ -e - nδ 2 n c 0 σ 2 -c 1 e -c2 ∥f ⋆ ∥ 2 k nδ 2 n σ 2 .Putting the pieces together, for some constants c 0 , c 1 > 0, we haveEQUATIONOverall, ( 116) and (119) collectively yields the desired result.To prove claim (118), we introduce two new intermediary and technical events. For some positive constant c 0 , define the event 5 when in-sample prediction error of f full,λ is appealingEQUATIONThe second intermediary event, denoted as E ∆KT (t), is the intersection of ( 128) and ( 129), which we do not elaborate here due to its technical nature-event E ∆KT (t) plays an analogous role to A c KT ∩ B c KT defined in ( 82) and ( 95) respectively. Our goal here is two-folds: first is to show{E full,n (t) ∩ E KT,δ ∩ E ∆KT (t)} =⇒ E KT,n (t)and second is to prove the following boundP E full,n (t) ∩ E KT,δ ∩ E ∆KT (t) | S in ≥ 1 -δ -e -nt 2 c 0 σ 2 -c 1 e -c2 n∥f ⋆ ∥ 2 k t 2 σ 2, from which (118) follows. Note that Wainwright [28, Thm. 13.17] showP(E full,n (t)) ≥ 1 -c 1 e -c2 n∥f ⋆ ∥ 2 k t 2 σ 2for some constants c 1 , c 2 > 0 and that P(E KT,δ | S in ) ≥ 1 -δ. So it remains to bound the probability of event E ∆KT (t), which we show below.Given f , define the following quantitiesL n (f ) ≜ 1 n n i=1 (f 2 (x i ) -2f (x i )y i ) + 1 n n i=1 y 2 i and L nout (f ) ≜ 1 nout nout i=1 (f 2 (x ′ i ) -2f (x ′ i )y ′ i ) + 1 n n i=1 y 2 i .In the sequel, we repeatedly make use of the following fact: on event E KT,δ defined in (23) , we haveEQUATION)The claim of (121) is deferred to the end of this section. Given f , we can show with some algebra thatEQUATIONwhere Z ≜ (f (x 1 ) -f ⋆ (x 1 ), . . . , f (x n ) -f ⋆ (x n )) and ξ ≜ (ξ 1 , . . . , ξ n ) are vectors in R n . Define the shorthands∆ KT ≜ f KT,λ ′ -f ⋆ and ∆ full ≜ f full,λ -f ⋆ .In the sequel, we use the following shorthands:Z f ull ≜ ( ∆ full (x 1 ), . . . , ∆ full (x n )) and Z KT ≜ ( ∆ KT (x 1 ), . . . , ∆ KT (x n )).∥ f KT,λ ′ -f ⋆ ∥ 2 n . When ∥ ∆ KT ∥ n < t, we immediately have ∥ ∆ KT ∥ 2 n < t 2, which implies (118). Thus, we may assume that ∥ ∆ KT ∥ n ≥ t. Note that∥ ∆ KT ∥ 2 n (122) = L n ( f KT,λ ′ ) + 2 n ⟨Z KT , ξ⟩ -1 n n i=1 ξ 2 i = L n ( f full,λ ) + L n ( f KT,λ ′ ) -L n ( f full,λ ) + 2 n ⟨Z KT , ξ⟩ -1 n n i=1 ξ 2 i .Given the optimality of f full,λ on the objective (4), we haveL n ( f full,λ ) ≤ 1 n n i=1 ξ 2 i + λ ∥f ⋆ ∥ 2 k -∥ f full,λ ∥ 2 k ≤ 1 n n i=1 ξ 2 i + λ∥f ⋆ ∥ 2 k ,where the last inequality follows trivially from dropping the -∥ f full,λ ∥ 2 k term. Thus,EQUATION)Using standard arguments to bound the term 2 n ⟨Z KT , ξ⟩, we claim that on the event E ∆KT , we haveEQUATIONfor some positive constants c, c ′ , and that P(E ∆KT | S in ) ≥ 1 -e -nt 2 2σ 2 . We defer the proof of claim (124) to the end of this section. 124)-first observe the following decomposition:L n ( f KT,λ ′ ) -L n ( f full,λ ) in (L n ( f KT,λ ′ ) -L n ( f full,λ ) = L n ( f KT,λ ′ ) -L nout ( f KT,λ ′ ) + L nout ( f KT,λ ′ ) -L n ( f full,λ ) .On the event E KT,δ (23) , the first term in the display can be bounded byL n ( f KT,λ ′ ) -L nout ( f KT,λ ′ ) (121) ≤ (∥ f KT,λ ′ ∥ 2 k + 2) η n,k . Note that f KT,λ ′ is the solution to the following optimization problem, min f ∈H(k) L nout (f ) + λ ′ ∥f ∥ 2k , so the second term in the display can be bounded by the following basic inequalityL nout ( f KT,λ ′ ) + λ ′ ∥ f KT,λ ′ ∥ 2 k ≤ L nout ( f full,λ ) + λ ′ ∥ f full,λ ∥ 2 k so that L nout ( f KT,λ ′ ) -L nout ( f full,λ ) ≤ λ ′ ∥ f full,λ ∥ 2 k -∥ f KT,λ ′ ∥ 2 k .Thus, on event E KT,δ , we haveL n ( f KT,λ ′ ) -L n ( f full,λ ) ≤ (∥ f KT,λ ′ ∥ 2 k + 2) η n,k + λ ′ ∥ f full,λ ∥ 2 k -∥ f KT,λ ′ ∥ 2 k = 2η n,k + λ ′ ∥ f full,λ ∥ 2 k + {η n,k -λ ′ } • ∥ f KT,λ ′ ∥ 2 k (i) ≤ 2η n,k + λ ′ ∥ f full,λ ∥ 2 k (ii) ≤ λ ′ (∥ f full,λ ∥ 2 k + 1)where steps (i) and (ii) both follow from the fact that λ ′ ≥ 2η n,k (see assumptions in Thm. 6). To bound ∥ f full,λ ∥ 2 k , we use the following lemma: Lemma 11 (RKHS norm of f full,λ ). On event E full,n (120), we have the following boundEQUATIONfor some constant c 0 > 0.See App. F.2 for the proof. Putting things together, we haveL n ( f KT,λ ′ ) -L n ( f full,λ ) ≤ cλ ′ (∥f ⋆ ∥ k + 1) 2for some constant c-substituting this bound into (124) yields∥ f KT,λ ′ -f ⋆ ∥ 2 n ≤ ct 2 (∥f ⋆ ∥ k + 1) 2 + c ′ λ ′ (∥f ⋆ ∥ k + 1)2 , for some constants c, c ′ , which directly implies (118), i.e. implying{E full,n (t) ∩ E KT,δ ∩ E ∆KT (t)} =⇒ E KT,n (t).Proof of claim (121). Given f , define the functionEQUATIONand note thatL n (f ) -L nout (f ) = 1 n n i=1 ℓ ′ f (x i , y i ) -1 nout nout i=1 ℓ ′ f (x ′ i , y ′ i ).We first prove a generic technical lemma: Lemma 12 (KT-COMPRESS++ approximation bound using kEQUATIONlies in the RKHS H(k RR ). Moreover, on event E KT,δ , we haveP in ℓ f1,f2 -Q out ℓ f1,f2 ≤ (|a| • ∥f 1 ∥ k ∥f 2 ∥ k + |b| • ∥f 2 ∥ k ) • η n,k .uniformly for all non-zero f 1 , f 2 ∈ H(k).See App. F.3 for the proof. Applying the lemma with f 1 ≜ f, f 2 ≜ g and a = 1, b = -2, we haveP in ℓ ′ f -Q out ℓ ′ f ≤ (∥f ∥ 2 k + 2) • η n,k, which combined with the observation (127) yields the desired claim.Proof of claim (124). Case I: First suppose that ∥ ∆ KT ∥ k ≤ 1. Recall that ∥ ∆ KT ∥ n ≥ t ≥ ε n by assumption. Thus, we may apply [28, Lem. 13 .12] to obtain1 n ⟨Z KT , ξ⟩ ≤ 2∥ ∆ KT ∥ n t w.p. at least 1 -e -nt 2 2σ 2 (128)Plugging the above bound into (123), we have with probability at least 1 -e -nt 2 2σ 2 :∥ ∆ KT ∥ 2 n ≤ 4∥ ∆ KT ∥ n t + λ∥f ⋆ ∥ 2 k + L n ( f KT,λ ′ ) -L n ( f full,λ ) .We can solve for ∥ ∆ KT ∥ n using the quadratic formula. Specifically, if a, b ≥ 0 and x 2 -ax -b ≤ 0, then x ≤ a + √ b. Thus, we have with probability at least 1 -e -nε 2 n 2σ 2 :∥ ∆ KT ∥ n ≤ a + √ b, where a ≜ 4t and b ≜ λ∥f ⋆ ∥ 2 k + L n ( f KT,λ ′ ) -L n ( f full,λ ) .Using the fact that (a + √ b) 2 ≤ 2a 2 + 2b, we have with probability at least 1 -e -nt 2 2σ 2 :∥ f KT,λ ′ -f ⋆ ∥ 2 n ≤ 32t 2 + 2λ∥f ⋆ ∥ 2 k + 2 L n ( f KT,λ ′ ) -L n ( f full,λ ) (114) ≤ ct 2 (∥f ⋆ ∥ k + 1) 2 + 2 L n ( f KT,λ ′ ) -L n ( f full,λ )Case II: Otherwise, we may assume that ∥ ∆ KT ∥ k > 1. Now we apply [28, Thm. 13.23] to obtainEQUATIONfor some universal positive constants c 1 , c 2 . Plugging the above bound into (123) and collecting terms, we have with probability at least 1 -c 1 e-nt 2 c 2 σ 2 : 7 8 ∥ ∆ KT ∥ 2 n ≤ 4t∥ ∆ KT ∥ n + 4t 2 ∥ ∆ KT ∥ k + λ∥f ⋆ ∥ 2 k + L n ( f KT,λ ′ ) -L n ( f full,λ ) .Solving for ∥ ∆ KT ∥ n using the quadratic formula, we have with probability at least 1 -c 1 e-nt 2 c 2 σ 2 : ∥ ∆ KT ∥ n ≤ a + √ b, where a ≜ 32 7 t and b ≜ 32 7 t 2 ∥ ∆ KT ∥ 2 k + 8 7 λ∥f ⋆ ∥ 2 k + 8 7 L n ( f KT,λ ′ ) -L n ( f full,λ ) .Using the fact that (a + √ b) 2 ≤ 2a 2 + 2b, we have with probability at least 1 -c 1 e-nt 2 c 2 σ 2 : ∥ f KT,λ ′ -f ⋆ ∥ 2 n ≤ 42t 2 + 10t 2 ∥ ∆ KT ∥ 2 k + 2.3λ∥f ⋆ ∥ 2 k + 2.3 L n ( f KT,λ ′ ) -L n ( f full,λ ) (i) ≤ 42t 2 + 10t 2 ∥ ∆ KT ∥ 2 k + 4.6t 2 ∥f ⋆ ∥ 2 k + 2.3 L n ( f KT,λ ′ ) -L n ( f full,λ ) (114),(125) ≤ c 3 t 2 (∥f ⋆ ∥ k + 1) 2 + c 4 L n ( f KT,λ ′ ) -L n ( f full,λ )for some positive constants c 3 , c 4 , where step (i) follows from that fact that λ = 2ε Given the optimality of f full,λ on the objective (4), we have the following basic inequalityL n ( f full,λ ) + λ∥ f full,λ ∥ 2 k ≤ 1 n n i=1 ξ 2 i + λ∥f ⋆ ∥ 2 k =⇒ ∥ f full,λ ∥ 2 k ≤ ∥f ⋆ ∥ 2 k + 1 λ 1 n n i=1 ξ 2 i -L n ( f full,λ ) .Since ∥ f full,λ -f ⋆ ∥ 2 n ≥ 0, we also have the trivial lower boundEQUATIONand it remains to bound 2 n ⟨Z f ull , ξ⟩. Case I: First, suppose that ∥ ∆ full ∥ k > 1. Then we may apply [28, Lem. 13.23] to obtain1 n ⟨Z f ull , ξ⟩ ≤ 2ε n ∥ ∆ full ∥ n + 2ε 2 n ∥ ∆ full ∥ k + 1 16 ∥ ∆ full ∥ 2 n w.p. at least 1 -c 1 e - nε 2 n c 2 σ 2 .Combining this bound with (130), we have with probability at least 1 -c 1 e- nε 2 n c 2 σ 2 : ∥ f full,λ ∥ 2 k ≤ ∥f ⋆ ∥ 2 k + 2ε 2 n λ ∥ ∆ full ∥ k + 2 λ 2ε n ∥ ∆ full ∥ n + 1 16 ∥ ∆ full ∥ 2 n (i) ≤ ∥f ⋆ ∥ 2 k + 2ε 2 n λ (∥ f full,λ ∥ k + ∥f ⋆ ∥ k ) + 2 λ 2ε n ∥ ∆ full ∥ n + 1 16 ∥ ∆ full ∥ 2 n (114) = ∥f ⋆ ∥ 2 k + ∥ f full,λ ∥ k + ∥f ⋆ ∥ k + 2 λ 2ε n ∥ ∆ full ∥ n + 1 16 ∥ ∆ full ∥ 2 n ,where step (i) follows from triangle inequality. Solving for ∥ f full,λ ∥ k using the quadratic formula, we have∥ f full,λ ∥ 2 k ≤ 2 + ∥f ⋆ ∥ 2 k + ∥f ⋆ ∥ k + 2 λ 2ε n ∥ ∆ full ∥ n + 1 16 ∥ ∆ full ∥ 2 n .On the event E full,n (120), we have ∥ ∆ full ∥ n ≤ c∥f ⋆ ∥ k ε n for some positive constant c, which implies the claimed bound (125) after some algebra.Case Combining this bound with (130), we obtain∥ f full,λ ∥ 2 k ≤ ∥f ⋆ ∥ 2 k + 1 λ ε 2 n (114) = ∥f ⋆ ∥ 2 k + 1 2, which immediately implies the claimed bound (125).Case II(b): Finally, assume ∥ ∆ full ∥ k ≤ 1 and ∥ ∆ full ∥ n > ε n . Applying [28, Lem. 13 .12] with u = ε n , we have1 n ⟨Z f ull , ξ⟩ ≤ 2ε 2 w.p. at least 1 -e -nε 2 n 2σ 2 .Combining this bound with (130), we obtainEQUATION= ∥f ⋆ ∥ 2 k + 2, which immediately implies the claimed bound (125).F.3 Proof of Lem. 12: KT-COMPRESS++ approximation bound using k RR By Grünewälder [12, Lem. 4] , ℓ f1,f2 lies in the RKHS H(k RR ), which is a direct sum of two RKHS:H(k RR ) = H(k 1 ) ⊕ H(k 2 ),where k 1 , k 2 : Z × Z → R are the kernels defined byk 1 ((x 1 , y 1 ), (x 2 , y 2 )) ≜ k 2 (x 1 , x 2 ) and k 2 ((x 1 , y 1 ), (x 2 , y 2 )) ≜ k(x 1 , x 2 ) • y 1 y 2 . Applying Lem. 2 with Z = X × Y, k ALG = k RR , and ϵ ⋆ = (∥k∥ 1/2 ∞,in +Ymax) 2 nout, yields the following bound on event E KT,δ (23) :sup h∈H(kRR): ∥h∥ k RR ≤1 |(P in -Q out )h| ≤ 2ϵ ⋆ + ∥kRR∥ 1/2 ∞,in nout • W kRR (n, n out , δ, R in , ϵ ⋆ ).We claim thatEQUATIONfor some positive constant c, where N † kRR is the cardinality of the cover ofB † kRR ≜ ℓ ′ f /∥ℓ ′ f ∥ kRR : f ∈ H(k)for ℓ ′ f defined by (126). Proof of the claims (131) and (132) are deferred to the end of this section. By definition of W kRR , we haveW kRR (n, n out , δ, R in , ϵ ⋆ ) (28) ≤ √ c • W k (n, n out , δ, R in , ∥k∥ 1/2 ∞,in +Ymax nout ) ≜ W ′ k .On event E KT,δ , we have sup h∈H(kRR):EQUATIONwhere step (i) follows from the fact that (∥k∥1/2 ∞,in + Y max ) 2 ≤ 2(∥k∥ ∞,in + Y 2 max ).Since f 1 , f 2 are non-zero, we have ∥ℓ f1,f2 ∥ k > 0. Thus, the function h ≜ ℓ f /∥ℓ f ∥ kRR ∈ H(k RR ) is well-defined and satisfies ∥h∥ kRR = 1. Applying (133), we obtain|P in h -Q out h| ≤ ∥k∥ ∞,in +Y 2 max nout • (2 + W ′ k ) on event E KT,δ .Multiplying both sides by ∥ℓ f ∥ kRR and noting that∥ℓ f,g ∥ 2 kRR = ∥a • f 1 f 2 ∥ 2 H⊙H + ∥b • f 2 ⊗ ⟨•, 1⟩ R ∥ 2 H⊗R ≤ a 2 ∥f 1 ∥ 2 k ∥f 2 ∥ 2 k + b 2 ∥f 2 ∥ 2 k ≤ (|a| • ∥f 1 ∥ k ∥f 2 ∥ k + |b| • ∥f 2 ∥ k ) 2 ,we have on event E KT,δ , ).P in ℓ f1,f2 -Q out ℓ f1,f2 ≤ (|a| • ∥f 1 ∥ k ∥f 2 ∥ k + |b| • ∥f 2 ∥ k ) • ∥k∥ ∞,inDefine a kernel on R by k R (y 1 , y 2 ) ≜ y 1 y 2 . When sup y∈(Sin)y |y| ≤ Y max , we haveN k R ([-Y max , Y max ], ϵ) = O(Y 2 max /ϵ) for ϵ > 0.Similarly, note that log N † k2 (S in , ϵ ⋆ ) ≤ log N k (S in , ϵ ⋆ /(∥k∥ Justification: Abstract and introduction give clear outlines on our contributions, and we present our contributions in the main text accordingly. We have also included pointers in the introduction that would link to the referred main text containing specific contributions.1/2 ∞ + ∥k R ∥ 1/2 ∞ )) + log N k R (S in , ϵ ⋆ /(∥k∥Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.We note that while KRR approach(15) does require k to be reproducing, the NW approach (3) in full generality is valid even when k is a not a valid reproducing kernel.In the sequel, we use "KT" and "KT-COMPRESS++" interchangeably since the underlying algorithm (kernel halving [10, Alg. 1a]) and associated approximation guarantees are the same up to small constant factors.One might expect the ridge regression kernel to include a term that accounts for y 2 . However, the generalization bounds turn out to be essentially the same regardless of whether we include this term when defining kRR.Note that when k is finite-rank, this condition is automatically satisfied.Since the input points in Sin are fixed, the randomness in f full,λ originates entirely from the randomness of the noise variables ξ. 


Decision: This paper introduces a meta-algorithm based on Kernel Thinning for non-parametric regression. The authors apply their meta-algorithm to Nadaraya-Watson regression and Kernel Ridge Regression and prove computational speedups and statistical efficiency in both cases. Overall the paper could benefit from more experimental evidence. In the original set of results there was not much empirical advantage over RPCholesky. Moreover the paper was missing comparisons to other nonparametric regression methods, there were no experiments on large datasets, and the results were limited to kernels with fast decaying tails (excluding Matern kernels). The authors have clarified all of these points in the rebuttal, including new experimental results and a claim (without proof) that their results extend to heavy tailed kernels. I believe that this paper should be accepted for publication despite these omissions. Nevertheless, I urge the authors to include all new results in the camera-ready, as well as a detailed proof for non sub-Gaussian kernels.

Reviews: review: paper_summary: paper_summary: After the rebuttal, I have updated my score from 3 to 5. 

----
This paper provides a meta-algorithm based on Kernel Thinning for non-parametric regression, in particular, the Nadaraya-Watson (NW) regression, and the Kernel Ridge Regression (KRR). 

The idea is to run NW or KRR on a thinned coreset by KT. The core of the method is to choose a suitable kernel for which the KT assumptions are satisfied. The paper introduces to meta-kernels that takes a base kernel of features x and produces the kernel over (x,y) suitable for regression tasks. 

The paper provides the theoretical guarantees for kernel thinned NW and kernel thinned KRR. 

In empirical study, the paper includes simulation results and real data results which demonstrate the proposed method's strength over full computation and standard thnnning in accuracy and time; however, the SOTA RPCholesky method seems to have a better empirical performance than kernel thinned KRR.
strength_weakness: strength_weakness: - The empirical advantage of the proposed kernel thinned KRR over RPCholesky is not clear. In fact, from figure 5, RPCholesky has slighly higher training time compared to kernel thinned KRR, however, a much lower test MSE. 

- Some technical background and discussion can be made more clear. See questions. 

- Writing can be improved. See questions.
questions: questions: - There is not an algorithmic description for Kernel Thinning algorithm. 

- Eq 4 and 6: why do you choose the kernel in this way? In particular, can you comment on the part $k(x_1,x_2) y_1, y_2$. Does the performance of the kernel rely on the relationship of $y$ and $x$, and if so, how? Also, in eq 6, is this a typo $k^2(x_1,x_2)^2$? 

- Notation suggestion: instead of $S_{KT}$ for KT-NW, could use $S_{NW}$ to match the use of $K_RR$ in KT-KRR part. 

- Sec 4.1: are the kernels listed in 1.2.3 used as the base kernel to construct the kernel for thinning, or they are the kernel for thinning?

- Fig 2: Right panel: why not include RPCholesky?

- Fig 2: why si the full not performing the best?

- Fig 3 and Fig 4: why not use the function in (11)? 

- Fig 3 and Fig 4: what is the unit of y axis? Are they all in the same unit?

- Sec 4.2: what is the conclusion from Fig 5?

- Sec 4.3: line 283-284: "In this example (and ....)" incomplete sentence. 

- It'd be good to include Fig 5 in the main text instead of appendix.
limitations: limitations: - KT-KRR does not demonstrate an empirical advantage over RPCholesky. Missing a discussion on this point. 

- Missing a discussion on potential failure cases of the method, though there is  a discussion on kernel mis-specification and potential way to deal with it.
paper_summary: paper_summary: The authors speed up two non-parametric regression estimators, Nadaraya-Watson (NW) and Kernel Ridge Regression (KRR), by using a kernel thinning (KT) technique to compress the input data in a way that preserves important statistical properties. They include a theoretical analysis proving that that KT-based regression estimators are computational more efficient and also exhibit improved statistical efficiency over i.i.d. subsampling of the training data. They also include empirical validation of their method on simulated and real-world data.
strength_weakness: strength_weakness: The empirical analysis is limited. The shortcomings on the real-world data are not adequately addressed (see also limitations).
questions: questions: I have no particular questions.
limitations: limitations: The comparison with other large-scale kernel methods is limited to RPCholesky; a more comprehensive comparison and overview would be nice.

Performance/utility of the proposed method on real-world data remains unclear. The results obtained on real-world dataset in section 4.2 are not convincing, and the authors acknowledge the problem of kernel mis-specification but do not elaborate or provide any resolution. Btw, there is an incomplete sentence in lines 283-284.
paper_summary: paper_summary: This paper applies the kernel thinning approach to non-parametric regression, with two kinds of estimators (NW and KRR) considered. This approach speeds up the computational efficiency by using carefully chosen coresets for approximation. Theoretical results on the approximation error are established. Numerical experiments include synthetic data and the well-known California housing dataset.
strength_weakness: strength_weakness: 1. Despite the well developed theory, I am not fully convinced yet by the significance/contribution of this method - I am willing to stand corrected if more comparisons can be made with other existing methods for efficient non-parametric regression. 
2. The application are only to toy examples or standard small datasets - e.g. for the Californian housing dataset with n = 20,640, doing non-parametric regression without thinning also seems fine.
3. In the experiments, the baseline is chosen to be NW / KRR with naive thinning - some more competitive, state-of-the-art baselines could have been chosen.
questions: questions: 1. In the beginning of Sec 3, the NW and KRR estimators are restricted to kernels taking the form of (4) and (6). What are the specific reasons to make these restrictions? And to see how restrictive/general this is, could you give some examples of commonly used kernels that fall under this category?
2. The theory makes the sub-Gaussian tail decay assumption on the kernel, which would exclude many commonly used kernels such as the Matern kernel. How would the theory work out, and how would the rates change if this assumption is generalized to sub-Weibull tail decay, or say a bound on the Orlicz-norm?
3. The method (or at least the theory) is developed for one-dimensional response data. How could this be generalized to multi-dimensional response data, where issues with computational efficiency becomes more significant?
4. In the discussion section, it is mentioned that the rates on the excess risk could be sub-optimal. Do you have a sense of what rate could be minimax optimal?
limitations: limitations: See weaknesses and questions.

Comments: title: Further response
Response to further comments
followup questions
comment: Thank you to the authors for the detailed reply. Most of my concerns and questions have been addressed, and I will increase my rating accordingly.
Thank you for your detailed response! Best of luck to the authors.
We thank Reviewer QHUc for your comments, which we believe have been addressed in our response. Please let us know if you have any other questions that we can address.​⬤
Thank you for your follow-up questions.

1. The values in Table A and B are test loss while the values in Figure 2 are excess risk. Hence a test loss of 0.9980 becomes an excess risk = test loss - 1 = -0.0020, which is not a valid value on log scale; hence the missing marker for RPCholesky in Figure 2.

2. Since the Tables are test loss, a lower value is better. So the opposite trend to what you wrote is true. In particular, in Table A, RPCholesky performs **worse** than our proposed method. On the other hand, RPCholesky is tends to do better than our proposed method In Table B. Overall the main takeaways from Table A/B (and Figure 2) is that **our Kernel Thinning based estimator**
   - _improves upon standard thinning for both NW and KRR estimators across all coreset sizes_
   - _improves upon RPCholesky for NW estimator across all coreset sizes_
   - _improves upon RPCholesy for small coreset sizes for KRR, and remains competitive (but a bit worse) in large coreset sizes for KRR._

RPCholesky’s superior performance over KT for KRR is not too surprising since RPCholesky a spectral method tuned to perform better in such settings, while NW estimator does not benefit immediately from an improved spectral approximation.

We remark that these results are consistent with our original motivation of designing a unified framework for speeding up general non-parametric regression estimators using a single recipe.
Thank you for your clarification. 

- I don't get why the loss of 0.9980 becomes undefined when y is in log2 scale? 

- Additionally, can you provide some reasoning between the performance comparison of Cholesky and the proposed method? In table A, cholesky performs better than the proposed for most of cases, but the trend is reversed in table B.
1. Thank you for the catch. RPCholesky-KRR achieves an estimated test loss of 0.9980 when $n=2^{14}$  (see Table B below)—close to the Bayes optimal loss of 1. We have estimated the population test loss by sampling 10,000 points from our data generating distribution (as a proxy for integrating), so the estimated test loss can actually be less than 1 with some probability. When we use log-scaling on the y-axis, negative values of excess risk become undefined, so they don't appear on the plot. Note that Full-NW and Full-KRR for $n=2^{14}$ also face this issue. We will make sure to add a note about this in the revision.
2. In Fig. 2, Full-NW and Full-KRR are the black box plots (always the first in each group when reading from left to right). Full is deterministic, so the box plots appear as lines. _Note that Full-NW and Full-KRR indeed perform the best across all $n$._ For additional clarity, we have attached tables of the test loss (mean ± standard deviation across 100 trials). (Subtract each value by 1 to get the excess risk values used in Fig. 2.)

**Table A: Test loss for NW estimators in Figure 2a of the attached PDF**
|                                                       | $n=2^{8}$       | $n=2^{10}$      | $n=2^{12}$      | $n=2^{14}$      |
|:------------------------------------------------------|:----------------|:----------------|:----------------|:----------------|
| Full                                                  | **1.1423 ± 0.0000** | **1.0494 ± 0.0000** | **1.0189 ± 0.0000** | **0.9997 ± 0.0000** |
| RPCholesky                                            | 3.0737 ± 0.2548 | 2.8583 ± 0.3150 | 2.0772 ± 0.2920 | 1.4720 ± 0.1026 |
| ST                                                    | 3.1181 ± 0.1813 | 2.8844 ± 0.3063 | 2.4720 ± 0.3058 | 1.7632 ± 0.2364 |
| KT w/ $\mathbf{k}_{\mathrm{RR}}((x_1,y_1),(x_2,y_2))$ | 2.9203 ± 0.4216 | 2.2230 ± 0.4160 | 1.3952 ± 0.0955 | 1.1403 ± 0.0274 |
| KT w/ $\mathbf{k}(x_1,x_2)$                           | 3.1040 ± 0.1928 | 2.4149 ± 0.2585 | 1.6458 ± 0.1386 | 1.3679 ± 0.0705 |
| KT w/ $\mathbf{k}((x_1\oplus x_2), (y_1\oplus y_2))$  | 3.0814 ± 0.1648 | 2.6633 ± 0.3603 | 1.7260 ± 0.1791 | 1.4282 ± 0.0825 |
| KT w/ $\mathbf{k}_{\mathrm{NW}}((x_1,y_1),(x_2,y_2))$ | 2.9838 ± 0.2056 | 2.2423 ± 0.3736 | 1.3799 ± 0.0808 | 1.1363 ± 0.0267 

$\newline$

**Table B: Test loss for KRR estimators in Figure 2b of the attached PDF**
|                                                       | $n=2^{8}$       | $n=2^{10}$      | $n=2^{12}$      | $n=2^{14}$      |
|:------------------------------------------------------|:----------------|:----------------|:----------------|:----------------|
| Full                                                  | **1.3750 ± 0.0000** | **1.0584 ± 0.0000** | **1.0121 ± 0.0000** | **0.9981 ± 0.0000** |
| RPCholesky                                            | 2.7112 ± 0.3636 | 1.7153 ± 0.3839 | 1.0112 ± 0.0029 | 0.9980 ± 0.0000 |
| ST                                                    | 2.9433 ± 0.3129 | 2.5855 ± 0.3570 | 2.1328 ± 0.3374 | 1.5634 ± 0.1437 |
| KT w/ $\mathbf{k}_{\mathrm{NW}}((x_1,y_1),(x_2,y_2))$ | 2.3068 ± 0.4381 | 1.9520 ± 0.2783 | 1.4859 ± 0.1146 | 1.2394 ± 0.0538 |
| KT w/ $\mathbf{k}(x_1,x_2)$                           | 2.8103 ± 0.1677 | 2.4297 ± 0.2559 | 1.7077 ± 0.1433 | 1.4974 ± 0.0885 |
| KT w/ $\mathbf{k}((x_1\oplus x_2), (y_1\oplus y_2))$  | 3.0312 ± 0.2431 | 2.6189 ± 0.3281 | 2.1086 ± 0.3178 | 1.4590 ± 0.1206 |
| KT w/ $\mathbf{k}_{\mathrm{RR}}((x_1,y_1),(x_2,y_2))$ | 2.0422 ± 0.3406 | 1.8856 ± 0.2575 | 1.4336 ± 0.1142 | 1.2316 ± 0.0499 |
I would like to thank the authors for their detailed reply. Most of my concerns are addressed and I have some followup questions:

- In Figure 2b in the attached pdf, why is RPCholesky not included at n=2^14? It is because it is not feasible to run RPCholesky at that point?

- Can you also address this question: Fig 2: why is the full not performing the best?