Kernel Quadrature with Randomly Pivoted Cholesky
Authors: Ethan Nicholas Epperly
Elvira Moreno Ferreira
Venue: NeurIPS 2023 spotlight (2023)
Keywords: kernel quadrature
Nyström approximation
reproducing kernel Hilbert space
randomly pivoted Cholesky
TL;DR: We develop new kernel quadrature schemes based on the randomly pivoted Cholesky sampling algorithm

Abstract: This paper presents new quadrature rules for functions in a reproducing kernel Hilbert space using nodes drawn by a sampling algorithm known as randomly pivoted Cholesky. The resulting computational procedure compares favorably to previous kernel quadrature methods, which either achieve low accuracy or require solving a computationally challenging sampling problem. Theoretical and numerical results show that randomly pivoted Cholesky is fast and achieves comparable quadrature error rates to more computationally expensive quadrature schemes based on continuous volume sampling, thinning, and recombination. Randomly pivoted Cholesky is easily adapted to complicated geometries with arbitrary kernels, unlocking new potential for kernel quadrature.

Kernel Quadrature with Randomly Pivoted Cholesky
Ethan N. Epperly and Elvira Moreno"
Department of Computing and Mathematical Sciences
California Institute of Technology
Pasadena, CA 911 125
( emoreno2)@caltecch edu
Abstract
This paper presents new quadrature rules for functions in a reproducing kernel Hilbert space using nodes drawn by sampling algorithm known randomly pivoted Cholesky. The resulting computational procedure compares favorably previous kernel quadrature methods, which either achieve low accuracy require solving a computationally challenging sampling problem. Theoretical and numerical results show that randomly pivoted Cholesky fast and achieves comparable quadrature error rates to more computationally expensive quadrature schemes based continuous volume sampling. thinning, and recombination. Randomly pivoted Cholesky is easily adapted to complicated geometries with arbitrary kernels, unlocking new potential for kernel quadrature.
Introduction
Quadrature is one of the fundamental problems in computational mathematics, with applications
Bayesian statistics [35), probabilistic ODE solvers (27). reinforcement learning (32), and model-
based machine learning 30 The task is to approximate an integral of a function ff by the weighted
sum of f's values at judiciously chosen quadrature points $1,..., Sn:
[ f(z)g(x) dulz) s Σ w.f(s:).
(I)
1=1
Here, and throughout, X denotes topological space equipped with Borel measure μ., and E (2)) denotes a square-integrable function. The goal of kernel quadrature is to select quadrature weights
W...... W., E I and nodes s...... 8, X which minimize the error in the approximation (I) for all
drawn from a reproducing kernel Hilbert space (RKHS) H of candidate functions.
The ideal kernel quadrature scheme would satisfy three properties:
1. Spectral accuracy. The eтror of the approximation (I) decreases at rate governed by the eigenvalues of the reprodueing kernel of H. with rapidly decaying eigenvalues guarantee- ing rapidly decaying quadrature error.
2. Efficiency. The nodes . Sn and weights Wi...., Wn can be computed by an algorithm which is efficient in both theory and practice.
ENE acknowledges support by NSF FRG 1952777, Carver Mead New Horizons Fund, and the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0021 100. EM was supported part by Air Force Office of Scientifie Research grant F49550-22-1-0225.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

The first two goals may be more easily achieved if one has access to Mercer decomposition
k( v) Σ (x)(())))))
(2)
where G1, C..... form an orthonormal basis of L²(()) and the eigenvalues X1 2 X2 are decreasing. Fortunately, the Mercer decomposition is known analytically for many RKHS's H. on simple sets X such as boxes [0, 1| However, such a decomposition is hard to compute for general kernels k and domains X. leading the third of our desiderata:
Mercer-free. The quadrature scheme can efficiently implemented without access to an explicit Mercer decomposition 2
Despite significant progress in the probabilistic and kernel quadrature literature (see, e.g., B 18 21 23 26 41), ongoing.
the search for kernel quadrature scheme meeting all three criteria remains
Contributions. We present new kernel quadrature method based on the randomly pivoted Cholesky (RPCHOLESKY) sampling algorithm that achieves all three of the goals. Our main contributions
1. Generalizing RPCHOLESKY sampling (introduced [0] for finite kernel matrices) to the continuum setting and demonstrating its effectiveness for kernel quadrature.
2. Establishing theoretical results (see theorem D which show that RPCHOLESKY kernel quadrature achieves near-optimal quadrature error rates.
3. Developing efficient rejection sampling implementations (see algorithms 2 and of RP- CHOLESKY in the continuous setting, allowing kernel quadrature to be applied to general spaces X, measures μl, and kernels & with ease.
The remainder of this introduction will sketch our proposal for RPCHOLESKY kernel quadrature. comparison with existing kernel quadrature approaches appears in section g
Randomly pivoted Cholesky. Let be an RKHS with kernel that is integrable the diagonal
[ M(z, x)) dy(z)) < +00.
RPCHOLESKY uses the kernel diagonal k(x, x)) as a sampling distribution to pick quadrature nodes. The first node is chosen to be sample from the distribution k(2; dp(x), properly normalized:
k(0, x) du(x))
81
Jx k(x x) du(æ))
Having selected $1, remove its influence the kernel, updating the entire kernel function:
k k(x v) k(x,)) ((2,81)k(111.))
8()1,81)
Linear algebraically, the update can interpreted as Gaussian elimination, eliminating "row" and "column" 81 from the "infinite matrix" k. Probabilistically. if we interpret l: as the covariance function of a Gaussian process, the update 4 represents conditioning on the value of the process at 81. We use the updated kernel & to select the next quadrature node:
1e' (x, d)(x))
82
Ix Wx, x) du(x))
whose influence is then subtracted off k'' as in (4) RPCHOLESKY continues along these lines until nodes have been selected. The resulting algorithm is shown in algorithm T Having chosen the nodes choice of weights wi,..., W,, is standard and is discussed in section 3.2
$1....;
8n, .......,n.....
RPCHOLESKY sampling is more flexible than many kernel quadrature methods, easily adapting to general spaces X, measures μu, and kernels k. To demonstrate this flexibility, we apply RPCHOLESKY to the region fig. eqleppeeed with the Matérn 5/2-kernel with bandwidth 2 and the measure
du(æ, 4) (x² + +²) dæ dy.

---

Algorithm I RPCHOLESKY: unoptimized implementation
Input: Kernel k X x X + I and number of quadrature points v
Output: Quadrature points $1,.... Sn
1: for i 1,2,...... n do
2: Sample s; from the probability measure k(z, 2) dp(x)/ Ix k(, 2) 2( du(a)
3: Update kernel t k(, 8.)(ss 8;)-"((s,,,
4: end for
set of n 20 quadrature nodes produced by RPCHOLESKY sampling (using algorithm 2 shown in fig. Ta The cyan cyan-pink pink shading shows the diagonal of the kernel after updating for the selected nodes. We see that near the selected nodes, the updated kernel is very small, meaning that future steps of the algorithm will avoid choosing nodes in those regions. Nodes far from any currently selected nodes have much larger kernel value, making them more likely to chosen in future RPCHOLESKY iterations.
0.035
Mean relative quadrafure eтτor
0.03
100
0.025
10
0.5
0.02
10
-055
Monte Carlo
0.01
10-3
IID
RPCholesky (Ours)
-1.5
0.005
Ob 4a-4
10
-2
10
10?
-2
-1
0
2
Number of nodes
(a)
(b)
Figure RPCHOLESKY kernel quadrature on oddly shaped region. Left: Black dots show 20 points selected by RPCHOLESKY ues k((z, v), (I, after removing influence of selected points. Right: Mean relative error for RPCHOLESKY kernel quadrature, iid kernel quadrature, and Monte Carlo quadrature for f(z. v) = (x) exp(y) with 100 trials. Shaded regions show 10%/90% quantiles.
crescent-shaped region. Shading shows kernel val-
The quadrature enror for f(r,y) sin(x) exp(v),
and different numbers of quadrature nodes for RPCHOLESKY kernel quadrature, kernel quadrature with nodes drawn iid from // (XX. and Monte Carlo quadrature methods would be difficult to implement in this setting because they would require an explicit Mercer decomposition (projection DPPs and leverage scores) volume sampling). A comparison of RPCHOLESKY with more kernel quadrature methods on a benchmark problem is provided in fig.
shown in fig. [b Other spectrally accurate kernel quadrature
expensive sampling procedure (continuous
Related work
Here, we overview past work on kernel quadrature and discuss the history of RPCHOLESKY sampling.
2.1 Kernel quadrature
The goal of kernel quadrature is to provide systematic means for designing quadrature rules RKHS's. Relevant points of comparison are Monte Carlo [36] and quasi-Monte Carlo 14] methods, which have (lI/Vn) and (ppoyylog(n)/n) convergence rates, respectively.
The literature on probabilistic and kernel approaches to quadrature vast, so any short summary necessarily incomplete. Here, we briefly summarize some of the most prominent kernel quadrature

---

methods, highlighting limitations that we address with our RPCHOLESKY approach. We refer the reader [22] Tab. 1] for helpful comparison of many of the below-discussed methods.
iteratively select quadrature nodes using greedy approach. These methods face two limitations. First, they require the solution of a (typically) nonconvex global optimization problem at every step, which may be computationally costly. Second, these methods exhibit slow (1/nn) quadrature error rates [[T] Prop. 4] (or even slower [16]..
Herding, Kernel herding schemes
26
Optimally weighted herding is known as sequential Bayesian quadrature [26]
Thinning. Thinning methods 17
try to select n good quadrature nodes from larger set of N >> n candidate nodes drawn from simple distribution. While these methods have other desirable theoretical properties, they do not benefit from spectral accuracy.
18 38
Leverage scores and projection DPPs. With optimal weights (section 3.2 quadrature with nodes
sampled using a projection determintal point process (DPP) [7] (see also 15 21 or iid from the (ridge) leverage score distribution [3] achieve spectral accuracy. However, known efficient sampling algorithms require access the Mercer decomposition (2). limiting the applicability these schemes to simple spaces X, measures 14, and kernels k where the decomposition is known analytically.
Continuous volume sampling. Continuous volume sampling is the continuous analog of n-DPP sampling (28], providing quadrature nodes that achieve spectral accuracy [8] Prop. 5]. Unfortunately, continuous volume sampling computationally challenging. In the finite setting, the best-known algorithm
9 for exact n- n-DPP sampling requires a costly O([X| n.5 + nn.5) operations. Inexact samplers based on Markov chain Monte Carlo (MCMC) (e-g.. 0 34 may be more competitive, but
the best-known samplers in the continuous setting still require an expensive O(1 log MCMC steps 34 Thms. 1.3-1.4]. RPCHOLESKY sampling achieves similar theoretical guarantees to continuous volume sampling (see theorem D and can be efficiently and exactly sampled (algorithm
Recombination and convex weights. The paper (see also 23 proposes two ideas for kernel quadrature when is a probability measure and y = 1. First, they suggest using a recombination algorithm (e.g.. (29) to subselect good quadratures from N >> n candidate nodes iid sampled from μ.. All of the variants of their method either fail to achieve spectral accuracy require an explicit Mercer decomposition 22 Tab. 1]. Second, they propose choosing weights (W1,.... ww,) that are convex combination coefficients. This choice makes the quadrature scheme robust against misspecification the RKHS 전 f. among other benefits. It may be worth investigating combination of RPCIIOLESKY quadrature nodes with convex weights in future work.
2.2 Randomly pivoted Cholesky
RPCHOLESKY was proposed, implemented, and analyzed in [00] for the task of approximating an M x M positive-semidefiiie matrix A. The algorithm is algebraically equivalent to applying an earlier algorithm known as adaptive sampling [12 13 to 1/2 Despite their similarities, RPC- HOLESKY and adaptive sampling are different algorithms: To produce rank-n approximation to an M x M matrix, RPCHOLESKY requires On MM) operations, while adaptive sampling requires much larger ( operations. See [10 $4.1] for further discussion on RPCHOLESKY'S history. In this paper, introduce continuous extension of RPCHOLESKY and analyze its effectivenees for kernel quadrature.
Theoretical results
In this section, prove main theoretical result for RPCHOLESKY kernel quadrature. We first establish the mathematical setting (section 3.1) and introduce kernel quadrature (section 3.2. Then, present main theorem (section 3.3) and discuss its proof (section 3.4
3.1 Mathematical setting and notation
Let μ be a Borel measure supported on a topological space X and let H be RKHS on with continuous kemel & X X R. We assume that + k(x 2) is integrable (3) and that H

---

is dense in (2). These assumptions imply that H. is compactly embedded in (2(,,, the Mercer decomposition 2 converges pointwise, and the Mercer eigenfunctions form an orthonormal basis of (2())) and an orthogonal basis of H (39 Thm. 3.1].
Define the integral operator
Tf(x) |. k(x, uvf(w) du(y).
(5)
Viewed operator T ( + (M), T is self-adjoint, positive semidefinite, and trace-class.
One final piece of notation: For function 8 X and an ordered (multi)set {si..... .... sn we let 5(S) be the column vector with ith entry ((s)). Similarly, for bivariate function h X X + R, h(,,5 S) (resp. h(s,.) -)) denotes the row (resp. column) vector-valued function with ith entry h( 8:) (resp. h(s;,)) (), and S) denotes the matrix with ij entry h(;;,85)
3.2 Kernel quadrature
Following earlier works (e.g., 3 7), let us describe the kernel quadrature problem more precisely.
Given function g (H), we seek quadrature weights w (wi..... Wnn E R and nodes
(s1...., 3n C that minimize the maximum quadrature error over all f H with I[ S
Err(S, w: = max
(())(()))x(z) - Σ wi [
(6)
II S1
short derivation
3 Eq. (7)) or appendix
Cp yields the simple formula
Err(S, w;gg)
Tg Σ w;k(- 8:)
(7)
This equation reveals that the quadrature rule minimizing Err(S, w;g) is the least-squares approxi- mation of Tg E H by a linear combination of kernel function evaluations E wi,k(,,8))
If we fix the quadrature nodes {s1,... ................................ sn), ) the optimal weights u. = (W+1,.... W.n)
RR
minimizing (rr(S, w: g) are the solution of the linear system of equations
kS,,Sww. Tg(S).
(8)
We (8) to select the weights throughout this paper, and denote the error with optimal weights by
Err(S:g) = Err(S, w.;g).
3.3 Error bounds for randomly pivoted Cholesky kernel quadrature
Our main result for RPCHOLESKY kernel quadrature is as follows:
Theorem Let S (s1..... sne be generated by the RPCHOLESKY sampling algorithm. For any function (p),,, nonnegative integer 2 0, and real number 5 > 0, we have
E Er: s |1 (olll
N
1=++1
provided
2\1
n 2 r rlog
+
(10)
3o Disr+1
To illustrate this result, consider an RKHS with eigenvalues X O(i-2x) for > 1/2. An example is the periodic Sobolev space H (00, 1]) (see section 5) By setting i = 1/r, we see that
E (s: - Or [|g (ي)
for RPCHOLESKY with = 2( logr).
The optimal scheme requires n (rr) nodes to achieve this bound g Prop. $2.5). Thus, to achieve an error of -2*), RPCHOLESKY requires just logarithmically more nodes than
5

---

the optimal quadrature scheme for such spaces. This compares favorably to continuous volume sampling, which achieves the slightly better optimal rate but much more difficult sample.
Having established that RPCHOLESKY achieves nearly optimal error rates for interesting RKHS's ss make general comments theorem [ First, observe that the error depends the Eier+1 A; of the tail cigenvalues. This tail sum is characteristic of spectrally accurate kernel quadrature schemes g The more distinctive feature in (I0) the logarithmic factor. Fortunately, for double precision computation, the achieveable accuracy is bounded by the machine precision 2-32 so this logarithmic factor effectively modest constant:
2
log
& log(2"2) < 37.
3.4 Connection Nyström approximation and idea of proof
We briefly outline the proof of theorem 1
Following previous works (e.g.. 3 23 ) we first utilize the
connection between kernel quadrature and the Nyström approximation 31 42
the kernel
Definition 2. For nodes {81,..., 8n) C X, the Nyström approximation the kernel k is
ks(æ, v) k(x, S) S(t k(S, v)-
(11)
Here, 1 denotes the Moore-Penrose pseudoinverse. The Nyström approximate integral operator is
Tsf == / ksl k)(,)))(o)) (d(()
This definition leads to formula for the quadrature rule Ei=1 w,k(,, with optimal weights ( Proposition 3. Fix nodes fs.....,sh} With the optimal weights w. E R" we have
Σ Wink(, 82) = =Tsg-
(12)
Consequently.
Erra = || Ts)gl S (g. (T - -T)))))
To finish the proof of theorem ! we develop and solve a recurrence for an upper bound on the largest eigenvalue of B[T Tss. See appendix A for details and for the proof of proposition 3
Efficient randomly pivoted Cholesky by rejection sampling
To efficiently perform RPCHOLESKY sampling in the continuous setting, we can use rejection sampling. (A similar rejection sampling idea is used in the MCMC continuous volume sampler of 34 Assume for simplicity that the measure μ is normalized that
[ ((, T
We assume two forms of access the kernel k:
1. Entry evaluations. For any x, y X, evaluate k(z, y).
2. Sampling the diagonal. We can produce samples from the measure k( x) du(x))
To sample from the RPCHOlESKY distribution, we use k(x,z) (( du(()) as proposal distribution and accept proposal s with probability
As 8)
k(s,8 s)
(Recall kg from
1I
The resulting algorithm for RPCHOLESKY sampling is shown in algorithm 2

---

Algorithm 2 RPCHOLESKY with rejection sampling
Input: Kernel k X x X + I and number of quadrature points v
Output: Quadrature points $1, .p......
1: Initialize L t- Onxn. i < 1-1,5
> L stores Cholesky decomposition of k(s,
2: while S do
3: Sample from the probability measure k(æ, x) du(x)
4: (d, c) KESIDUALLKENNEL(S;; L, k)
Helper subroutine algorithm
5: Draw uniform random variable U ~ UNIF(O, 1]
6:
if U < d/k(s,,,.) then
Accept with probability d/ (+; 8;)
7:
Induct 8; into the selected set: S SU {3i). i it +i11
8: L(i, 1: :1 x- [c_ Vd]
9: end
10: end while
Algorithm 3 Helper subroutine to evaluate residual kernel
Input: Point E X, set S C X, Cholesky factor L of k(S, and kernel A: Output: d = k(s, 8) ks (s, 8) and e k)))(S,,)
I: procedure RESIDCAL
2: C L kksss)
3: d+ ||||
d = k(8,3) -ks(s,
4: return c)
5: end procedure
Theorem 4. Algorithm producces exact RPCHOLESKY samples. Let tl denote the trace-error of the best rank-i approximation to
1/ A
j=++1
The expected runtime of algorithm 2 is most O (/T))) O(n /nn-))
This result demonstrates that algorithm 2 suffers from the curse of smoothness: The faster the cigenvalues 'ι, A2,.... decrease, the smaller 7h-1 will be and, consequently, the slower algorithm will be in expectation. While this curse is an unfortunate limitation, it is also a common one. The curse of smoothness affects all known Mercer-free spectrally accurate kernel quadrature schemes. In fact, the situation is worse for other algorithms. The CVS sampler of [34), for example, requires as many On log n) MCMC steps, each of which has cost O(n? //n-)) According to current analysis, algorithm
2 is O( log n)) -times faster than the CVS sampler of 34
Notwithstanding the under minimal assumptions and can reach useful accuracies r 10-5 while sampling n = 1000 nodes on laptop. Algorithm 2 may also be interesting for RPCHOLESKY sampling for large finite set X > 10S, since its runtime has no explicit dependence on the size of the space X.
of smoothness, algorithm 2 is useful in practice. The algorithm works
To achieve higher accuracies and blunt the curse of smoothness, we can improve algorithml 2 with optimization. Indeed, the optimal acceptance probability for algorithm would be
p(s;; a:) k(8;, 8;) kg(St, s()
k<x,, x) ks(x.
where a = max aEX
k(z,z)
This suggests the following scheme: Initialize with a = and run the algorithm with acceptance probability p(s;; a). If we perform many loop iterations without acceptance, then recompute the optimal a by solving an optimization problem. The resulting procedure is shown in algorithm 4
If the optimization problem for is solved to global optimality, then algorithm 4 produces exact RPCHOLESKY samples. The downside of algorithm 4 is the need to solve (typically) nonconvex global optimization problem to compute the optimal a. Fortunately, in our experiments (section 5).
only small number of optimization problems (S 10) are needed to produce sample of 1000 nodes. In the setting where the optimization problem is tractable, the speedup of algorithm 4 can

---

Algorithm 4 RPCHOLESKY with optimized rejection sampling
Input: Kernel k X x X + I and number of quadrature points v
Output: Quadrature points $1,.... Sn
1: Initialize L t- Onx.. i < 1, S t 0, trials + 0. t
2: while S do
3: trials + trials +
4: Sample s; from the probability measure k(x,, x) du(x))
5: cd < RESTDUAL S, L, k)
Helper subroutine algorithm
6:
Draw uniform random variable U UNIF(O,
7:
U < (1/a) d/k(s;, 8;) then
Accept with probability d/ //Kksi, 8:)
8:
S + (8:) ) i i 1, trials + L - c vd
9:
10:
end if
11: if trials 2 TRIALS MAX then
P We use TRIALS_ [25, 1000]
12: a+ maX RESIDUAL (,, trials +
13:
14: end if
15:
end while
be immense. To produce n. 200 RPCHOLESKY samples for the space [0, 1] with the kernel 13) algorithm&require just 0.19 seconds compared to 7.47 seconds for algorithm
Comparison of methods on a benchmark example
Having demonstrated the versatility of RPCHOLESKY for general spaces X, measures μı, and kernels in fig. [ we now present a benchmark example from the kernel quadrature literature to compare RPCHOLESKY with other methods. Consider the periodic Sobolev space Hper (l0, 1]) with kernel
Σ m 2- cos(2Tm(r v)) = 1 + (-1)a-1 (2m)) B2. - v}), (2x)
k(x, v) = 1+
m==1
where B2s is the (2s)th Bernoulli polynomial and {-} reports the fractional part of real number 3 7]. We also consider [0, 1] equipped with the product kernel
@8 (x 22, va), (y1, 92, ya)) k( 91 (2)2222))(r ( va)-
(13)
We quantify the performance of nodes S and weights uo using Err(S, g):g), which can be computed in closed form [22 Eq. (14)). We set == UNIF(O, 1]d and = 1.
We compare the following schemes:
Monte Carlo, IID kernel quadrature. Nodes s Sn iid
μ with uniform weights wy 1/n
(Monte Carlo) and optimal weights (8) (IID).
Thinning. Nodes $1...., Sn thinned from n² iid samples from μ = UNIF(O, using kernel
thinning [17 [8] with the COMPRESS+ algorithm [38] with optimal weights (8).
Continuous volume sampling (CVS). Nodes $1..... S, drawn from the volume sampling distribution Def. 1] by Markov chain Monte Carlo with optimal weights (8).
RPCHOLESKY. Nodes 81,...,8, sampled by RPCHOLESKY using algorithm with the optimal weights 8
Positively weighted kernel quadrature (PWKQ). Nodes and weights computed by the Nys- tröm-empirical-opt method as described on 22 p. 9].
See appendix D for more details about numerical experiments. Experiments MacBook Pro with a 2.4 GHz 8-Core Intel Core CPU and 64 GB 2667 MHz DDR4 RAM. Our code is available httptttt//ggiiiiie
com/ cyt.m//eeeeppeee/ooooolllleyyeee
Errors for the different methods are shown in fig. 2 (left panels). We see that RPCHOLESKY consistently performs among the best methods every value of numerically achieving the rate

---

100
10
10
10
Eay w.
Monte Carlo
Time
10-4
1
Thinming
CVS
RPCholesky (Ours)
PWKQ
10
10
10
10a
100
10 Number nodes
10
10a
nodes 12
(a)
3, d
101
10
Ear(S ():
10
Time (soc)
100
10
1
10
10.4
10
102
100
10a Number of nodes
10a
Number of modes
(b) = d
Figure 2: Benchmark example: Sobolev space. Mean quadrature error Err(S, w;g) (left) and sampling time (right) for different methods (100 trials) for (bottom). Shaded regions show 10%/90% quantiles.
1, d = 3 (top) and d = 3
of convergence of the fastest method for each problem. Particularly striking is the result that RPCHOLESKY sampling's s performance is either very close to or better than that continuous volume sampling, despite the latter's slightly stronger theoretical properties.
Figure 2 (right panels) show the sampling times for each of the nonuniform sampling methods. RPC- HOLESKY sampling is the fastest by far. To sample 128 nodes for and RPCHOLESKY was 17x faster than continuous volume sampling, 52x faster than thinning, and 170x faster than PWKQ.
To summarize our numerical evaluation (comprising fig. and fig. from the introduction), we find that RPCHOLEEEE EEKKY is among the most accurate kernel quadrature methods tested. To us, the strongest benefits of RPCHOLESKY (supported by these experiments) are method's speed and flexibility. These virtues make it possible to apply spectrally accurate kernel quadrature in scenarios where would have been computationally intractable before.
Application: Analyzing large chemistry datasets
While our focus thusfar has been on infinite domains X, the kernel quadrature formalism can also be applied a finite set X of data points. Let μ UNIF X be the uniform distribution and set g 1. Here, the task is to exhibit n nodes S = {s1, sn, C X and weights w E IR" such that the average of every "smooth" function X over the whole dataset is well-approximated by sum:
mean ( =
Σ zEX
f(x)
~ Σ w.f(())
(14)
X
Here is an application to chemistry. Let X denote large set of compounds of interest, and let X + I denote target chemical property such specific heat capacity. We interested

---

in computing mean(f). Unfortunately, evaluating f on each E X requires an expensive density functional theory computation, it can be prohibitively expensive to evaluate f on every z E X. Fortunately, we can obtain fast approximations mean()) using (4), which only require evaluating on a much smaller of size [S] [
Our experimental setup is as follows. For X, we use 2 x 10 randomly selected points from the QM9 dataset [33 37 40 We represent each compound
vector in 1 using the many-body tensor computed with the DScribe package (24]. Choose k to be a Gaussian kernel with bandwidth chosen by median heuristic [20] on a further subsample 10y random points. We omit continuous volume sampling, thinning, and positively weighted kernel quadrature because of their
representation 22
computational cost. their place, add the greedy Nyström method [19 with optimal weights (8)
Isotropie polarizability error
10
10
Ber(S un
2
Greody
1
Monte Carlo
3
IID
RPChollssk (Ours)
10
10
10
10
Number of nodes 12
Number of nodes n
(a)
(b)
Figure 3: Application: chemistry. Worst-case quadrature error (left) and mean relative error for estimation of the average value of the isotropic polarizability function f(x) (right) for different methods (100 trials). Shaded regions show 10%/90% quantiles.
Results shown in fig. Figure shows the worst-case quadrature error (9) For this metrie, IID and randomly pivoted Cholesky are the definitive winners, with randomly pivoted Cholesky being slightly better. Figure 3h shows the mean relative error for the kernel quadrature estimates of the mean isoropic polarizability of the compounds in X. For sufficiently large n, randomly pivoted Cholesky has the smallest error, beating IID and Monte Carlo by factor of three at 512.
7 Conclusions, Limitations, and Possibilites for Future Work
In this article, we developed continuous RPCHOLESKY sampling for kernel quadrature. Theorem I demonstrates RPCHOLESKY kernel quadrature achieves near-optimal error rates. Numerical results (fig. (2) hint that its practical performance might be even better than that suggested by our theoretical analysis and fully comparable with kernel quadrature based on the much more computationally expensive continuous volume sampling distribution. RPCHOLESKY supports performant rejection sampling algorithms (algorithms) and which facilitate implementation for general spaces X, measures μ., and kernels k with ease.
We highlight three limitations of RPCHOLESKY kernel quadrature that would be worth addressing in future work. First, given the comparable performance of RPCHOLESKY and continuous volume sampling in practice, it would be desirable to prove stronger error bounds for RPCHOLESKY sampling or find counterexamples which demonstrate a separation between the methods. Second, it would be worth developing improved sampling algorithms for RPCHOLESKY which avoid the need for global optimization steps. Third, all known spectrally accurate kernel quadrature methods require integrals of the form fx k<x,, (y) du(y), which may not be available. RPCHOLESKY kernel quadrature requires them for the computation of the weights (8). Developing spectrally accurate kernel quadrature schemes that avoid such integrals remains major open problem for the field.
10

---

Acknowledgements
We thank Lester Mackey, Eliza O' O'Reilly, Joel Tropp, and Robert Webber for helpful discussions and feedback.
This report was prepared as an account of work sponsored by an agency of the United States Govermment. Neither the United States Govermment nor any agency thereof, nor any of their employees, makes any warranty, express or implied. or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.
Disclaimer:
References
[1] Nima Anari, Yang P. Liu, and Thuy-Duon Vuong. Optimal sublinear sampling of spanning trees and determinantal point processes via average-case entropic independence. 63rd Annual Symposium on Foundations of Computer Science, pages 123-134. IEEE, October 2022.
[2] Nima Anari, Shayan Oveis Gharan, and Alireza Rezaci. Monte Carlo Markov chain algorithms for sampling strongly Rayleigh distributions and determinantal point processes. In Proceedings of the 29th Conference on Learning Theory, pages 103-115. PMIR, June 2016.
[3] Francis Bach. On the equivalence between kernel quadrature rules and random feature expan- sions. Journal of Machine Learning Research, 18(1):714-751, January 2017.
[4] Francis Bach, Simon Mucoste-Julien, and Guillaume Obozinski. On the equivalence between herding and conditional gradient algorithms. Proceedings of the 29th International Conference on Machine Learning, March 2012.
[5) Rémi Bardenet and Adrien Hardy. Monte Carlo with determinantal point processes. The Annals of Applied Probability, 30(1):368-417, February 2020.
[6] Ayoub Belhadji. An analysis of Lemakov-Zolotukae quadrature using kernels. In Advances in Neural Information Processing Systems, volume 34, pages 27278-27289. Curran Associates, Inc., 2021.
[7] Ayoub Belhadji, Rémi Bardenet, and Pierre Chainais. Kernel quadrature with DPPs. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[8] Ayoub Belhadji, Rémi Bardenct, and Pierre Chainais. Kernel interpolation with continuous volume sampling. In Proceedings of the 37th International Conference Machine Learning, pages 725-735. PMLR, November 2020.
[9] Daniele Calandriello, Michal Derezinski, and Michal Valko. Sampling from k-DPP without looking at all items. In Advances in Neural Information Processing Systems, volume 33, pages 6889 6899. Curran Associates, Inc., 2020.
[10] Yifan Chen, Ethan N. Epperly, Joel A. Tropp, and Robert J. Webber. Randomly pivoted Cholesky: Practical approximation of a kernel matrix with few entry evaluations, February 2023. arXiv:2207.06503
[11] Yutian Chen, Max Welling, and Alex Smola. Super-samples from kernel herding. In Proceedings of the Twenty-Sixhh Conference on Uncertainty in Artificial Intelligence, page 109-116. AUAI Press, 2010.
[12] Amit Deshpande, Luis Rademacher, Santosh Vempala, and Grant Wang. Matrix approximation and projective clustering via volume sampling- In Proceedings of the seventeenth annual SCAM symposium on Discrete algorithm, pages 1117-1126. ACM, 2006.

---

[13] Amit Deshpande and Santosh Vempala. Adaptive sampling and fast low-rank matrix approxi- mation. Josep Díaz, Klaus Jansen, José D. P. Rolim, and Uri Zwick, editors, Approximation, Randomization, and Combinatorial Optimization. Algorithus and Techniques, pages 292-303. Springer Berlin Heidelberg. 2006.
[14] Josef Dick and Friedrich Pillichshammer. Quasi-Monte Carlo integration, discrepancy and reproducing kernel Hilbert spaces, page 16-45. Cambridge University Press, 2010.
[15] Tobin A. Driscoll, Nicholas Hale, and Lloyd N. Trefethen. Chebfun guide, 2014.
[16] J. C. Dunn. Convergence rates for conditional gradient sequences generated by implicit step length rules. SIAM Journal on Control and Optimization, 18(5):473 (87)):433-887 1980.
[17] Raaz Dwivedi and Lester Mackey. Kernel thinning. In Proceedings Thirty Fourth Conference on Learning Theory, pages 1753-1753, PMLR, August 2021.
[18] Raaz Dwivedi and Lester Mackey, Generalized kernel thinning. In Tenth International Confer- ence on Learning Representations, April 2022.
[19] Shai Fine and Katya Scheinberg. Efficient SVM Training training using low-rank kernel representations. Journal of Machine Learning Research, 2(Dec):243-264, 2001.
[20] Damien Garreau, Wittawat Jitkrittum, and Motonobu Kanagawa. Large sample analysis of the median heuristic, October 2018. 1eXXv:1707.07266
[21] Guillaume Gautier, Rémi Bardenet, and Michal Valko. On two ways to use determinantal point processes for Monte Carlo integration. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[22] Satoshi Hayakawa, Harald Oberhauser, and Terry Lyons. Positively weighted kernel quadrature via subsampling. In Advances Neural Information Processing Systems, volume 35, pages 6886-6900) 6900. Curran Associates, Inc., 2022.
[23] Satoshi Hayakawa, Harald Oberhauser, and Terry Lyons. Sampling-base Nyström approxima- tion and kernel quadrature, January 2023, arXiv:230.00511.
[24] Lauri Himanen, Mare O. J. Jäger, Eiaki V. Morooka, Filippo Federici Canova, Yashasvi S. Ranawat, David Z. Gao, Patrick Rinke, and Adam S. Foster. DScribe: Library of descriptors for machine learning in materials science. Computer Physics Communications, 247: 247:106949, February 2020.
[25] Haoyan Huo and Matthias Rupp. Unified representation of molecules and crystals for machine learning. Machine Learning: Science and Technology, 3:045017, November 2022.
[26] Ferene Huszár and David Duvenaud. Optimally- herding is Bayesian quadrature. In Proceedings of the Pventy-Eightte Conference on Uncertainty Artificial Intelligence, page 377-386. AUAI Press, 2012.
[27] Hans Kersting and Philipp Hennig. Active uncertainty calibration in Bayesian ODE solvers. In Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence, page 309-318. AUAI Press, 2016.
[28] Alex Kulesza and Ben Taskar. k-DPPs: Fixed-size determinantal point processes, In Proceedings of the 28th International Conference on Machine Learning, pages 1 1190-1200. Omnipress, January 2011.
[29] C. Litterer and T. Lyons. High order recombination and an application to cubature on Wiener space. The Annals of Applied Probability, 1204):13 August 2012.
[30] Kevin P. Murphy. Machine Learning: A Probabilistic Perspective. The MIT Press, 2012.
[31] E. J. Nyström. Über Die Praktische Auflösung von Integralgleichungen mit Anwendungen auf Randwertaufgaben. Acta Mathematica, 54(0):185-204, 1930.
12

---

[32] Supratik Paul, Konstantinos Chatzilygeroudis, Kamil Ciosek, Jean-Baptist Baptiste Mouret, Michael A. Osbore, and Shimon Whiteson. Alternating optimisation and quadrature for robust control. In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), pages 3925-3933, AAAI Press, February 2018.
[33] Raghunathan Ramakrishnan, Pavlo O.
Dral, Matthias Rupp. and O. Anatole von Lilien- feld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific Data, ((1 1(1):140022, December 2014.
[34] Alireza Rezaei and Shayan Oveis Gharan. A polynomial time MCMC method for sampling from continuous determinantal point processes. In Proceedingy of the 36th International Conference on Machine Learning, pages 5438-5447. PMLR, May 2019.
[35] Christian P. Robert. The Bayesian Choice: From Decision-Theoree Foundations to Computa- lional Implementation. Springer Texts Statisties. Springer, second edition, 2007.
[36] Christian Robert and George Casella. Monte Carlo Statistical Methods. Springer Texts in Statistics. Springer, second edition, 2004.
[37] Lars Ruddigkeit, Ruud van Deursen, Lorenz C. Blum, and Jean-Louis Reymond. Enumeration of 166 billion organic small molecules in the chemical universe database GDB-17. Journal of Chemical Information and Modeling, 52(1 I 1):2864-2888 November 2012.
[38] Abhishek Shetty, Raaz Dwivedi, and Lester Mackey, Distribution compression in near-linear time. In Tenth International Conference on Learning Representations, April 2022.
[39] Ingo Steinwart and Clint Scovel. Mercer's theorem on general domains: On the interaction between measures, kernels, and RKHSs. Constructive Approximation, 35(3):363 June 2012.
[40] Annika Stuke, Milica Todorović, Matthias Rupp, Christian Kunkel, Kunal Ghosh, Lauri Hima- nen, and Patrick Rinke. Chemical diversity in molecular orbital energy predictions with kernel ridge regression. The Journal of Chemical Physics, 250(20):20412,,, May 2019.
[41] Kazuma Tsuji, Ken ichiro Tanaka, and Sebastian Pokutta. Pairwise conditional gradients without swap steps and sparser kernel herding. In Proceedings of the International Conference on Machine Learning, 2022.
[42] Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. In Proceedings of the 13th International Conference on Neural Information Processing Systems, pages 661-667 MIT Press, January 2000.
13

---

A Proof of theorem I
We start off by proving proposition 3
Proof of proposition 3 Equation
follows from a direct calculation:
w.k(+,s;) = k(,, S) / ((.5))(S..))''k((s..))) -
dp(()) Tsg.
3 $2.1] that 7-1/2 L2())) is an isometry and To establish the second claim, use the fact
invoke the definition (9h of Err(S;
T 1/2 Ts)g 2
Err Err?(S: (S:g)
(T Ts)gll
(12(K)
(T (=)))) S (g. (T Ts) (()))
The last inequality follows because T Ts and Ts T (T- are both positive semidefinite.
we employ the following notation. Let < denote the Loewner order for bounded, linear operators on L2 (m)). That is, A B B A is positive semidefinite (psd). Let k; denote the residual
For the rest of our proof of theorem 1
k k(s)......;)
for the approximate kemnel associated with the first nodes. (See proposition 8 below.) The associated integral operator
Tif = f kc(. du(r).
To analyze RPCHOLESKY, begin by analyzing single step of the procedure. The expected value of the residual kernel after one step, k1, is
Sx kr, s) du(s)
(x, = k(r.u)
Ix kk(z d(((2)
Consequently, the expected value the integral operator Ti is
2
(TT) = E
=
(15)
T
The map & enjoys several useful properties 10
Lem. 3.2].
Proposition 5. The map & is positive, monotone, and concave with respect to the order 1 on the set of trace-class psd operators on L2 (p).. That is. for trace-class psd operators A, B,
Σ A - 0 3 ( 3
(16)
and, for 0 E [0.
のΦ () + (1 ( W + (0. + - B)
(17)
Proof. We reproduce the proof from 10 Lem. 3.2] for completeness.
Positivity. Assume A 0. Observe that Id 2 A/ A, where Id denotes the identity operator. Thus
A
((A) A ld
0,
Tr
This establishes that 4 is positive: A > 0 implies ((A) > 0.
Concavity. For E [0, 1). denote 0 = 1 0. Then
00
Tr
Tr A V Tr
(0.A 0B) ( 00 (3) 0 A 0 T
1
Tr
This establishes the concavity claim 17
14

---

Monotonicity. Suppose 0 3 A = B. Observe the map + is positive homogeneous: For > 0, (0AA) (A). Thus,
A (B A) A )
( = ( + (B - A)) = 2
2
Invoke concavity
17 obtain
A + (B - A)
( 24
( + 2 - A) = ((A).
2
In the last inequality, use the positivity ( A) 0. This completes the proof of (16).
We can use the map & to bound the residual integral operator T.: Proposition 6. For i 0, have the bound
[ 3 46 (T);
(18)
where कh denotes the i-fold composition of 4.
Proof. The ith step of RPCHOLESKY applies step of RPCHOLESKY to k;-1- Thus,
E[T.] E SI...., s; ill E[
Now, we can apply concavity 17 to obtain
[7: E[ ([E(T:-)))
Applying this result twice and using monotonicity [6 yields:
E[T:] 3 & (E[T.-1] and ETT;-1] = (T B[T]] = ((((EET)-)))
Iterating this argument, we obtain the desired conclusion
(88)
We now have the tools to prove theorem We actually prove the following slight refinement:
Theorem 7. Let 5 = {si,.... 8n) be generated by the RPCHOLESKY sampling algorithm. For any function g L (f)), nonnegative integer O, and real numbers a, >> 0, we have
E EEE < (191l2a())
aat8 Σ A
1=+#1
provided
A1
1
P log
(19)
Theorem immediately follows from theorem by taking = = 8/2 and a = (6/2) Ei do We
now prove theorem 7
Proof of theorem 7 In the standing notation, 7 - Ts corresponds to Tn- Thus, by proposition 6 we have
Bgg, (T - T5) = (g, ((T,))) S (g, (T) S ( Amax (D"(T))
Thus, it is sufficient to show that
Amax (T)(T) < +
X
(20)
i=++1
holds under condition (19
Using the definition 15 the following recurrence
of कb we see that the eigenvalues ( 2 1(0) 2 ... (T() are given by
2001 = 1 (i-1)
for i,j = 1,2, ...
(21)
Exe (0-1)
15

---

with initial condition 2/ Ay for all j. From this formula, we immediately conclude that each 215 is nonnegative and nonincreasing in In addition, the ordinal ranking X (8) for each Indeed, if X 1 2 (i-1) 2 0, then
2 2 is preserved
2 241 = =Aft1 - 1 ( )
EE=1 2/ 11-11)
(1
A 1 + 1 x
j+1 - -
2 0, Et=l A (1-1)
and the claim follows by induction on i. To bound Amax Famxx(D'(11 (T'(T) 1 bound the sum of eigenvalues appearing in the recurrence
Σ t=1
= x
11-1) + Σ x
(1-1) sra (1-1)
de-
(22)
(=++1
8=++1
To bound the first terms, we used the ordering of the eigenvalues X 1(1-1)
( To bound
the tail, we used the fact that for each E {r + +, ..... the eigenvalue A (1- 1) in and thus bounded by 2 = dę. Substituting (2 21). we obtain
is nonincreasing
(1/ (1-11))
1 < 11
(i-1)
for i 1,2,
(23)
(( Etert1
To bound the solution to the recurrence
23 we pass to continuous time. Specifically, for any / E Z4p.
X is bounded by the process x+)) solving the initial value problem
xt)²
(1) r1( + oo
xx0) = A1.
di
Xe"
The bound 1 < xtt) holds for all E 24 because + -22/(rx + Li-r+1 de) is decreasing
for v E (0, oo). We are interested in the time = l, at which x(44) = a t Ei=rt1 de which we obtain by separation of variables
tts Ei de
TI + 2t=r+1 de
Σx 21sr41 de
s=a+e Eierts d
dx = -rlogæ +
t,
x²
A1
Eler41 de
Eler41 de
d1 <rlog
rloog & Eie de
a a Eesr41 de.
A1
Since z is decreasing, the following holds for satisfying 10
Amaxe (T)) = ((m) < x(n)= S xtt.) a t Σ N
(=1+1
This shows
20
completing the proof.
B Proof of theorem
To see why algorithm 2 produces samples from the RPCHOLESKY sampling distribution, observe that each effective iteration the matrix L 1,1: 1) is the Cholesky factor of ((S, Therefore, d is equal to the residual kernel kges
= k ks evaluated at 8;
d = kg rrs st).
16

---

At each effective iteration i, the candidate point 8; is sampled from k(z, 2) du((a) and accepted with probability (s;, (), s:). Therefore, conditional on acceptance and previous steps of the algorithm, is sampled from the probability measure
Ks (2,
ks (z, x) du(x)
k(x, x) du(x))
P P(Aceptaance S}
k(,,))
P(Acceptance S}
Since this is probability measure which must integrate one, the probability of acceptance is
PP PPAcceptance [} / h (x.
Therefore, conditional on acceptance and previous steps of the algorithm, & is sampled from the probability measure (z, dp(z)/ [x (x, du(x). This corresponds exactly with the RPCHOLESKY sampling procedure described in algorithm ! (see proposition teatiiionnnbloooo beloow.
Now, we bound the runtime. For each iteration i, the expected acceptance probability is bounded from below by the smallest possible trace error achievable by a rank (i 1) approximation
[
P(Acceptance) [[cceeaaaa [}]] EE
kge (z. x) dd(a)) 77-1-
Therefore, the expected number of rejections at step i is bounded by the expectation of a geometric random variable with parameter V/i-1, i.e., by 1--1- Since each loop iteration takes O operations (attributable to line 4) and at most n; expected iterations are required to accept at step i, the total expected runtime most O i= 82 //:-11) 011//m1-1)
C Derivation of 7
For completeness, we provide a derivation of 7 This result is standard; see, e.g., 3 definition 6 for Err(S, w:
Eq. (7)). Recall
( (()))()))paa)
Err(S, w;gg)
max |
Σ 1=1
wi (a]((())
Use the reproducing property and linearity to write
/ (U.K du(z)
Err(S, w;g) max Illl $1
Σ wil k())
Now apply linearity of (1,)) obtain
( (((,z))9())) dp(z))))) dp(x) - Σ (u:k(,,s))) si)
Err(S, w: g) = max
|
1=1
Use the dual characterization of the RKHS norm ||-1| and the definition
5 of I'' to conclude
Err(S, w:g) = Tg D wPk( 8;)
D Details for numerical experiments
In this section, we provide additional details about our numerical experiments.
Continuous volume sampling. We are unaware of provably correct n-DPP/continuous volume distribution sampler that is sufficiently efficient to facillitate numerical testing and comparisons. The best-known sampler (34) for the continuous setting inexact and requires O(n log n) MCMC steps, which would be prohibitive for us. To produce samples from the continuous volume sampling distribution for this paper, we use a continuous analog of the MCMC algorithm (2. Given a set
S {s,,..., } and assuming is probability measure, one MCMC step proceeds follows:
17

---

1. Sample Snew ~ μ and sold ~ UNIF S.
2. Define SU {Snew) {sola}.
3. With probability
det k (') det k(S,S) 5 accept and set + S'. Otherwise, leave unchanged.
2 min {
In our experiments, we initialize with , drawn iid from and run for 10n MCMC steps.
Running for just 10n MCMC steps is aggressive, even in the finite setting. Indeed, the theoretical analysis of 2 Thm. 2, Eq. (1.2)] for μ = 2, [X]]] proves the sampler mixes to stationar- ity in at most (x) MCMC steps when properly initialized. Even with just 10n MCMC steps, continuous volume sampling is dramatically more expensive than RPCHOLESKY in our experiments for n 2 64.
Finally, notice that the performance of continuous volume sampling degrades that of iid sampling for sufficiently large n. We believe that the cause for this are numerical issues in evaluating the determinants of the kernel matrices k(S, S). This suggests that, in addition to being faster, RPCHOLESKY may be more numerically robust than continuous volume sampling.
Thinning and positively weighted kernel quadrature. Our kernel quadrature experiments with thinning and PWKQ based on code from the authors of [22] available online at httttps://gttuub com/ satoshi Saaoyyaawaa/kereeeee quarratuuuurruurrr which uses the goodpoints package (https http://coooooooooooooorott/ gowptt/iooooiinntsoooo by the authors of [8 38 to perform thinning. We
//github. com tttp:/miccoottttteeeeeecooo
use 4, 0.5, and four bins for the COMPRESS+ algorithm.
Miscellaneous.
Our code available online
To evaluate exact integrals
Tg(x) = / k(x, uy) du(y)
for fig. L we use the 2D functionality in ChebFun [15]. To compute the optimal weights 8 we add a small multiple of the identity to regularize the system:
W reg = (k(S,S) S) 10 trace(k(S, -)) I)g())
Here, Emach = 2-
52 2-52 is the double precision machine epsilon.
E Randomly pivoted Cholesky and Nyström approximation
In this section, we describe the connection between the RPCHOLESKY algorithm and Nyström approximation. We hope that this explanation provides context that will be helpful in understanding the RPCHOLESKY algorithm and the proofs of theorems I and 4
To make the following discussion more clear, we begin by rewriting the basic RPCHOLESKY algorithm (see algorithm D in algorithm 5 avoid overwriting the kernel & at each step:
From this rewriting. we see that in addition selecting quadrature nodes s1...., Sp. the RPC- HOOLESKY algorithm builds an approximation hn to the kernel kt Indeed, one can verify that the approximation kn produced by this algorithm is precisely the Nyström approximation (definition 2
Proposition 8. The output kn of algorithm is equal to the Nyström approximation k().......)
Proof. Proceed by induction n. The base case (S1,..., Sn s' = Sat1, and S' = {si,..., Sn, s'}. Write the Nyström approximation ky. as }
is immediate. Fix and let
S(
k(S, k(s',8)
ks = [k(.)) k(,
k(8' S)
18

---

Algorithm 5 RPCHOLESKY: unoptimized implementation without overwriting
Input: Kernel k X x X + I and number of quadrature points v
Output: Quadrature points 81, 81,-.., Sn and Nyström approximation kin = (n to k
Initialize kn <- kig +-
2: 1or1=1,2 ................................ do
3: Sample *; from the probability measure k(x, a) du(x)/ Ix k( x) du(a))
4: Update Nyström approximation k < k + ki-1(: 8, 8,) - (86,)
5: Update kernel k; f k 1 ( ())1()
6: end for
Compute block LDL"
factorization of the matrix in the middle of this expression:
((S,S) k(s,
()s,,)) !]
I
()))
I
S(s ((s,,)]
0
(s'
0
1
Here, we've used the induction hypothesis
kn'd', = k's', - ks 80) - - -
(24)
Thus,
k(S,S)
k(S, 22
k((S
:]
ks. = S):(,)) k(),8')
S)
k(6',5
A(s',-)
I
S), k(S, $')
0
S)( s(),
O
1
0
kn bn(s'ss')~1
[-b(d, (),((,
I
0
R((S :)
1
k(8)
(K(S,S)-1
A(S, kn(-,
()) (=
(25)
0
kn kn's',
kn + (. s') kkn (s', s')-''n (s'. knt1.
For 25
used the induction hypothesis (4 again. Having established that Ks- kn+1. the proposition is proven.
19


Reviews:
[{"review_id": "5nMGXLQ5ph", "reviewer": "Reviewer_dsgV", "summary": "This article studies a family of quadrature rules suited for functions that belongs to an RKHS. The proposed construction is based on an algorithm called randomly pivoted Cholesky (RPCholesky). The latter was studied for the task of approximating positive-semidefinite matrices. In particular, the authors extended RPC to continuous domain, and investigated its effectiveness for kernel quadrature. More specifically, the authors studied the quadrature rule obtained by taking the nodes to be the configuration of points that RPCholesky and by taking the weights that minimises the worst-case integration error on the unit ball of the RKHS (WCE). In particular, they proved that the expected value of the squared WCE of this quadrature rule converges to 0, with a rate of convergence that depends on the eigenvalues of the integration operator. Moreover, they studied the expected runtime of the proposed algorithmes.\n\nThe contributions of the article may be summarised as follows: \n* An extension of the RPCholesky algorithm to a continuous domain (Algorithm 1) along with a detailed implementation of it using rejection sampling (Algorithm 2).\n* An implementation of RPCholesky with optimized rejection sampling (Algorithm 4).\n* A generic result (Theorem 1) gives the upper bound of the expected squared WCE for the proposed quadrature rule.\n* The study of the expected runtime of Algorithm 2 (Theorem 4). \n\n", "questions": "- The empirical comparison in section 5 was limited to Sobolev spaces of order s=3. It may be interesting to compare RPCholesky to other algorithms on other RKHSs. \n- The empirical experiments make use of a continuous analog of a ’basis-exchange’ MCMC algorithm proposed in [2] instead of the MCMC algorithm based on Gibbs sampling proposed in [1]. Up to my knowledge, only the latter was studied theoretically, and it is possible that the convergence of ’basis-exchange’ MCMC algorithm is slower than the Gibbs sampling.", "limitations": "None", "rating": 4, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- The article is very well written. \n- The article expands the application of the randomly pivoted Cholesky algorithm into the continuous domain, establishing another bridge between the numerical linear algebra community and the kernel quadrature community.\n- The proof of the main result (Theorem 1) is original and may be beneficial for the community of kernel quadrature in the future.\n", "weaknesses": "At first sight, the proposed algorithm seems to bring a significant algorithmic improvement upon continuous volume sampling. This is at least what we may understand reading the comments of Figure 2 (line 235-line 238). While, the relative improvement is undeniable, RPCholesky suffers from the same limitation of continuous volume sampling that is poor scalability for large values of $N$, especially for smooth kernels. Figure 2 of the article is corroborating this claim: the dependence of the empirical expected time as a function of $N$ seems to be polynomial. This is particularly the case in dimension 1,  where we have $\\lambda_m = \\mathcal{O}(m^{-2s})$ so that $1/\\sum_{m=N+1}^{+\\infty}\\lambda_{m} = \\mathcal{O}(N^{2s+1})$. In other words, the smoother is the kernel the longer it will takes to sample from the RPCholesky, which is exactly the same problem with Gibbs sampling proposed in [1], which provably approximate continuous volume sampling. It turns out that Algorithm 1 is nothing but the algorithm proposed in [1] to initialise an MCMC algorithm that approximate continuous volume sampling; see Algorithm 1 in [1].  The difference between the two works is in the subroutine that consist in sampling from the probability measure $x \\mapsto k(x,x) \\mathrm{d}\\mu(x)/\\int_{\\mathcal{X}}k(u,u)\\mathrm{d}\\mu(u)$: [1] proposed to use a ‘stochastic greedy’ algorithm using rejection sampling, while this article proposed to use RPCholesky with rejection sampling. Interestingly, the two algorithms have an expected time that depends on $1/\\sum_{m=N+1}^{+\\infty}\\lambda_{m}$, where is the number of nodes in the quadrature rule, see Lemma 5.1 in [1]. In other words, this work does not break the ‘curse of smoothness’ that hinders the scalability of continuous volume sampling. \n\n[1] Rezaei, A., & Gharan, S. O. (2019, May). A polynomial time MCMC method for sampling from continuous determinantal point processes. In International Conference on Machine Learning (pp. 5438-5447). PMLR.\n\n[2] Anari, N., Gharan, S. O., & Rezaei, A. (2016, June). Monte Carlo Markov chain algorithms for sampling strongly Rayleigh distributions and determinantal point processes. In Conference on Learning Theory (pp. 103-115). PMLR.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's very incisive comments. We are delighted to hear that they found the article well-written and that they enjoyed the further connections our article draws between kernel quadrature and numerical linear algebra. Below, we address your concerns and questions:\n\n## Curse of Smoothness\n\nWe agree that Algorithm 2 does not break the curse of smoothness. However, RPCholesky has a substantially better scaling with $n$ compared to CVS. [1, Thm. 1.3] suggests that $O(n^5 \\log n)$ MCMC steps are needed to sample from the CVS process. For kernels with polynomial spectral decay, RPCholesky provably produces near-optimal quadrature nodes with just $n$ steps. **The improvement in runtime over [1] is log-quartic $O(n^4 \\log n)$.**\n\nAs the reviewer suggests, RPCholesky is identical to the \"warm start\" used to initialize the CVS MCMC sampler of [1]. We believe the connections between RPCholesky and DPP/CVS to be a feature, not a bug. One way to frame the results of this paper is that the initialization phase of [1] is already enough to get near-optimal quadrature nodes, and the costly $O(n^5 \\log n)$ MCMC steps required to sample from the CVS distribution can be avoided.\n\nWe also bring the reviewer's attention to Algorithm 4. This algorithm blunts the curse of smoothness by incorporating a small number of global optimization steps to choose a better threshold for rejection sampling. There are certainly limitations to Algorithm 4: We do not have theory for this variant, and the requirement for optimization limits the range of applicability. However, in many settings where the curse of smoothness really starts to bite, we hope Algorithm 4 will be a useful tool in practice (see below for more discussion on this point).\n\n## Gibbs Sampler vs Basis Exchange\n\nOur paper used the basis exchange sampler rather than the Gibbs sampler because a direct implementation of the Gibbs sampler as described in [1] for our problem would be intractable. Indeed, using a direct implementation of [1, Alg. 3] and running for an optimistic $10n$ MCMC steps (rather than the theoretically sanctioned $O(n^5 \\log n)$ steps [1, Thm. 1.3] or the empirically suggested $O(n^2)$ steps [1, Fig. 1]), **the Gibbs sampler required 30 minutes just to produce a single sample of size 20 for the $s=3$, $d=1$ Sobolev space. Algorithm 4 requires just 0.1 seconds to produce an RPCholesky sample of this size.** We hope this example provides significant evidence to the reviewer that, while we have not provably conquered the curse of smoothness, Algorithm 4 can improve the state of affairs significantly in practice.\n\nWe can improve the runtime of the Gibbs sampler of [1] by the ideas of Algorithm 4 to sample the conditional 1-DPP subproblems needed to implement the algorithm. Specifically, at each MCMC step, we compute an upper bound on the residual kernel using a global optimization step and use this to determine the optimal \"$M$\" parameter for rejection sampling. We again run for $10n$ MCMC steps. Results are shown in Fig. 2 in the attached PDF. The quadrature error with CVS sampled by the Gibbs sampler is worse than both RPCholesky and for CVS sampling using basis exchange. **Even with the optimistic choice to run for only $10n$ MCMC steps and the algorithmic improvements adopted from Alg. 4, the Gibbs sampler takes 120$\\times$ longer than RPCholesky to sample $n=64$ nodes**\n\n## Conclusion\n\nWe appreciate the reviewer's thoughtful comments, which demonstrate both a close reading of our paper and a lot of knowledge of the related literature. We hope that the above discussion is convincing that RPCholesky results in a substantial acceleration over the existing state-of-the-art spectrally accurate kernel quadrature schemes in both theory and practice."}, {"review_id": "SCacvA8pvR", "reviewer": "Reviewer_4xjP", "summary": "The recently proposed Randomly Pivoted Cholesky procedure is an efficient and accurate method for approximating symmetric PSD (e.g. kernel) matrices. It builds a Nystrom approximation via randomizing the choice of pivot in the classical Cholesky decomposition procedure.\nThe authors of the present submission transport the same sampling idea to the continuous case and develop new kernel quadrature rules.\nThe proposed algorithm is based on rejection sampling, it's simple to implement and it's much faster than many other kernel quadrature procedures. Rigorous theoretical guarantees and experimental analysis with synthetic data accompany the algorithm.", "questions": "1) The analysis and experiment focus on and minimize the worst case error (line 140). What, if anything, would change if you considered some expected error (loss) when $f$ is drawn from some distribution?\n2) Since the experiments in Figure 2) are rather low dimensional ($d=1$ and $d=3$) could you also include a Quasi-Monte Carlo baseline with equal and optimal (eq (5)) weights similarly to current Monte Carlo and IID? Since QMC is not spectrally accurate one would still expect RPCholesky to have much lower error, the gap might shrink though. Could you also consider running experiments with much higher $d$?\n", "limitations": "Theory paper, societal impact is not applicable. The authors do a thorough job in disclosing and discussing technical limitations.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1) Solid theory. The error of the quadrature depends on the eigenvalue distribution of the kernel (it's spectrally accurate).\n2) The rejection sampling algorithm is practical to implement and it's efficient.\n3) Empirical evaluation supports the claims.", "weaknesses": "1) As the authors themselves point out in Section 6 all spectrally accurate quadratures, including this, currently require the exact computation of the Hilbert–Schmidt integral of function $g$, which may limit the applicability of the method. E.g. the experiments are all based on $g\\equiv 1$.\n2) Although the first sentence of the introduction in lines 12-13 refers to several applications all the experiments in Figures 1 and 2 rely on computing the worst case error of kernel integrals in low dimensions on a given domain with a particular (complicated) kernel. Albeit this style of synthetic evaluation seems common in the cited prior work I believe it would increase the impact of the paper if the authors demonstrated end to end improvements on an actual machine learning tasks that requires computing these integrals. Including a concrete application would boost the accessibility and motivation for non experts reading the paper as well in my opinion.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks to the reviewer for their considered remarks. We are very glad to hear the reviewer found the theory and algorithms of the paper to be solid and practical. Below, we address the reviewer’s questions and suggestions:\n\n## Application \n\nAt the suggestion of this review and others, we have included a new experiment to the paper (Fig. 1) featuring a concrete application with high (1,500)-dimensional chemistry data. We agree with the reviewer that including this example will help make the paper more interesting and accessible to a wider audience. We thank the referee for this suggestion.\n\n## Choice of Loss Function\n\nThe reviewer asks what would change if considered the expected loss over a random function $f$. Fortunately, there is a quite satisfactory answer. Let $f$ be a Gaussian process with covariance $k$ and consider the loss$$\\tilde{\\mathrm{Err}}^2(\\mathsf{S},\\boldsymbol{w};g) \\coloneqq \\mathbb{E} \\left( \\int f(x) g(x) \\, \\mathrm{d}\\mu - \\sum_{i=1}^n w_i f(s_i) \\right)^2.$$Then $\\tilde{\\mathrm{Err}}^2(\\mathsf{S},\\boldsymbol{w};g) = \\mathrm{Err}^2(\\mathsf{S},\\boldsymbol{w};g)$ is exactly the same as the worst-case loss defined in the paper. Here is a brief sketch of the proof. The Gaussian process $f$ has a Karhunen–Loève expansion $f = \\sum_i Z_i (\\sqrt{\\lambda_i} e_i)$ in terms of the $\\langle\\cdot,\\cdot\\rangle_{\\mathcal{H}}$-orthonormal basis $(\\sqrt{\\lambda_i} e_i)$, where the $Z_i$ are iid standard normal. Writing $\\int f(x) g(x) \\, \\mathrm{d}\\mu - \\sum_{i=1}^n w_i f(s_i) = \\langle f, Tg - \\sum_{i=}^n w_i k(\\cdot,s_i) \\rangle_{\\mathcal{H}} \\eqqcolon \\langle f, \\ell \\rangle_{\\mathcal{H}}$, we can expand $\\ell = Tg - \\sum_{i=1}^n w_i k(\\cdot,s_i) = \\sum_j a_j (\\sqrt{\\lambda_j} e_j)$ in this orthonormal basis as well. Thus, by orthonormality of  $(\\sqrt{\\lambda_i} e_i)$,$$\\tilde{\\mathrm{Err}}^2(\\mathsf{S},\\boldsymbol{w};g) = \\mathbb{E}\\langle f, \\ell\\rangle_{\\mathcal{H}}^2 = \\mathbb{E} \\left( \\sum_i a_i Z_i \\right)^2= \\sum_i a_i^2 = \\lVert\\ell\\rVert_{\\mathcal{H}}^2 = \\mathrm{Err}^2(\\mathsf{S},\\boldsymbol{w};g).$$\n\n## Quasi-Monte Carlo\n\nSee the attached Fig. 3 (left panel) for the $s = d = 3$ example with the quasi-Monte Carlo baselines included. With the optimal weights, the performance is quite similar to randomly pivoted Cholesky and continuous volume sampling. Our results agree with [19, Fig. 1], which show optimally weighted QMC to be among the best performers for their 2D and 3D experiments.\n\nTo concur with one of the referee’s other remarks, we agree that these low-dimensional Sobolev examples are an odd choice to have become such standard benchmarks in the kernel quadrature literature, as existing quadrature methods (optimally weighted QMC, Gaussian, Clenshaw–Curtis, etc.) are already quite effective for these low-dimensional examples. Quasi-Monte Carlo integration is well-known to become less effective for dimensions larger than, say, six to ten (Morokoff and Caflisch), and would thus have no hope of addressing the 1,500-dimensional chemistry application in Fig. 1. \n\nMorokoff, W. J., & Caflisch, R. E. (1995). Quasi-Monte Carlo Integration. Journal of Computational Physics, 122(2), 218–230. https://doi.org/10.1006/jcph.1995.1209"}, {"review_id": "OUAzj7UeLa", "reviewer": "Reviewer_yiWH", "summary": "The authors provide theoretical analysis and show that RPCholesky has near optimal error rates.\nMeanwhile, the authors propose to combine RPCholesky with rejection sampling in order to use it in the continuous setting. Experiments show that the proposed method is compared favorably to continuous volume sampling which is much more difficult to sample.", "questions": "Can the authors include additional experiments and show that the proposed method can also perform well in practice, e.g. real-world experiments? Is this something difficult to implement?", "limitations": "The limitations are discussed in Section 6.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is clearly written.\n\nBoth assumptions and theorems are well presented and discussed.\n\nThe experiments demonstrate the good performance of the proposed method.", "weaknesses": "For ''Err'' under line 140 and ''Err'' under line 141, it’s better to include the derivation in the appendix to make the paper self-contained.\n\nThe description in Line 226 and Line 227 is unclear. For instance, “opt” is imprecise.\n\nThe proof of the theorem in Section 3.4 can be moved to the appendix. It seems not necessary to include it in the main paper.\n\nI think the biggest weakness of the paper is: the experiments are very limited. I would like to see more experiments as the authors claim that their method is very general.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are glad the reviewer found the paper well-written and found the experiments to “demonstrate the good performance of the proposed method”. We thank the reviewer for their suggestions for improvement, all of which we have accepted and incorporated into the main text:\n\n- We now provide a derivation of the Err formulas on lines 140–141 in the appendix.\n- We provide a precise description of the exact PWKQ method we use (lines 226–227).\n- We have moved the proof in section 3.4 to the appendix.\n- We have added an application experiment, demonstrating the utility of RPCholesky kernel quadrature for a large-scale data analysis problem.\n\nTaken as a whole and with the support of the new experiments (Fig. 1 in the attached PDF) suggested by the reviewer, our numerical results show that RPCholesky is consistently one of the most accurate kernel quadrature methods. RPCholesky kernel quadrature can be orders of magnitude faster than existing state-of-the-art spectrally accurate kernel quadrature schemes, and these benefits persist to a high (1,500)-dimensional chemistry data analysis task. We believe that, with the addition of the reviewer’s suggested application example, our paper makes a convincing demonstration that RPCholesky kernel quadrature represents a meaningful contribution to the kernel quadrature field."}, {"review_id": "xCvs3SIB6r", "reviewer": "Reviewer_D4rs", "summary": "This paper introduces a new method for kernel quadrature which relies on a randomly pivoted Cholesky decomposition.\nThe sampling of the quadrature points is done thanks to a rejection sampling algorithm. To improve performance, an optimized rejection sampling procedure is also introduced, so that the acceptance probability is updated if the number of rejections becomes too large. Numerical simulations show that this method is fast and often as accurate as other performant methods.\nAn error bound for the expectation of the quadrature method is also derived.", "questions": "- Concerning the simulations reported in Fig 2 (a) (left-hand side), I wonder if there is a numerical issue in the computation of the curve corresponding to continuous volume sampling (CVS). It seems to me that CVS should have at least a similar performance as RPCholesky.\nHowever, for the last two data points, CVS suddenly has a worse performance (following the IID case). This phenomenon is not happening in Fig 2 (b) (left-hand side). Hence, I wonder if there could be an issue in the implementation of CVS such as a numerical instability. Could you comment about that ?\n", "limitations": "The authors mention possible limitations.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- This paper is very well written. It mentions other important methods for kernel quadrature.\n- The main idea of this research work is simple and elegant. \n- The theoretical results are well introduced and the proofs are clear although non-trivial.\n- From the practical perspective, the numerical simulations are convincing. There is a clear gain in sampling time, in particular with respect to continuous volume sampling which is a closely related method.", "weaknesses": "There are very few weaknesses in this paper.\n- It seems that the proof techniques are mainly inspired from [10] to which the authors refer. This is not an issue and does not limit the significance of this contribution.\n- The end of the proof of proposition 3 is not very explicit. However, the authors provide the main elements to finish the proof so that the reader can complete the argument on a sheet of paper.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thoughtful feedback, and we are delighted that they appreciated the paper. Below, we address the reviewer’s comments and questions:\n\n## Proof of Proposition 3\n\nAt the suggestion of yiWH, the proof of proposition 3 has been deferred to the supplementary material. We provide a more explicit proof to make this result easier to read and understand.\n\n## CVS Numerical Performance\n\nThe reviewer very astutely points out that CVS sampling degrades to the error of IID sampling in Fig. 2(a). We observed this too and found it very curious. Due to space reasons, we commented on this in Appendix C in the supplementary materials, writing:\n\n> Finally, we notice that the performance of continuous volume sampling degrades to that of iid sampling for sufficiently large $n$. We believe that the cause for this are numerical issues in evaluating the determinants of the kernel matrices $k(\\mathsf{S},\\mathsf{S})$. This suggests that, in addition to being faster, RPCholesky may be more numerically robust than continuous volume sampling.\n\nNumerical issues with DPP samplers, continuous and finite, inexact and exact, seem to be a widespread issue. In our experience with both our implementation and DPPy, DPP samplers frequently fail or suffer performance degradations due to numerical issues. It is possible these numerical issues could be resolved with further research attention."}, {"review_id": "viYqMB39qa", "reviewer": "Reviewer_PPmk", "summary": "The kernel quadrature problem computes quadrature weights and nodes that minimize the maximum quadrature error over all functions in RKHS with the norm bounded by 1. This paper proposes new quadrature rules for functions in RKHS using nodes drawn by randomly pivoted Cholesky algorithm. The randomly pivoted Cholesky used in the paper is an extension of the finite variant introduced in [10]. The authors provide the unoptimized implementation followed by the variant with efficient rejection sampling.", "questions": "1. For the right side of the Figure 2 (time vs number of nodes): what is the error level considered for all methods?\n\n2. Could the authors comment on how the results could change for other Mercer kernels?\n\n3. How about the effect of dimensionality on the empirical results? The reference [19] does offer more parameter settings used for evaluation.\n\n4. Seeing that Algorithm 2 and Algorithm 3 stores Cholesky factors, could the authors comment on scalability of the method as the number of quadrature points is increased?", "limitations": "There are no results on real datasets and even on synthetic datasets, the results are limited to one type of kernel. Although it is clear the proposed method has wide applications in machine learning, the authors should demonstrate a usecase in which the proposed method can be applied (e.g. scaling kernel machines).", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Quadrature problem is one of the fundamental problems in computational mathematics and the paper proposes an efficient method for it.\n\n2. The experiments offer comprehensive comparisons against other state-of-the-art methods.\n\n3. Theorem 1 provides error bounds which achieve spectral accuracy with just logarithmic number of nodes.\n\n4. The proposed method (both the unoptimized and the optimized variants) is easy to implement from the description.", "weaknesses": "1. There are no empirical results on real datasets despite earlier work such as [19] having results on real datasets.\n\n2. Discussion on the application of the proposed method in machine learning is largely missing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful reading and their many insightful questions and suggestions. Our responses are as follows:\n\n## Applications\n\nThe new version of the paper discusses applications more and includes a new experiment demonstrating an application to analyzing large chemistry datasets. (See the response to all reviewers for a description of our experiment and Fig. 1 in the attached PDF for the results.) We hope this example is both simple (so that it may be appreciated by a wide audience) and convincing (the QM9 dataset we use is considered a standard, but challenging, benchmark in chemical machine learning). Many of the other applications of kernel quadrature to machine learning (such as the applications we cite to reinforcement learning [27]) are sufficiently intricate that it would be difficult to briefly summarize, but we mention them and provide references.\n\n## Questions\n\nBelow, we answer each of the reviewer’s questions:\n\n> For the right side of the Figure 2 (time vs number of nodes): what is the error level considered for all methods?\n\nWe are not certain exactly how to interpret the reviewer’s question. We will do our best to answer based on our interpretation of the question (and hope that we can provide a better answer, if needed, in the discussion period). For a given number of nodes, Fig. 2 (left) shows the error and Fig. 2 (right) shows the runtime. We do not specify a particular error level for Fig. 2 (right). \n\n> Could the authors comment on how the results could change for other Mercer kernels?\n\nAs we demonstrate in Thm. 1, the performance of RPCholesky kernel quadrature depends on the eigenvalues of the kernel (as revealed by its Mercer decomposition (2)). For generic domains equipped with the Matérn (Fig. 1) or Sobolev (Fig. 2) kernels, the eigenvalue decay is polynomial and thus the quadrature error decays at a polynomial rate as well. For a Gaussian kernel (Fig. 1, attached PDF) the generic eigenvalue decay is exponential and thus the quadrature error decays at an exponential rate. In high dimensions, the asymptotic convergence regime may take exponentially large numbers of quadrature nodes to reach, and the kernel eigenvalues can depend sensitively on the ground space $\\mathsf{X}$. Whatever the eigenvalues are, RPCholesky achieves near-optimal quadrature error and is, in our empirical tests, the fastest spectrally accurate kernel quadrature scheme.\n\n> How about the effect of dimensionality on the empirical results? The reference [19] does offer more parameter settings used for evaluation.\n\nWhile the reference [19] considers a few more parameter settings for the Sobolev example, they are all restricted to the same values $s,d\\le 3$ that we consider in the present paper. It is difficult to compare all methods for large dimensions $d$, as larger values of $n$ are needed before the error starts to show any decrease. These larger values of $n$ are hard to reach using more expensive sampling methods like thinning and PWKQ. Indeed, we did some experiments with $d=5$ and $s=8$ using a subset of the methods; see the middle and right panels of Fig. 3 in the attached PDF. RPCholesky is noticeably—if modestly—better than CVS and IID. Given the time constraints of the NeurIPS response window, we were unable to include the more time-expensive thinning and PWKQ examples for this problem.\n\nWe also draw the reviewer’s attention to our new application example (Fig. 1 in the attached PDF). This demonstrates that RPCholesky kernel quadrature extends gracefully to a 1,500-dimensional application example.\n\n> Seeing that Algorithm 2 and Algorithm 3 stores Cholesky factors, could the authors comment on scalability of the method as the number of quadrature points is increased?\n\nAll optimally weighted kernel quadrature methods require the solution of a linear system of equations $k(\\mathsf{S},\\mathsf{S})\\boldsymbol{w}_\\star = Tg(\\mathsf{S})$, requiring a Cholesky factorization of $k(\\mathsf{S},\\mathsf{S})$ at $O(n^2)$ space and $O(n^3)$ time. So, by this metric, RPCholesky’s resource requirements are no more than any optimally weighted kernel quadrature scheme. In our review of the kernel literature, existing investigations have focused on $n\\lessapprox 1000$, so the requirements of storing and factoring $k(\\mathsf{S},\\mathsf{S})$ have not come to be a true bottleneck.\n\nShould we be interested in pushing to larger values of $n$, we suspect that RPCholesky kernel quadrature (and indeed other optimally weighted kernel quadrature schemes) could be accelerated by using fast solvers developed in the kernel algorithms community. In place of inversion or full Cholesky factorization of $k(\\mathsf{S},\\mathsf{S})$, natural approaches in low dimensions are hierarchical low-rank decompositions (Ambikasaran et al.) or sparsified Cholesky factorizations (Schäfer, Katzfuss, and Owhadi). For high dimensions, Nyström preconditioning (Frangella, Tropp, and Udell) could be an effective tool to compute the action of $k(\\mathsf{S},\\mathsf{S})^{-1}$ (which is all that is needed to implement RPCholesky kernel quadrature).\n\nAmbikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. (2016). Fast Direct Methods for Gaussian Processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(2), 252–265. https://doi.org/10.1109/TPAMI.2015.2448083\n\nFrangella, Z., Tropp, J. A., & Udell, M. (2021). Randomized Nyström Preconditioning. ArXiv:2110.02820 [Cs, Math]. http://arxiv.org/abs/2110.02820\n\nSchäfer, F., Katzfuss, M., & Owhadi, H. (2021). Sparse Cholesky Factorization by Kullback–Leibler Minimization. SIAM Journal on Scientific Computing, A2019–A2046. https://doi.org/10.1137/20M1336254. "}]

Author rebuttal:
We thank the reviewers for their very thoughtful comments. We have revised the paper to address any questions and concerns raised and have accepted all specific recommended changes by all reviewers. We have responded to the questions of each reviewer in the individual responses, and we look forward to addressing any lingering concerns in the discussion period.

Several reviewers asked us to add an additional example, with different reviewers requesting an example with a higher dimension, a different kernel, or an application. In the interest of conciseness, we have addressed these concerns jointly by adding a new section to the paper featuring an application to large chemistry datasets.

The experimental setup is as follows: Following inspiration from [19], we apply the kernel quadrature formalism to a finite dataset, using kernel quadrature to estimate the average value of a chemical property over a subsample of 20,000 molecules from the QM9 dataset. Each molecule is represented as a vector in **1,500-dimensional space** using the many-body tensor representation, a standard featurization in chemical machine learning. We use a Gaussian kernel; with this new application, our paper now includes the Matérn, Sobolev, and Gaussian kernels. We show both the worst-case error and, following the suggestion of 4xjP, the error for a specific chemical property, the isotropic polarizability. We omit some of the methods from the previous comparison because of their high computational cost, and add in the greedy Nyström method (Fine and Scheinberg) with optimal weights (8).

Fig. 1 in the attached PDF shows the results. RPCholesky improves over all other methods we tested by a factor of three for isotropic polarizability and has the lowest worst-case error.

We thank the reviewers for the suggestion to include such an example. We hope that this example is both simple enough to be widely appreciated, and demonstrates the flexibility and utility of RPCholesky kernel quadrature for high-dimensional data problems.

Fine, S., & Scheinberg, K. (2001). Efficient SVM training using low-rank kernel representations. Journal of Machine Learning Research, 2, 243–264.

Meta-review: This is a high-quality submission that makes a major contribution to kernel cubature and is presented excellently.  I am delighted to recommend its acceptance at NeurIPS!
Decision: Accept (spotlight)