ParK: Sound and Efficient Kernel Ridge Regression by Feature Space Partitions
Authors: Luigi Carratino
Stefano Vigogna
Daniele Calandriello
Lorenzo Rosasco
Venue: NeurIPS 2021 Poster (2021)
Keywords: Large-scale kernel methods
kernel ridge regression
random projections
partitions

Abstract: We introduce ParK, a new large-scale solver for kernel ridge regression. Our approach combines partitioning with random projections and iterative optimization to reduce space and time complexity while provably maintaining the same statistical accuracy. In particular, constructing suitable partitions directly in the feature space rather than in the input space, we promote orthogonality between the local estimators, thus ensuring that key quantities such as local effective dimension and bias remain under control. We characterize the statistical-computational tradeoff of our model, and demonstrate the effectiveness of our method by numerical experiments on large-scale datasets.

ParK: Sound and Efficient Kernel Ridge Regression by Feature Space Partitions
Luigi Carratino" MaLGa DIBRIS, University of Genova carraaiirrati unige. bagge.1 it
Stefano Vigogna"
MaLGa DIBRIS, University of Genova vigognaédibris. unige.it it
Daniele Calandriello
Lorenzo Rosasco
DeepMind Paris
MaLGa DIBRIS, University of Genova IIT. CBMM MIT
deallandrielll@googlee. com
1rosascodmit.edu
Abstract
We introduce ParK, new large-scale solver for kernel ridge regression. Our approach combines partitioning with random projections and iterative optimization to reduce space and time complexity while provably maintaining the same statistical accuracy. In particular, constructing suitable partitions directly in the feature space rather than in the input space, promote orthogonality between the local estimators, thus ensuring that key quantities such as local effective dimension and bias remain under control. We characterize the statistical statisticalccomputationaaa tradeoff of our model, and demonstrate the effectiveness of our method by numerical experiments on large-scale datasets.
1 Introduction
The development of provably accurate and efficient algorithms for learning is key to tackle modern large-scale applications. Kernel methods [31, 32] provide natural ground to develop this research direction. On the one hand they have sound statistical guarantees [7, 32, 33]. but on the other hand their basic implementations
limited to sample size of only few tens of thousands of points [32, Chapter 11]. Recent years have witnessed a growing literature introducing algorithmic solutions to improve efficiency, but also theoretical guarantees that quantify how accuracy is affected.
We next recall few lines of work relevant to our study. first line of work is based on exploiting ideas from optimization and numerical analysis. This includes for example gradient methods [37], as well as their accelerated [3], stochastic [10]. preconditioned [13] and distributed [27] variants. second line of work is based on using randomized approaches to reduce the size of the problem to be solved. This includes Nyström approximations [36], random features [25] and more generally sketching methods [1]. The theoretical properties of these methods have been recently characterized in terms sharp statistical bounds [28, 30]. Finally, a third line of work considers different partitioning strategies to divide the estimation step in smaller subproblems. This approach is based splitting the input space in regions where local estimators are defined [23, 35, 34, 11, 24, 5]. In this context, allowing the estimation of larger classes of functions. Another form of
the emphasis is typically partitioning, called divide-and-con -ovide-and-conquer, is instead based on randomly splitting the training data to then obtain global estimator by averaging [38, 20, 14]. In this approach, the focus is primarily computational saving. Notably, a number of works have considered combinations of these ideas, see for example [6, 29, 8, 24, 19].
equal contribution
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

In this paper we propose and study a local kernel algorithm, called ParK, combining partitioning with iterative optimization and sketching. Our goal is to provide an efficient and accurate approximation to a global kernel ridge regression estimator. The main novelty in ParK is in the form of the considered partition, that we define in the feature space, rather than in the input space as in traditional partitioning methods. This allows to promote orthogonality between the local estimators, and thus to control the local effective dimension and the local bias. Given a partition, local kernel ridge estimators are computed using sketching and preconditioned conjugate gradient iterations [29]. From theoretical point of view, our main contribution characterizing the statistical properties of ParK, terms of conditions the partition and the choice of the ayper-parametert Borrowing ideas from subspace clustering [12]. we show that the minimal angle between suitable subspaces induced by the partition plays crucial role. Indeed, our analysis shows that, if such an angle is sufficiently large, ParK can achieve the same accuracy as global kernel ridge regression estimators, with only a fraction of computations, Our theoretical results complemented with numerical experiments very large datasets, which show that ParK can indeed provide excellent performances, on par and often better than the best available large-scale kernel methods.
The rest of the paper is organized as follows. In Section we state the problem and recall the basics of kernel ridge regression. In Section 3 we illustrate our algorithm. In Section 4 we analyze the prediction eπror of our method. In Section 5 we present the results of our numerical experiments. In Section 6 we draw some conclusions and report the main limitations of our work. Additional proofs and details are collected in Appendix A.
Background
Let (Iz, Us) with i E [n] = (1,...,,7 n be pairs of points in X' x V where X' C Rd with d E N and V S R. We assume the relation between input points 2; and output points ys to be determined by the noisy evaluations of an unknown function S. X + y as
y = f. (zx) + E i [n].
(1)
Based on the samples (t), we want estimate the function f., searching for solutions in suitable hypothesis space H as detailed below.
Let H. be reproducing kernel Hilbert space (RKHS), that is, a Hilbert space of functions with inner
product 2y and symmetric positive definite kernel K : X' x X + IR such that K2 = K K(x,-) E H and ((x (J, Ka)y for all ʃ E H. E X. We recall that, for every RKHS H. there exist Hilbert feature space ff and a feature map Φ X + F such that K(x, = (o(x), (((x))> for all X. The feature map not unique; in particular, one may take, as we do in all that follows, I, F= H and 6(x)= Kg. in which case H. span o(X'), where ((X) = (ő(x) x E X}.
Kernel ridge regression (KRR) corresponds minimizing
Σ (f(x.) vil + 11
min JEH
where À 0 and |l (f. fyM. By the representer theorem [31], the (unique) solution to problem (2) can be written as
Aa(x) Σ o ((;,a))) a a= (Kn Ann)
where a = (a......................... anl Y [yy....., yn JT R", and K,, E Raxx is the kernel matrix defined by (Kn),, ((x,, .,) for i, E [n]. As consequence, the estimator (3) can be derived restricting the minimization problem (2) to the finite-dimensional subspace H. = span( Computing (3) for large is prohibitively expensive, space and time complexities are, respectively, Onn²) and O(n²). The goal of this paper is to provide an algorithm to compute an efficient approximation to (3).
i E [n]]-
3 Algorithm
Our method combines diverse techniques, including partitioning, sketching and preconditioned iterative optimization. We begin focusing on partitioning. While classical partitioning methods

---

construct partitions in the input space, the main novelty of our approach is that we construct partitions the feature space. Note that, the case universal kernel on compact imput space, every feature map is injective [32, Lemma 4.53], hence every partition of the input space defines corresponding partition of the feature space. Thus, we may see our approach as generalization of classical input space partitioning approaches. As will become apparent from our analysis, the performance of partitioned kernel estimator depends crucially on two main quantities: the local biases and the local effective dimensions. Since both quantities are strictly related to the RKHS of choice, constructing partitions feature space allows for more direct control. In particular, promoting orthogonality in the RKHS metric will generate feature space partitions which tend minimize both the local biases and the local effective dimensions, In the next section we start discussing how such partitions can be defined.
3.1 Learning feature space partitions
For Q E N, we define a partition of d(2) as a family (Volee of subsets Vg S o(1) such that
(xx) = U Va
Vonv = a k.
(4)
q=1
The partition (4) induces a local subsampling of the training set and a local hypothesis space. Namely, define
[mlg {i E [m] d(x,) E VgJ. Hy span VoJ.
Also, we denote by mg # [m]] the local subsampling rate.
Voronoi partitions. Notice that far V, is an arbitrary subset of H. and therefore computing the set [r], can be arbitrarily difficult (e.g., V,, could be defined using an infinite number of constraints and be non-computablle. For this reason, although our theoretical analysis holds for any partition defined as in (4), we focus on the special case of Voronoi partitions, where the subsets (also called cells) are induced by set of Q centroids {Ø(cg)) with Cg X points in the input space. Then, each cell Vg uniquely defined as
Vq {Ø(z) arg min |lő(x) (ee))
ke(q)
with ties broken arbitrarily (c.g.. by assigning the point to the cell with the smaller q). It is now possible to identify the set indices [ni], using the RKHS distance
[|<(()) <(x') K(x, 2) + (('','') 2K(x,
(5)
computing the distance each centroid and taking the minimum.
We remark that our approach based on directly partitioning the feature space has quite different implications compared previous approaches that partition the input space. For example, Voronoi partition of the feature space is very different from a Voronoi partition of the imput space, since the pre-image (x: X' Φ(x) Vg) does not need to follow any Voronoi shape. Moving from input to feature space partitions also opens new computational challenges. For example, we choose to explicitly represent the cell centroid d(c,) so that computing the distance and the assignment of each point to a centroid is a (11) operation. If instead we chose a more complex centroid, such as a cluster barycenter generated by kernel k-means, or an eigenvector computed by kernel PCA, this complexity might be much larger. As an example, the barycenter of a cluster of m points in H. might not correspond any single point in X, and therefore cannot be explicitely represented, but only implicitly as an average of m points in H. Therefore, computing distance to such a centroid O(1). These and more subtle pitfalls appear only when we
would be an O(m) operation rather than a consider the more flexible framework of feature space partitions.
Minimal principal angle. Focusing on partitions of Voronoi type, constructing a good partition is equivalent to choosing set of centroids that preserve the learning accuracy as much as possible. As we rigorously show in Section 4, this can be guaranteed by choosing centroids that maximize the minimal principal angle between subspaces. This quantity frequently appears in the analysis of
3

---

subspace clustering [12], and will be important for us to control both the bias and the variance of our estimator. The first principal angle between two linear subspaces and W of Hilbert space of inner product ( and norm is defined as
LU, W) min mirnaoccos(( (u, w)) u U, w W, ||||| = || = 13.
We call the minimal first principal angle between the subspaces Hg. that is,
min (Hzg H.x) qzt
(6)
Once again, for computational reasons cannot use direct optimization of this quantity in H as objective, since the optimal centroid placement might be impossible to express using points from X. Instead, to promote large principal angles and obtain centroids that computationally easy to handle, we consider the following greedy iterative procedure. Let X = {2s á E [m]]. Then (c),
C1 arg max K ((c, c),
Cq+1 arg
max CEX CEX (C
eEX
where
=
K((,),,,,,,,,,,,,
K(c,
is the Schur complement of candidate centroid with respect to the already selected centroids {cI,..., C. and K, E R2x4 is defined by (K,),, be efficiently computed using rank-1 updates. This strategy has been originally proposed by [9] with the goal of maximizing the volume spanned by the points in the feature space, which is achieved when the angle between all points selected large as required by our condition. Crucially, it is also easy to apply to RKHS's, since computing Schur complements involves only inner products. Beyond promoting large volume and orthogonality, the Schur complement also has important links with uncertainty estimation and spectral approximation. In particular, SCy (c) is also equivalent to the posterior variance of in a Gaussian process [26], and to the leverage score of e w.r.t. the already selected point in the context of randomized linear algebra [21].
Note that the inversion of K2 can
3.2 Learning local KRR estimators by sketched preconditioned conjugate gradient
For each cell of partition, local estimator Fg can be defined as the solution to the local KRR problem
min FEH, Thy kE(rig
Σ 1f(x;) - val + Ag 11fl
(8)
with A Given the local estimators Fap we then define a global estimator f by
= fy if ((x) E Vg-
(9)
Note that the evaluation of the global estimator at point needs only one local estimator.
Guidance on how to pick the values Ag in (8) will follow from our theoretical analysis. Meanwhile, focus on how to efficiently solve the minimization problems (8). Let X, {x, E X : [mly] Ro,xd and Y, {yr Y : (mlg) E R be the local subsets of input/output points, and let Kn, E Xn, be the local kernel matrix with entries (Kn,)i. = K(x;, for i, E [m]y-
Following the same ideas to derive (3), one could compute fa by
[x(x)=
Σ (ag), VEimly
(4qq (Kn, + T) "
(10)
This would already result in a smaller computational burden compared to the vanilla KRR estimator (3): the space and time complexities are now ((maxye(Q) n,,2) and O(EnE(a) 3) potentially with
My R. Moreover, additional saving in time can be obtained by distributing each task (10) over Q different machines, leading to (max, (Q) n2y time complexity. However, the scaling in n, is still quadratic and cubical. To improve these dependencies, we solve (10) only approximately, using the
FALKON algorithm proposed in [29]. To this end, we first need to introduce several key ingredients. While the following constructions hold in general for any set of points, here we adapt them to the partition setting outlined in the previous section.

---

Local Nyström subsampling, For each [Q], we consider a subset of my S ny points
X, fTg.1.... Tgmm, Xg sampled uniformly random from Xg- We then define Km, Rm by (Km,)i kis K (Tgg, (9g) for i, E [mg]] and Kn,m, E 1R Xma by (Kn,m,)is = K (Zi, for [m]g:]] [mg]]
(11)
Local Preconditioner. For each g E [Q]. we define the local (sketched) preconditioner Bg E Rmyxma as
Bg =
Bg Kmm + Xmn
mg
Conjugate gradient descent. We let Bg E Rm be the t-th iteration of conjugate gradient minimiz- ing
1
(3) = - -I ( B X33 (B, Km, Bg)
rgll2 -
(12)
T4
Finally, we define the local FALKON estimator
(z) = > By (
q [Q].
(13)
3.3 ParK
We are now ready to present ParK. Let (p(c)))g=1 with ey X be the centroids of the cells (V₂) selected greedily according to (7), We define the ParK estimator as
=x) = Fz) if of)()) Vgg
(14)
The algorithm to train the above estimator (see Algorithm 1) consists of three main parts. The first such
one greedily identifies the representative points (cy) such that (o(c))) are the centroids of the
cells; the second one identifies the subsets of points Xy. associated to each cell: the third one uses the FALKON algorithm to solve the local minimization problem (12) for each X,, Y, with g E [Q]. thus deriving the Q local estimators (13). prediction time, the algorithm first identifies which cell the test point belongs, and then proceeds using the local estimator of the selected cell to predict output (see Algorithm 2). Note that the RKHS distances in llne of Algorithm and line I of Algorithm 2) are computed using the polarization identity (5).
Algorithm ParK: Train
Require: Training set = (a:) E Raxd y (yi) E IR" numbers of local Nyström centers {mg]g=1 E Nº, local regularization parameters ( E R_, number of local iterations (ta)==1 (t ENQ
1
1: Initialize [n], = 0 for all q
[Q]
2: Greedily select (o(c,)) according to (7)
3: for f 1,
do
4: Compute ar argmin Φ [(((z,) o(c,)11
5: Update [n], [n], u (i)
6: end for
7: for 1,..., do
8: Select X, = (z, X ( and Y, {v E Y i [m]p]
9: Compute Fq.,, as in eq. (13) using X,, Y,
end for
11: Collect the local estimators Jan, and return the ParK estimator 7, as in eq. (14)
5

---

Algorithm 2 ParK: Predict
Require: Test point RR, local estimators ( " representatives of partition C' = (cg)g=1 (l((x) [l((x)
1: Select 7 argminge(q)
((c,)ll
2:
Evaluate (z,)(x)
The time complexity of training Park is log"n) to compute the centroids, ((Qn) to com- pute the indices [n]]g. and (gg) to compute each local estimator. Putting these quantities together we get OQQ'n n log(n) + Enero fg))))g)) in time, and (mmax m2²) in space. If we parallelize the training of the local estimators over Q machines, the time complexity further reduces to ( log(n) m QQ] ng)) In many practical scenarios, we can think Q as (11) For example, in all our experiments we take Q 32 (see Section 5). We compare the complexity of several KRR solver Table 1.
Table I: Computational complexity of some KRR solvers (up constants). For D&C and ParK, report the time complexity on Q parallel machines and the space requirement for each machine.
naive
iterative [37]
Nyström/RI [36, 25] FALKON [29]
D&C [38]
Park
space
n²
n²
m²
m²
(n/()²
max, my 2
time n
fn
m2n
tmn
(n/Q)$
Qonlog(n) + max, fy m, + max, my
test
Space partitioning vs data splitting
We conclude this section commenting on a different yet related distributed approach. As briefly recalled in the introduction, straightforward way to decompose the KRR problem is by simple split of the training data. For example, one can divide the samples uniformly at random into Q disjoint subsets of cardinality n, n/Q. Methods performing such step are known divide-and-conquer [38, 20, 14J. Consisting essentially in a block diagonal approximation of the kernel matrix, the resulting final estimator is an average of globally subsampled models. Divide-and-conquer methods are appealing due to the extreme simplicity the splitting procedure and the direct control of the subsampling rates ng. However, they suffer from worse approximation error (see discussion in [34]), and be expensive at test and evaluation time. On the other hand, partitions present several potential benefits. First, data splitting is byproduct of a geometric partition. This opens to the opportunity of exploiting the structure of the space, for instance enforcing notions of locality or orthogonality. Consistently, the final estimator is union of local estimators, opposed to average of global ones. Hence, partitioning may enhance the approximation power of the model, capturing relevant local correlations [23]. As another consequence, at evaluation time only one local estimator, instead of the average of all estimators, needs to be called, yielding further computational saving. These nice properties have motivated fruitful line of research, notably [23, 34, 24], where the advantage in the partitioning approach has been studied both in statistical and in computational terms. In this paper we concentrate on the computational aspects, expanding on theoretical tradeoffs outlined in [34, 24] and developing [24] with new algorithmic ideas.
4 Theory
To simplify the analysis and better highlight the new ideas in play, we consider the problem (1) in fixed design setting [2, 15], where the 24 are deterministic and the E; are independent and identically distributed random variables.
Let 1.2 Tp)(p) with p Em &..-- We may identify 1,2 with RD endowed with the inner product
(u, w) Lي = fu T w. We define the excess risk of an estimate π of f. in problem (1) as
1 Σ [xx - (x T
R() = 115 -.
(15)

---

We are interested in studying the performance of the estimator J, defined in (14) given partition (4), as measured by the excess risk (15). Our theory will suggest how to construct the partition and tune the regularization in order to get the best learning rate.
4.1 Definitions and assumptions
We start by defining some relevant operators global and local variants. view of (3) (and the fixed design setting), we assume without loss of generality that H = Hn. We define the covariance operator T H - as 1 Eie(m) d(x,) & ((x)), where, for e, w E H, v @ w denotes the operator u E H + (u, H) E H. The operator T is standard in the analysis of kernel methods [7]. We now define the local version of the covariance operator conditioned on the partitioning (4). Thanks to (10) (and the fixed design setting), we can assume without loss of [mlg]- The local covariance operator Tç H. + H. is defined as T. 1 LieM., ((x;)) @ ((x)). We denote with P. H. + H. the orthogonal projection onto the Ri subspace H.g For all E [Q]. we let P. = nn/n. Recall that we denote by the minimal principal angle between the subspaces Hg. as defined in (6).
generality that Hg = (244)
To mcasure the capacity of the hypothesis spaces, we will use the standard notion of effective dimension [7].
Effective dimension. The (global) effective dimension of the space H is given by ((A Tr((T + IT)
A 0.
Consistently, we define the local effective dimension of each space Hg
N₂,(,,) = (T((T, + A,,) >,,) -T,)
Ag > 0.
We also define the local maximal degrees of freedom [2] as
o sup (ó(x), (T, + (x=)) =
Ag ~000
which gives the bound (An) S Noc.a (2e)) S sup, K xx), The effective dimension is related to the spectrum decay of the covariance operator, and thus it provides a way to quantify how many important eigenfunctions the RKHS contains. In this sense, serves as an implicit number of parameters for the nonparametric model represented by the RKHS. The interplay between global and local effective dimensions, hence between global and local model complexity, will play major role in our analysis. Similarly, there exist an interplay between local and a global version of the maximal degrees of freedom (A), which also connected to the coherence of the T operator, and to the concept of maximal leverage score [2].
We will need few basic assumptions.
Assumption 1. f. E H.
Assumption 2. x² = supzex K(x, x) < co.
Assumption 3. The noise variables E; are i.i.d. sub- nub-Gaussian of variance proxy 0² 00, i [m].
Assumption 3 is standard in the analysis of any regression model. In particular, uub-Gaussianity allows control the tails of the noise, and therefore to establish bounds in high probability. Bounded and Gaussian noise are examples, but any variable with sub- Sab-Gassiin tail is covered. Assumptions 1 and 2 are instead typical of kernel methods. With Assumption 1, we suppose that the RKHS is a well specified model. We stick to Assumption 1 for simplicity, but we could easily relax it assuming the existence of function in the RKHS with same excess risk as fa, or considering the excess risk with respect to the best in class. Assumption 2 allows to provide explicit bound for kernel related quantities, and ensures in particular that functions the RKHS are bounded.
4.2
Main results
Our first proposition generalizes the classical bias-variance tradeoff of KRR estimators incorporating iterative optimization, random projections and feature partitioning. The result a high probability bound for the excess risk of our ParK estimator.
7

---

Proposition 1. Let 8 E (0,1). Under the regression model (1) and the assumptions of Section 4.1,
let Jt, be the Park estimator as defined in (14). If for each q E [Q), 0 < Ag S kn?
SN²
40 (1)
my 2 5[1 + 14N (Aa)) log(
2,
b 2 2log
then, with probability at least 1 46,
Q Ng + vNg(A,) log(
1/8) 1ogg(1/))
R(T,)) S 16
Σ q=1
PPaJalla Ap4 + 0²
The proof of Proposition is given in Appendix A.2. The bound consists of bias and variance term. The bias term is an average of local biases, measured by the projection the target function onto the local hypothesis spaces, regularized by local penalization. The variance term is essentially the ratio between the sum of local effective dimensions and the global sample size. We are going to control bias and variance in the next two propositions, whose proof is postponed Appendix A.3. For the bias, we prove the following generalized Bessel inequality. Proposition 2. With the definitions of Section 4.1, we have
Pat. $ Q2 cos co)(l)
Proposition bounds the possible redundancy of the local projections by the minimal principal angle between the local subspaces. In particular, if the local subspaces are an orthogonal decomposition of the global space, the partitioned estimator has no additional local bias. On the other hand, lack of orthogonality results in a larger bias. Turning to the variance, we obtain the following bound on the local effective dimensions.
Proposition 3. With the definitions of Section 4.1. 2g = =Ap= we have
Σ NG(W (1+ n2 9-1
A (0) N(A).
Once again, the minimal principal angle controls the ratio between local and global quantities. Where there perfect orthogonality, splitting the hypothesis space does not increase the effective dimension; otherwise, a price proportional to the minimal principal angle is paid. With the above results in hand, we can now control the excess risk the ParK estimator in terms of the global norm of the target function and the global effective dimension. This allows to compare the performance of our partitioned method to that of a typical global method.
Theorem 4. Let E (0, 1). Under the same assumptions of Proposition I, let q = argmin |E|Q] Pq. for 0 < X S paR2 when Xq Xp= for each E [Q], with probability at least - 46,
402
(1+12000=(0)
R(F.) & 16(1 + Q² cos cos( 1
((A)
If we consider a model where H. is the orthogonal sum of the subspaces Hep as [24], then (os = O, and the bound of Theorem 4 simplifies to R(T.) = OAI]s,111 (1))m-1) In particular, setting À ((1/v1)) optimal rate, in the minimax sense, for global KRR models [7]. Note that, in the orthogonal case, the constraint S Py translates to the minimal local point requirement n, Z vn for all g, and hence to bound on the partition size, namely Q S Vn. On the other hand, when the subspaces H. are perfectly orthogonal, our bound manifests Maaiisical-computationn tradeoff, which is however quantified by the minimal principal angle. Further, the constraints on the local number of Nyström centers my. iterations tg of Proposition 1, and the choice Xg = Ap with À = 1/vn to achieve the minimax rate, allows recover time complexity of O (Qen log(n) + Eve(Q) may vn
obtain the learning rate ((1//n) This is known to be the
log(
nl vn

---

Analyses of (input space) partitioned kernel estimators have been conducted within different models, such as Gaussian SVM's on Voronoi partitions [23]. general kernels on clusters [34], and block- diagonal kernels on arbitrary partitions [24]. In these works, the bounds are established in random design, for plain [23, 34] or Nyström [24] local KRR estimators. Our result in fixed design, but compared to [24] incorporates the additional algorithmic ingredient of iterative optimization. For perfectly orthogonal model (cos = 0), we recover the result in [24] as a special case (although in fixed design). In [34], the bias is controlled choosing same X, for all g. while the crucial bound of Proposition 3 is made as an assumption. Note however that, at least in our proof of Proposition 3, is important to choose differently scaled X for each cell. Furthermore, our analysis and numerical tests motivate that partitioning the feature space is key to control both local bias and local effective dimension. Rather than computational aspects, [23] focuses extending statistical optimality for functions of local smoothness. This theme is also explored in [34, 24]. However, since the proposed partitioning step is either unsupervised [23, 34] unspecified [24], improved rates be obtained only under oracle assumptions, that is, assuming that the smoothness of the target function is localized right on the cells of the chosen partition. Partitions adapting to the unknown local smoothness of the target function can arguably be learned only in a supervised manner. This has been done for piecewise polynomial regression drawing ideas from multi-resolution analysis [4, 18]. An application of these ideas for kernel methods is not straightforward due to the usual computational constraints, but could be subject of future work.
Experiments
In this section we study the performance of ParK on some large-scale datasets (n z 10 107, 10%). In particular we consider dataset where at the moment, because of their cardinality, only a few solvers can efficiently learn from. For this reason we compare to the global large-scale kernel method FALKON which has so far being the method that performs the best in terms of time and accuracy on these datasets [22]. A standard divide-and aeeee-nnd-couuee method can not run on these datasets (for the high space complexity), for this reason we compare with a version where each local estimator is sketched KRR estimator computed with FALKON. We run two different versions of this algorithm, D&C-FALK(vI) and D&C-FALK(v2), that differ only in their hyper-parameters choices as specified later in this section. We also consider a second version of ParK where the centroids of the partition's 's cells are chosen as {o(c,)] with Cg selected uniformly at randomly from the training data X (referred to ParK-Uni). For each experiment we report mean and standard deviation on 10 trials. The experiments are implemented in python using pytorch and the FALKON library [22]. The experiments run on machine with 2 Intel Xeon Silver 41 16 CPUs and 1 GPU NVIDIA Titan Xp. The ram of the machine is 256 GB. We perform experiments on the four large-scale datasets TAXI (n 10", d = 9, regression), HIGGS (n 10 d = 28, classification), AIRLINE (n 2 10 d = 8, regression), AIRLINE-CLS (n RV 105, d = 8, classification) with the same pre-processing and same random train/test split used in [22]. We do not cross validate hyper- prper-parameters of the local estimators of ParK. Instead we use the same used by FALKON in the paper [22] with the following exeptions: let be the global regularization parameters of FALKON and m the number of the Nyström points, the local estimators of ParK use regularization Ag Ap and my = mp, as suggested by the theory. D&C-FALK(v1) also follows the same rule for setting the hyper-parameters of its local estimators, while D&C-FALK centers which are 3m, in AIRLINE and AIRLINE-CLS, 5mg in HIGGS, and 6m in TAXI. The number of centroids used by ParK and D&C-FALK is Q 32 for all experiments. Performance for different Q values remains almost identical but worsen in time for higher values. Further, note that the local estimators of ParK and D&C-FALK are learned sequentially. We report in Table 2 the errors and times. In particular, for ParK(-Umi we report the initialization time that include the greedy algorithm to select the centroids (not for ParK- Uni) and the assignment of the training points to the corresponding cell, the sequential training times of the local estimators, and the total time of this pipeline.
uses the same of the (v1) version except the number of Nyström
We can see that ParK can match the accuracy of the global FALKON estimator with a smaller computational cost. ParK-Uni requires further less time, at the expense of some loss in accuracy, confirming that a worse partition can affect generalization, as suggested by our theory. The reason of the ParK-Uni speedup is twofold. First, the initialization step requires only to assign points to set of randomly selected centroids, and second, the local subsets of points have uniform cardinality (which is usually not the case for normal ParK). D&C-FALK(vI) is the algorithm with with smallest

---

Table 2: Accuracy and running time comparison on large-scale datasets.
TAXI
10"
HIGGS 10"
(MIN.)
ERROR (1-AUC)
(SEC.)
(RMSE)
INIT
TRAIN
TOTAL
INIT
TRAIN
TOTAL
PARK
312.0++0.2
25:11
39=13
64+13
0.182++0.001
30+2
474+172
504+172
PARK-UNI
315.7=00.6
5#1
13=1
18:1
0.192=00.000
3+1
67#7
70+7
FALKON
311.710.1
120.11
0.180.10.001
7151:
D&C-FALK(VI)
356.2100.
1411
0.21211.000
50:11
D&C-FALK(V2)
327.4+0.1
29:11
0.195=10.000
288+2
AIRLINE n 10
AIRLINE-CLS n 10
ERROR (MSE)
TIME (SEC.)
ERROR (C-BRR)
TIME (SBC.)
INIT
TRAIN
TOTAL
INIT
TRAIN
TOTAL
PARK
0.760:10.005
6111
71-99
77::10
31.5+0.2%
9+1
5516
6416
PARK-UNI
0.766=00.06
1+1
32=3
3343
31.6+022%
2+1
22422
24:2
FALKON
0.758++0.005
334+-2
31.5+002%%
391+5
D&C-FALK(Y 0.804
27+1
33.2+ %
20+1
D&C-FALK(V2)
0.799110.05
96111
32.210.1
7311
training time but achieve significantly worse performance using the same rule to choose the number of Nyström points of ParK. For this reason, in D&C-FALK(v2) we increase the number of centroid improve the performance, but the error of the method still results higher than the others with training time now higher than ParK.
6 Conclusions and limitations.
In this paper we have proposed a new algorithm for large scale kernel ridge regression. Our method integrates and jointly exploits three previously uncombined algorithmic strategies, namely partitions, sketching and (preconditioned) iterative optimization. Distinctively from traditional partitioned methods, have introduced the idea of partitioning the feature space, which allows to directly control and resolve the localization of the kernel model. We have presented a simple analysis that characterizes the statistical-comme conppuatocomtttinnnall trade-off of partitioned kernel estimator by the interplay of intuitive quantities. Moreover, we have demonstrated that our algorithm performs favourably against a sfate-ofthhee large scale global method.
The main theoretical limitation our work the lack of a result connecting the proposed partitioning algorithm the properties of the resulting partition. This seems to be a common gap the literature of partitioned kernel methods, where partitions are often assumed to be given or, if explicitly constructed, statistically characterized. While construction theoretically motivated by the analysis and practically validated by the experiments, an actual guarantee is missing. In particular, one could try to prove that the proposed greedy procedure would actually find a maximally orthogonal decomposition of the hypothesis space, under suitable assumptions. From an algorithmic point of view, we point out that the computational cost of the greedy algorithm limits the choice of the partition size. Indeed, large partitions accelerate the training step, but increase the initialization time. We remark, however, that our model is flexible enough include cheaper partitioning options. For example, our experiments show that uniformly chosen partitions can still produce good results.
A Acknowledgments and Disclosure of Funding
The authors thank Nicole Mücke for her useful feedback. This material is based upon work supported by the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF- ICF-1231216 and the Italian Institute of Technology. We gratefully acknowledge the support of NVIDIA Cor- poration for the donation of the Titan Xp GPUs and the Tesla k40 GPU used for this research. L. R. acknowledges the financial support of the European Research Council (grant SLING 819789), and AAA-FRL---OOOR-200
the AFOSR projects FA9550-18-15 18-1-7009, (Furopean Office of Aerospace Research and Development), and the EU H2020-MSCA-RISS project NoMADS DLV-777826.
10

---

References
and [1] H. Avron, V. Sindhwani, and D. Woodruff. "Sketching structured matrices for faster nonlinear
regression". In: Advances in neural information processing systems. 2013, pp. 2994-3002.
[2] F. Bach. "Sharp analysis of low-rank kernel matrix approximations". In: Proceedings of the 26th Conference on Learning Theory (COLT) PMLR 30 (2013), pp. 185-209.
[3] F. Bauer, S. Pereverzev, and L. Rosasco. "On regularization algorithms in learning theory". In:
Journal of Complexity 23.1 (2007), pp. 52-72.
[4] P. Binev, A. Cohen, W. Dahmen, and R. A. DeVore. "Universal Algorithms for Learning Theory Part II: Piecewise Polynomial Functions". In: Constructive Approximation 26.2 (2007),
pp. 127-152.
[5] I. Blaschzyk and I. Steinwart. Improved Classification Rates for Localized SVMs. arXiv:1905.01555 2019.
[6] R. Camoriano, T. Angles, A. Rudi, and L. Rosasco. "Nytro: When subsampling meets early stopping".. In: Artificial Intelligence and Statistics. 2016, pp. 1403-1411.
[7] A. Caponnetto and E. De Vito. "Optimal rates for the regularized least-squares algorithm". In: Foundations of Computational Mathematics 7.3 (2007), pp. 331-368.
[8] L. Carratino, A. Rudi, and L. Rosasco. "Learning with SGD and Random Features". In: Ad- vances Neural Information Processing Systems. Ed. by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett. Vol. 31. Curran Associates, Inc., 2018.
[9] L. Chen, G. Zhang, and E. Zhou. "Fast Greedy MAP Inference for Determinantal Point Process to Improve Recommendation Diversity". In: Advances in Neural Information Processing Systems. Ed. by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. R. Garnett. Vol. 31. Curran Associates, Inc., 2018.
[10] A. Dieuleveut and F. Bach. "Nonparametric stochastic approximation with large step-sizes".
In: The Annals of Statistics 44.4 (2016), pp. 1363-1399.
[11] F. Dumpert and A. Christmann. "Universal consistency and robustness of localized support vector machines" In: Neurocomputing 315 (2018), pp. 96-106.
[12] E. Elhamifar and R. Vidal. "Sparse Subspace Clustering: Algorithm, Theory, and Applications".
In: IEEE Transactions on Pattern Analysis and Machine Intelligence 35.11 (2013), pp. 2765- 2781.
[13] Gonen, Orabona, and S. Shalev-Shwartz. "Solving ridge regression using sketched preconditioned svrg". In: International Conference on Machine Learning. 2016. pp. 1397- 1405.
[14] Z.-. Guo, S.-B. Lin, and D.-X. Zhou. "Learning theory of distributed spectral algorithms", In: Inverse Problems 33.7 (2017), p. 074009.
[15] Hsu, M. Kakade, and T. Zhang. "Random Design Analysis of Ridge Regression", In: Proceedings of the 25th Annual Conference on Learning Theory (COLT) PMLR 23 (2012),
pp. 9.1.4.24
[16] D. Hsu, S. Kakade, T. Zhang, et al. "A tail inequality for quadratic forms of subgaussian random vectors". In: Electronic Communications in Probability 17 (2012).
[17] T. Klock, A. Lanteri, and Vigogna. "Estimating multi-index models with response- conditional least squares". In: Electronic Journal of Statistics 15.1 (2021), pp. 589-629,
[18] W. Liao, M. Maggioni, and S. Vigogna. "Multiscale regression unknown manifolds". In: Mathematics Engineering 4.4 (2022), pp. 1-25.
[19] J. Lin and V. Cevher. "Kernel Conjugate Gradient Methods with Random Projections". In: Applied and Computational Harmonic Analysis (2021).
[20] S.BB. Lin, X. Guo, and D.-X. Zhou. "Distributed Learning with Regularized Least Squares".
In: Journal of Machine Learning Research 18.92 (2017), pp. 1-31.
[21] M. W. Mahoney and P. Drineas. "CUR matrix decompositions for improved data analysis". In:
vol. 106. 3. 2009, pp. 697-702.
[22] G. Meanti, L. Carratino, L. Rosasco, and A. Rudi. "Kernel Methods Through the Roof: Handling Billions of Points Efficiently". In: Advances Neural Information Processing Systems. Ed. by H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin. Vol. 33. Curran Associates, Inc., 2020, pp. 14422-14422.

---

[23] M. Meister and I. Steinwart. "Optimal Learning Rates for Localized SVMs". In: Journal of Machine Learning Research 17.194 (2016), pp. 1-44.
[24] Mücke. "Reducing training time by efficient localized kernel regression". In: Proceedings
of the Twenty-Second International Conference Artificial Intelligence and Statisties. Vol. 89.
2019, pp. 2603-2610.
[25] A. Rahimi and B. Recht. "Random features for large-scale kernel machines". In: Advances in neural information processing systems. 2008, pp. 1177-1184.
[26] C. E. Rasmussen and C. K. I. Williams. Gaussian processes for muchine learning. Adaptive computation and machine learning. MIT Press, 2006. ISBN: 026218253X.
[27] D. Richards and P. Rebeschini. "Optimal Statistical Rates for Decentralised Non-Parametric Regression with Linear Speed-Up". In: Advances in Neural Information Processing Systems.
Vol. 32. 2019.
[28] A. Rudi, R. Camoriano, and L. Rosasco. "Less is more: Nyström computational regularization".
In: Advances in Neural Information Processing Systems. 2015, pp. 1657-1665.
[29] A. Rudi, L. Carratino, and L. Rosasco. "FALKON: An optimal large scale kernel method". In: Advances in Neural Information Processing Systems. 2017, pp. 3891-3901.
[30] A. Rudi and L. Rosasco. "Generalization Properties of Learning with Random Features". In: Advances in Neural Information Processing Systems 30. 2017, pp. 3215-3225.
[31] Schölkopf and A. Smola. Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT Press, 2002.
[32] I. Steinwart and A. Christmann. Support vector machines. Springer-Vertee New York, 2008.
[33] I. Steinwart, D. R. Hush, C. Scovel, et al. "Optimal Rates for Regularized Least Squares Regression." " In: COLT, 2009,
[34] R. Tandon, S. Si, P. Ravikumar, and I. Dhillon. Kernel Ridge Regression via Partitioning. arXiv: LeXiv:1608.000 1600.61966.666. 2016.
[35] P. Thomann, I. Blaschzyk, M. Meister, and I. Steinwart. "Spatial Decompositions for Large Scale SVMs". In: Proceedings of the 20th International Conference Artificial Intelligence and Statistics. Vol. 54. Proceedings of Machine Learning Research. PMLR, 2017, pp. 1329- 1337.
[36] C. K. Williams and M. Seeger. "Using the Nyström method to speed up kernel machines". In: Advances in neural information processing systems. 2001, pp. 682-688.
[37] Y. Yao, L. Rosasco, and A. Caponnetto. "On Early Stopping in Gradient Descent Learning". In: Constructive Approximation 26.2 (2007), pp. 289-315.
[38] Y. Zhang, J. Duchi, and M. Wainwright. "Divide and Conquer Kernel Ridge Regression:
A Distributed Algorithm with Minimax Optimal Rates". In: Journal of Machine Learning Research 16.10 6.1002 (2015), pp. 2390-
12


Reviews:
[{"review_id": "oGAyxqXRAJf", "reviewer": "Reviewer_J64c", "summary": "The paper proposes a computationally efficient algorithm for KRR, called \"ParK\". ParK combines three ideas in the literature: partitioning input samples, the FALKON solver for KRR, and Nystrom subsampling. After partitioning, FALKON and Nystrom are applied black-box based on the prior literature, while the partitioning is claimed to be new (I am not an expert on partitioning techniques). The partitioning is a Voronoi Diagram in the feature space whose centers are heuristically chosen to promote orthogonality between local estimators.\n\nA few experiments are run, with emphasis on improved runtime. Statistical guarantees are in the fixed design setting, and most of the theoretical effort is spent showing that if the partitions are nearly orthogonal in the feature space, then ParK has near-optimal statistical complexity.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper is both nice and approachable. I recommend it for publication.\n\nIt would benefit from minor revisions throughout, but those are just simple fixes for the sake of clarity.\n\nMy background is largely theoretical in KRR, so I will mostly discuss that aspect, and not focus on the experimental side so much. To me, the experiments seem reasonable and promising.\n\nThe ParK method is a combination of largely existing ideas. Nystrom and FALKON are both pre-existing, but the feature-space partitioning appears to be novel (I am not an expert on partitioning techniques). Since the feature space partitioning is the focus of the paper, the paper seems original enough and distinct from prior works. I feel the idea is significant enough to merit publication.\n\nRelated works are well cited throughout the paper.\n\nThe paper seems to be technically sound, and I tried to review all the proofs in it, but the mathematical notation can be frustrating at times. Further, the proofs lack clarity in several nontrivial steps. There are several reasonable-looking mathematical statements that I did not verify. I list concrete examples at the end of this text, and would strongly recommend clarifying much of this.\n\nAside from these mathematical frustrations, the paper is clear and interesting.\n\nThe theoretical results are sufficient, but not particularly exciting. The fixed design setting is okay, but I am admittedly a little concerned about how well the voronoi diagram will generalize near the edges of the diagram (i.e. on the edges of the partition).\n\nThe main stated limitation of not guaranteeing highly orthogonal subspaces is earnest and interesting. I am curious to see if future works can provably increase the angle between subspaces / decrease the incoherence of each partition's kernel. (see below for more on incoherence)\n\nSome smaller questions:\n- [Line 52] Is there any reason why the input values $x_1,...,x_n$ need to be real-valued? This seems like a strong constraint which is rarely needed in kernel problems.\n- [Line 178] Since ParK is not an average of local estimators, does it have discontinuities at the boundaries of the partitions? Will this hurt the out-of-sample performance of ParK, or limit it's statistical guarantees for $x\\\\notin\\\\{x_1,...,x_n\\\\}$?\n- [Lines 207, 230] The term $N_\\{q,\\infty\\}$ seems to be an incoherence parameter (specifically, it is the largest ridge leverage score). This means that Theorem 1 and Corollary 4 both have the number of samples in a partition depend on the incoherence of the kernel matrix of that partition.\n    - Is this circular dependence inevitable?\n    - This interpretation seems interesting enough of a context to be worth mentioning?\n\n\nUnclear points:\n- $\\\\|\\cdot\\\\|$ is never defined without subscript, for neither vectors nor matrices, despite the fact that both appear in the appendix. I understood these to be the L2 norm for vectors, and the spectral norm for matrices.\n- In Algorithms 1 and 2, the process of finding the corresponding centroid should be written in terms of the kernel function, since we generally don't expect algorithms to explicitly build vectors in the feature space.\n- [Line 98] Does finding a centroid not take $O(Q)$ time, not $O(1)$. $Q$ is assumed non-constant in other big-O notations later on.\n- [Line 191] The overall notation of $L^2 = L^2(\\\\rho)$ often confuses me. The inner product $\\\\langle u,w \\\\rangle_{L^2} = \\\\frac1n u^\\\\intercal w$ suggests that $u\\\\in L^2$ is a real vector, presumably with $n$ components. The following $\\\\|\\hat f - f_*\\\\|_{L^2}$ suggests that real-valued functions of real-inputs lie in $L^2$.\n    - The confusion propagates to the Appendix A.1, where the distinction between functions and vectors gets much more confusing. Does $Sf(x)$ denotes $\\[Sf\\](x)$ or $S\\[f(x)\\]$? If $Sf(x)=\\\\langle f, \\phi(x) \\\\rangle_{H}$, this seems to be a real number so how can it lie in $L^2$?\n- [Line 519] The W matrix should have different subscripts in the bottom left and top right corners, in both block matrices.\n- [Line 522] Specify $\\beta_q \\geq 0$, which is used a few lines lower.\n- [Line 525] In the last line of math, what inequality is used to pull this spectral norm of $M_q$ out?\n- [Line 526] The 5th line, I think, uses an adapted version of the Push-Through identity (specialized to use a scaled identity matrix)? I would clarify this.\n- [Line 527] Again, I think an adapted push-through identity was used? Clarify this.\n- [Line 528] What inequality is used in the final line? I get the triangle inequality to separate the $+1$ term, but how did $\\\\lambda_q$ get outside the norm?\n- [Line 528] How does $\\\\frac1n \\\\sum_\\{\\\\phi(x_i)\\\\notin V_q\\} \\\\|P_q\\phi(x_i)\\\\|$ relate to the right hand side? Was this bounded anywhere else in the paper?\n\n\nRecommended edits:\n- State somewhere early on that you take $Q=O(1)$ empirically, around 32, to help the audience understand the scale of the partitioning scheme.\n- [Line 74] The sentence \"Unlike classical ...\" is awkwardly phrased\n- [Line 117] This Schur complement also arises as the variance of a Gaussian Process at $c$ given $c_1,\\\\ldots,c_q$. Seems like another interesting interpretation for the heuristic, in case you want to mention that.\n- [Line 198] The $\\\\otimes$ notation is only used in some KRR communities. I recommend defining it in section 2.\n- [Line 230] It seems odd to describe Theorem 1 as a new theorem. The proof of it just plugs together theorems of prior works, and spends most of its time working with notation. I would at least elevate Corollary 4 to be a theorem, and considering demoting this to be a corrolary.\n- [Line 454] It may be unwise to define both the symbols $L$ and $L^2$, where one is an operator and the other is a space. I'd recommend using different symbols.\n- [Line 527] The inline equation's first step was unclear to me for a while. Simply adding \"Taking this inequality and setting...\" to the start of the sentence would be nice.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors are candid and clear with the limitations of the work. Most centrally, there lacks a guarantee connecting the partitioning algorithm to the orthogonality of the resulting partitions.\n\nWhile it would, of course, be better to have seen such a result in this paper, the empirical evidence for the algorithm is convincing enough. Further, provably promoting such orthogonality is a problem of large enough scope to merit a completely separate publication.\n\n\nSocietal Impact: No concerns; not relevant.", "review_text": "", "extra_scores": {"ethical_concerns": "No concerns; not relevant.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "7"}, "extra_text": {}}, {"review_id": "SSjSYUcdhnV", "reviewer": "Reviewer_o4bz", "summary": "The paper proposes a large-scale solver for kernel ridge regression. The main idea is to partition the input space using clustering in the kernel feature space and then solve local regression problems over the clusters using the FALKON method (Rudi et al., NIPS 2017). During inference, one first loads the partition to which a test point is assigned and then returns the prediction using the hypothesis learned over that partition/cluster. \n\nTheoretical guarantees are provided in the fixed design setting and rely on the previous analysis done by Rudi et al. (2017). The empirical results demonstrate that the approach can attain the accuracy of the FALKON method while making significant computational savings.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I find that the paper is clearly written and well-organized. It is quite easy to follow the presentation and pretty much all of the arguments (apart from lines 75-78: the reference to Lemma 4.55 in [28] is not sufficiently clear).\n\nPartitioning data in the kernel feature space can be computationally expensive. To avoid that the authors propose an iterative procedure for the selection of centroids such that they are maps of instances from the input space. Otherwise, the pre-image of a kernel feature space centroid might not exist and would complicate the inference. The proposed iterative procedure is a greedy volume maximization approach, similar to selection of landmarks in the Nystroem method based on determinantal point processes. In each step of this iterative procedure, one needs to invert a matrix of size q (Eq. 6), where q denotes the cardinality of the current set of selected centroids. While this is an effective strategy one could ask whether an improvement might be possible using an approach with a linear complexity (e.g., kernel k-means++ sampling scheme). \n\nThe hypothesis that is learned over a data partition uses the FALKON method and it would be good to add a discussion to the main part of the paper, clarifying in what way the proof of Theorem 1 is different from that work. In my brief pass over the appendix, I would say that there is no substantial extension of that work and that the main theoretical contributions of this paper are Propositions 2 and 3. The latter two results state that the selection of data partitions such that they are orthogonal is beneficial for the effectiveness of the local approximations. Would it be possible to link results in Propositions 2 and 3 a bit stronger to the selected data partitioning method?\n\nEmpirical study has been performed on four datasets with millions of points. The provided evidence clearly supports the argument for computational effectiveness of the approach. Now, what is lacking here is the dimensionality of the input space and how the proposed approach fares relative to that parameter. Namely, one needs to solve the partitioning problems which are affected by the dimension of the instance space and there should be a threshold where computational savings compared to the FALKON method become void or less relevant. Would it be possible to list the dimensions of these datasets and possibly add an experiment where the dimension increases ($d = 10,100,1000,10000$) -- this could be a synthetic problem purely to showcase the dependence on the input space dimension.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "This has been adequately addressed.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "A-mZJsjZSLr", "reviewer": "Reviewer_iZMk", "summary": "This paper deals with scaling kernel ridge regression, providing a new large-scale algorithm to solve kernel ridge regression called ParK. ParK is based on combining several ideas: partitioning the training set so training can be done on smaller datasets, using a geometric approach to build the partition (i.e. space partitioning as opposed to data splitting), building the split in the feature space, forming the global estimator as a \"stitch\" of local ones, computing the local estimators using preconditioned CG where the preconditioner is formed using Nystrom approximation.  Theoretical analysis in a fixed design setting is provided. Basic experiments are performed.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This is a well written paper that suggests a new method (ParK) for solving large-scale KRR problems. The idea of splitting the global estimator as a stitch of local-in-feature-space estimators is interesting. This idea is complemented with an interesting theoertical analysis.\n\nWhile the method has desirable computational properties, and generalizes well according to the theory (provided we have a good partition), they method and the associated theory has several shortcomings, not all of which are discussed by the authors:\n1)  The estimator that is found is discontinuous since it is a stitch of several local estimators. \n2) Since the partition is geometric, there is guarantee that it will be balanced. In fact, it is possible for almost all points to belong to a single local estimator, thereby ruining any computational gain.\n3) In the theory, the parameters m_q and t_q are set using quantities that are unknown (f_\\star and N_\\infty,q). For t_q, this is less an issue because we can simply check if the the method converged or not. However, for m_q we just need to hope the parameter is set high enough (but not too high so to waste computational resources) \n4) Generalization bounds are based on the minimal angle between the space spanned by the partition. There is a algorithm that partitions the points so to maximize this angle. However, it seems that this algorithm is heuristic, and no formal bound is given.\n\n(all the above points should be discussed in the paper)\n\nNovelty: of all the ideas combined to form ParK (see summary of paper), only the idea of partitioning using the feature space seems novel. All other ideas appeared elsewhere, e.g. the actual algorithm that is used to split the data seems to be from [8]? Preconditioning using Nystrom is from the FALKON paper. Much of the description in the main paper seems to be redundant and just recasting things in term of the local estimator. For example, what is essentially new in Section 3.2? Also: some of the theoertical results seems to based on adapting similar results from [24].\n\nExperiments: somewhat basic and not really convincing. \n\nMinor comments: \n1) Notation is sometimes ambiguous, e.g. H_q seems to denote both space spanned by q vectors and the one spanned by partition q.\n2) Might be useful to mention that N_q,infty is essentially the coherence of T (maximal leverage score). \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No societal impact. Limitations are discussed in main review. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "4YUcOuUN8Z9", "reviewer": "Reviewer_jHTA", "summary": "This paper studies the problem of kernel ridge regression. They propose a large-scale solver named ParK, which combines multiple techniques from the literature, including sketching and preconditioning. Their main conceptual contribution is to use a divide-and-conquer approach that partitions the feature space rather than the input space. The authors show that this makes the algorithm faster without affecting its statistical performance.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper provides a solid contribution by pushing the computational limits of kernel ridge regression solvers. Moreover, I believe the ideas in this paper could help improve the state-of-the-art software. Nonetheless, the proposed method is not particularly original as it primarily combines ideas from [8] and [19].\n\nTo the best of my knowledge, the technical contributions are simple and mathematically sound. The authors present a fair amount of theoretical and practical evidence to support their claims. In particular, they prove nontrivial error bounds, that are minimax optimal for some settings, and show numerically that their implementation is a constant times faster than FALKON [19] (a state-of-the-art solver). Broadly speaking, the paper is well-written and easy to read. However, there are some points I'd like the authors to address to improve this manuscript; see below.\n\n# Comments\n- I'd like to see a better comparison with existing methods. For example, I'd like to know how ParK compares in terms of sample complexity, risk bounds, and computational complexity. Maybe a table?\n- On a related note, the guarantees prove in this work are for a fixed design (i.e., fixed features x_i). However, most of the literature concerning partitions has guarantees assuming a random design. Therefore, it would be interesting to extend the analysis to random settings and see how they compare.\n- I would suggest the authors work a bit on the introduction. I felt the second paragraph is broad and hard to follow without more context. It will be better to start with the contributions of the paper rather than related work\n- The authors should specify the RAM of the machine used to run the experiments.\n- For me is a little bit of an overstatement to say that Corollary 4 shows a  statistical-computational tradeoff. This result proves an upper bound for one method. It doesn't give a lower bound for a meaningful class of polynomial-time methods.\n\n# Questions\n- The authors claim that using partitions in the feature space allows them to promote orthogonality among the estimators. This is intuitively true, but is there a way to formalize it and compare it with other approaches? Why isn't it possible to do something similar with previous estimators?\n- The following comment might be a result of using different machines, but I found it interesting and wondered if the authors have an explanation for it. The experiments used the same datasets as the ones used by FALKON [19]. Yet, the running times are not consistent with the ones obtained in [19]. For example for TAXI, [19] got a run time of 3628 sec (as opposed to 120 sec in the present manuscript), while for AIRLINE-CLS, the time reported by [19] was 186 sec (in comparison, the experiment in this paper gave 391 sec).\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes, they did address these issues. This work doesn't have a significant societal impact as it deals with a theoretical question.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The focus of the submission is the solution of kernel ridge regression (KRR) in the large-scale setting. Particularly, the authors propose to use a combination of partitioning (a Voronoi one in the feature space), preconditioned iterative optimization and sketching. The resulting KRR solver is shown to be able to achieve the minimax optimal global KRR rate in the fixed design setting in terms of the excess risk (14), while significantly reducing the computational load. The efficiency of the approach is demonstrated on 4 large-scale benchmarks.

Kernel methods are among the most successful machine learning tools of our times. By the construction of a novel KRR solver with optimal statistical performance and reduced computational complexity, the authors provide a valuable contribution to the field which is of both theoretical and practical interest.
Decision: Accept (Poster)