Positively Weighted Kernel Quadrature via Subsampling
Authors: Satoshi Hayakawa
Harald Oberhauser
Terry Lyons
Venue: NeurIPS 2022 Accept (2022)
Keywords: kernel quadrature
recombination
reproducing kernel Hilbert space
Nyström approximation

Abstract: We study kernel quadrature rules with convex weights. Our approach combines the spectral properties of the kernel with recombination results about point measures. This results in effective algorithms that construct convex quadrature rules using only access to i.i.d. samples from the underlying measure and evaluation of the kernel and that result in a small worst-case error. In addition to our theoretical results and the benefits resulting from convex weights, our experiments indicate that this construction can compete with the optimal bounds in well-known examples.

Positively Weighted Kernel Quadrature via Subsampling
Satoshi Hayakawa, Harald Oberhauser, Terry Lyons Mathematical Institute, University of Oxford (hayakawa, oberhause,,
( (lyons]@maths ox. ac uk
Abstract
We study kernel quadrature rules with convex weights. Our approach combines the spectral properties of the kernel with recombination results about point measures. This results effective algorithms that construct convex quadrature rules using only access to i.i.d. samples from the underlying measure and evaluation of the kernel and that result in small worst-case error. In addition our theoretical results and the benefits resulting from convex weights, our experiments indicate that this construction can compete with the optimal bounds in well-known examples.
Introduction
The goal of numerical quadrature to provide, for given probability measure μ on a space X, set of points 2,...., E X and weights Wi, wi....; W', IR such that
Σ (x.) (()) du(x)
holds for a large class of functions f X IR. Kernel quadrature focuses on the case when the function class forms a reproducing kernel Hilbert space (RKHS). What makes kernel quadrature attractive, is that the kernel choice provides simple and flexible way to encode the regularity properties a function class. Exploiting such regularity properties is essential when the integration domain X' is high-dimensional or the function class is large. Additionally, the domain X' does not have be Euclidean, but can any topological space that carries a positive semi-definite kernel.
1...., n) C R x X' denote with प = Σ w;dz, the resulting measure on X. We refer to Q [resp. n6 as a quadrature [resp. quadrature measure], to the points Z1...., Fn as the support of Q [resp. J. The aim of kernel quadrature is to construct quadrature measures & that have small worst-case error
More formally, given set Q = {ww;, 2;)
1
- f(x) f()) dp(x)
wee(Q; (cc((iHe,sH) = sup
(2)
Il $1
where He denotes the RKHS associated with positive semi-definite kernel k. If the weights positive and sum up to one, w; > 0. Σn wy = 1, then we refer to Q as a convex quadrature rule.
Contribution. The primary contribution of this article is leverage recombination (a consequence of Carathéodory's 's Theorem) with spectral analysis of kermels to construct convex kernel quadrature rules and derive convergence rates. We also provide efficient algorithms that compute these quadrature rules; they only need access to i.i.d, samples from and the evaluation of the kernel k. See Table [ for comparision with other kernel quadrature constructions.
CCode: httpp:://gichub..com/ssttttt hww.eaaaaaaaaaaaaaaaaaaeeee http::/hhh......oooooooo///aa/aa
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

The table is written by using and Tn, which represent a sort of decay of the kernel with respect to μ. Typical regimes are on nnn (e.g. Sobolev) or on exp( -yn) (e.g. Gaussian) depending on the 'smoothness' of the kemee [e.g. 16 3 (see also Section B.3 and in such regimes (with 3 2 2 or y > 0), on or ( noa) provide faster rates than wee? 1/n of the usual Monte Carlo rate. For
more examples including multivariate Sobolev spaces, see Bach 3
Section 2.3].
Limitation. Our proposed methods are based on either Mercer or Nyström approximation. Though baseer methods result in strong theoretical bounds, they require the knowledge of Mercer decomposition like B. 8 which not available for general (k, μ)). Our Nyström- based methods apply to much more general situations and outperform existing methods in experiments, but the n n/VE term makes their theoretical bound far from competitive. Further study needed to bridge the gap between theory and empirical results.
Method
Bound of squared wee
Computational complexity
M E
Herding 10
1/n Not found Om: (nggng))
n (n:: global optimization)
SBQ 2
n (n²: global optimization) Unavailable
Leveraged
DPP 7
Tn+1
n (rejection sampling)
CVS
a On+1 (1/n² + 1/N) polylog(N)
Unavailable
KT++
14
15
56
N Nlog* N
Ours:
Mercer
In Tn + N
nN₂ + C(n, No. n. + + log(N/n) neN, + nl2 + C(n, Nop néN + nz + + log((//)
M. + empirical
Nyström
n n + nn+11 Vi
N. + empirical#
nL + 11+1 + +
Comparison on n-point kernel quadrature rules. We are omitting the O notation throughout the table. Note that the assumption under which the theoretical guarantee holds varies from method to method, and this table displays just representative bound derived in the cited references. Here are remarks on the notation. (1) 7m is the m-th eigenvalue of the integral operator I., and F.. = 2i=m Tse (2) The symbols in the first line respectively mean C: convex, M: not using the knowledge of Mercer decomposition, and E: not using the knowledge of expectations such fx k(æ, y) duly)) (3) The (m: global optimization) is indicating the cost of globally optimizing function whose evaluation costs (mm) (t) Mercer/Nyström are the algorithms based on random convex hulls, see Section 2.4 and Appendix D ()) M./N + empirical are the algorithms discussed in the main text.
Table 1:
Why Convex Weights? There are several reasons why convex weights are preferable: (i) Posi- tive Integral Operator: Kernel quadrature provides approximation of the integration operator + I(f) = f(x) dp(x). Hence, a natural requirement is to preserve basic properties of this operator and positive weights preserve the positivity of this operator. (ii) Uniform estimates and Robustness: In applications, the RKHS H4 may be mis-specified if quadrature rule with negative weights is applied to function & H4. the approximation error ①I can get arbitrary bad; in contrast, a simple estimate shows that convex weights give uniform bounds, see Appendix [B.4 (iii) Itera- tion: Consider the m-fold product of quadrature formulas for approximating μ j on n This is a common construction for a multidimensional quadrature formulas (e.g., for polynomials) from
one-dimensional formulas Gil numeries for stochastic differential equations [E12] In doing so,
working with a probability measure is strongly preferred, since otherwise the total variation of their m-fold product gets exponentially large increases (I| [| = (
Related Literature. Roughly speaking, there have been two approaches to kernel-based quadrature formulas: kernel herding and random sampling. In kernel herding its variants, the points (x. ) are found iteratively, typically based on the Trank-Wolf gradient descent algorithm [10 25
In the random sampling approach, (x:)²= sampled and subsequently the weights optimized.
Generically, this results only in signed measure प but not probability measure. Bach [ and

---

Belhadji et al. 7 use the eigenvalues and the eigenfunctions of the integral operator K. : f H Sx k). du(y) to obtain a Mercer-type decomposition of k [59. Bach [3] then uses the eigenvalues and eigenfunctions of KC define an optimized measure from which the points (z:) are i.id, sampled. This achieves a near optimal rate, but the exact sampling from this measure is usually unavailable, although for special cases, can be done efficiently. In contrast, Belhadji et al. A proposes innindeeenddnt sampling based on the determinantal point process [DPP; 24] These two papers also treat the more general quadrature problem that includes a weight function E (2),,, i..,, approximating Jx (z) du(zz) for f E H4 C (2),,, which we do not discuss in this paper. Another recently introduced method is kernel thinning [1 15). which aims at efficient compression of empirical measures that can be obtained by sampling like our '+ empirical" methods. Its acceleation [56] makes it a competitive candidate in terms of compressing N ~~ n² points ('KT++' in Table
Finally, we emphasize that the kernel quadrature literature is vast, and the distinction between herding
and sampling is only a rough dichotomy, see e.g. 12 40 27
48
29
28
57 Beyond kernel
quadrature, our algorithms can also contribute to the density estimation approach in [4] which relies on recombination based on Fourier features although we do not pursue this further in this article.
Outline. Section 2 contains our main theoretical and methodological contribution. Section provides numerical experiments on common benchmarks. The Appendix contains several extensions of our main result, proofs, and further experiments and benchmarks.
2 Main Result
Assume we are given a se of n functions PI..... Pn-1 : X - E such that their linear combinations well approximate functions in Hy. Then our kernel quadrature problem reduces to the construction of an n-point discrete probability measure nΩn. Σ w;Or, such that
J (p
(((r)d((z) for every = n-
(3)
A simple way to approximately construct this Sº is to first, sample N >> n points, (U:) from y such that their empirical measure, jun N Σ oy. is good approximation to f in the sense that J (p; dun se J (i du for I,..., n 1, and secondly, apply a so-called recombination algorithm (Remark) I that takes as input (y:) and n functions (P,,...,PP: and outputs measure p.... Σ w.ddr, by selecting subset (x:) of the points (2:), and giving them weights (w:) such that o0 is a probability measure that satisfies the equation B with fu replaced by uN.
The challenging parts of this approach are (i) to construct functions PI...... P that approximately span the RKHS 74 for small n; (ii) to arrive at good quantitative bounds despite the (probabilistic) sampling error resulting from the use of the empirical measure NN. and the function approximation error via PP...... Pa-1- To address ii
we look for functions such that
k( Fo ko v) = Σ (xp,()))p:))
(4)
with some c; 2 0. Two classic ways to do this are the Mercer and Nyström approximations. The remaining. item (i) is our main contribution. Theorem I shows that the worst-case error, (5) is controlled by the sum of two terms: the first term stems from the kernel approximation 4 the second term stems from the sample error.
Theorem 1. Let μ be a Borel probability measure on X and k a positive semi- memi-definite kernel on X such that Sxx k(x, xr) du((x) < co. Further, let n be a positive integer and assume kn is positive semi-definite kernel on X such that
1. k- - kn is a positive semi-defiitte kernel on X, and 2. dim Hes < There exists a function KQuad such that if DN is set of N i.i.d. samples from μ., then Qn KQuad(D Dx) is a random n-point convex quadrature that satisfies
EDy [we He, 1)) S { (M(1,,)) (x)x +
203.
(5)
N
Thee number stems from Carathéodory's theorem, Remark 1 and leads to an n point quadrature rule.
3

---

where C,.u == Sx k(x, x)
(((,x))((()) du(x) - xfxxx k(z,v v) dp(x) du(y).
Moreover, the support of Qn is subset of DN and given functions Pi...., Pn-1 E L
(n)
with Hig C spauly!....., Pn-1), Qn = (Quad(Dx)) can be computed with Algorithm in OnNN= + n² Io)( computational steps.
The function KQuad is deterministic but since DN is random, the resulting quadrature rule Qn is random, hence also the resulting worst case error wee(Q: M4. and the expectation in 5 denotes the expectation over the N samples in Dy. The theoretical part of Theorem [T ffolows from more general results that we present and prove in the Appendix: Theorem//provee the inequality, essentially by comparing HA with Hoo: Theorem g proves the existence. The algorithmic part of Theorem | is discussed in Section 2. below. Theorem covers our two main examples for the construction of ko, resp. the choice of P1,..... P 1. and for which the error estimate gets quite explicit: the Mercer approximation, see Section 2.2 and the Nyström approximation, see Section 2.3 The former requires some knowledge about the spectrum of the kernel which is, however, known for many popular kernels; the latter works full generality but yields worse theoretical guarantees for the convergence rate. Finally, we emphasize that N and in Theorem can be chosen independently and we will see that from a computational point the choice N ~ n² is preferable in which case 5 is faster rate than Monte Carlo, also Table !
2.1 Algorithm
Algorithm 1 Kernel Quadrature with Convex Weights via Recombination KQuad
Input: A positive semi-defnitte kemel ls on X, a probability measure μ on X, integers N 2 n 2 1, another kernel ko functions ff,,..., P on X with Hs, C span(fi..... (ean..)......11)) and a set Dx of N i.i.d. samples from n..
Output: A set Qn == {(w,, 2;) Li.,,...,m n C R x X with w; 2 0, Σi Wi =
1: Apply a Recombination Algorithm (Remark 1 with ap = (4), ...., Pn- ding))) T to the empirical measure N EpED, by to obtain points (āi,.... [nt1] C DN and weights v (0+1)..,,++1)" 2 that satisfy 1Tv and (b)j N ΣN (ā.)). where (z)
[vb(z))..... (2)++ E (n+1)),,
2: Apply SVD with the matrix A = [Pt- E (=+(+ with go = 1 to find a nonzero vector
u E [R*+1 such that Au = and ding (2) u 0
3: Compute the smallest 2 0 such that au 2 0 and U/ auj for some j
4: Return (w:) + (2% - aru)) kEE and (x:);= + (**)ket, where I = {1...................... n + 1} {o}
Suppose we are given ko and ......, pn-1 E L (n) with Hk. C spanf (a.nff/1......... :..........1111.. pn-1) 1), and also N independent samples from μ denoted by D2 -yi....,,yy, Theorem in the Appendix shows that if we construct a convex quadrature Qn = (w), <;) satisfying
N
N
Σ (x;((x;)
=
1 N
Σ 1=1
((v)),
Σ wrki, (2x))
S
Σ
k1,diag(Vh) (V+),
(6)
N
1=1
where s0 (191, Pn and ks diidi (I) k(x, ko(x, x)), it satisfies the bound problem, we can use the so-called recombination algorithms:
5 For this
Remark 1 (Recombination). Given d functions (called test test functions) and probability measure supported on N > d points, there exists a probability measure supported on a subset of d points that gives the same mean to these d- functions. This follows from Carathéodory's theorem and is known as recombination. Efficient deterministic [38 43 67] as well as randomized [TT) algorithms exist to compute the new probability measure supported on d points; e.g. deterministic algorithmss perform the recombination, step D in O (C + d log(N/d)) time, where C.p is the cost of computing all the test functions at one sample. If each function evaluation is in constant time, Cup = (dd).
Let us briefly provide the intuition behind the deterministic recombination algorithms. We can solve the problem of "reducing (weighted) 2d points points in Rd while keeping the barycenter" by using linear programming or a variant of it. If we apply this to 2d points each given by barycenter

---

of approximately N 2d points, we can reduce the original problem of size N to problem of size d N By repeating this procedure log2( / times obtain the desired measure.
Although the recombination introduced here only treats the equality constraints in O we can satisfy
the remaining constraints just with points by modifying it. This is done in Algorithm I which works as follows: First, via recombination, find an (n+ popin convex quadrature Ra+1 (%) 31) n that exactly integrates functions .i...... ................................
1n di,diag with regard the empirical measure N Σ Bye- Second, reduce one point, find a direction (-u in the algorithm) in the space of weights on (5: )++1 that does not change the integrals of 99....., Pn-1 p and the constant function 1, and does not increase the integral of k didig Finally, move the weight from the above direction until an entry becomes zero, at t aru. Such whose entries sum up to one. Now we have convex weight vector with at most n nonzero entries, so it outputs the desired quadrature satisfying 6
0 exists,
must have positive entry since it is nonzero vector
2.2 Mercer Approximation
In this section and Section 2.3 we assume that & has a pointwise convergent Mercer decomposition
k(x,, y)
Im= CmCm (2)) (v) with d1 2 02 2 2 0 and (em)m= [²²(x) being orthonor-
mal [59] If we let IC be the integral operator (u) + (1)) given by f J fx k(, (v) du(y),
then (am: (mn)mn are the eigenpairs of this operator.
The first choice of the approximate kernel ko is just the trucation of Mercer decomposition.
Corollary 2. Theorem applied with ko( y) = Em=1 omem (z)em (m)m( (y) yields a random convex quadvature rule Qn such that
EDx wwee(Qn: Ho. 2 S Σ Fm
20k,u
N
m=R
Proof. It suffices to prove the result under the assumption Jx k(x, x) du(x)) Em on 00, as otherwise the right-hand side of 7 is infinity.
For hn k ko. we have that kn (xr, v) = E= om (z)em (y) and it is the inner product of +.....) .). and so positive semi-definite. Thus k and Ko satisfies the assumption of Theorem and fx hy (I, x) dp(z) Em om applied to
(xx) = (ymmem (x)))m)))) (r)) mm and ((y) in (2 +1,
yields the desired inequality.
2.3 Nyström Approximation
Although the Nyström method 668 1 34] is primarily used for approximating large Gram matrix by a low rank matrix, can also used for directly approximating the kernel function itself and this is how we use it. Given a set of l points Z =24)_= C X, the vanilla Nyström approximation of k(x,,) y) given by
k(x,8) v)
~ (P2k( x), P2k(: v)) =: k2 (x.
(8)
where P2 14 + H4 is a projection operator onto span [k((, 2:) In matrix notation, we have
k(21, v)
1² v) = k(x, Z)W (2),k((..))))) = [b(2,, 21), k(1, W+
(9)
k(ze, 4)
where W = (k(2,, z=)) is the Gram matrix for Z and W denotes its Moore-Penrose inverse. As k.7 is an Ê-dimensional kernel, there exists an (E + 1)-point quadrature formula that exactly integrates functions in Hzz. For quadrature formula, exactly integrating
We discuss the equivalence between x and (9) in Appendix B.5
the functions in Hyz is indeed equivalent to exactly integrating k(z;,) -) for all 1 < < É, as long as the Gram matrix k Z, Z) is nonsingular. Proposition I in the Appendix provides bound for the associated worst case error. From this viewpoint, the Nyström approximation offers a natural set of test functions.
5

---

The Nyström method has a further generalization with a low-rank approximation of k(Z, Z). Con- cretely. by letting W, be the best rank-s approximation of W k A(Z, Z) (given by eigendecomposi- tion), we define the following s-dimensional kernel:
1 (z, v) = k(z, Z)W," , k(Z, y).
(10)
Let W UAU be the eigendecomposition of W, where U [u1..... uel E Rexe is real
orthogonal matrix and A = diag(A)....,.,)) with x 2 ... 2 de 2 0. Then, Xe > we have
k₂ (x, v) = Σ
1 (w) SooTk(,,0 ((2 k(z,)))
(11)
i=1
So we can use functions u, k(Z, (i as test functions, which is chosen from a larger span(k(2; [k2 1)1-- Although closer to the original usage of the Nystöm method
1....,,))
dimensional space span
to obtain u, -(. as an approximation of i-th eigenfunction of the integral operator K. with Z appropriately chosen with respect to μ4, we have adopted an explanation suitable for the machine
learning literature 13 34
The following is a continuous analogue of Kumar et al.
34 Theorem 2] showing the effectiveness of the Nyström method. See also Jin et al. 26 for an analysis specific the case f.
Theorem 3. Let s < i be positive integers and 8 > O. Let Z be an E-point independent sample from g. If we define the integral operator K²² (p) + (²( by f + Sx kg (-. (y) duly), then we have, with probability at least 1 - 8, in terms of the operator norm,
|IK? K]ll S a++1 + 2supbEx k(x,x)
(1 (1+ 22 log 1
(12)
Ve
The proof is given in Appendix C.5 By using this estimate, we obtain the following guarantee for the random convex quadrature given by Algorithm aand the Nyström approximation.
Corollary 4. Let D be N N-point independent sample from μ and let Z be an l-point independent sam- ple from μ. Theorem ] applied with the Nyström approximation ko = knz yields an random n-point convex quadrature rute Qn such that, with probability at least 1 5 and kmax := suprex k(a, x.),
16(n - km
1 1 2 8
20R. N
EDD [weelQ Qu: Hk: 11)3 2 < 8 nan +
Σ m²n
Fm
+
Ve
Proof. From 2 (z. y) (x.) E ((,TT(Z)))) (u, ( k(2,))) v)) (ignore the terms with X; = 0 if necessary). and it is thus positive semi-definite. If we define Pz : Ha - Hs as the projection operator onto the orthogonal complement of span (ki( 2) } then, from ( we also have k(x, y) - h² (a, y) = (P (() k(,,),(),)))))) P k(),v)()))) He= so k - k² is also positive semi-definite. In particular, k 2 (k: 12) + (n22 1 is positive semi-definite. Also, it suffices to prove the result when _m On < co,
now apply Theorem
For hi == b/ 1- -1, we prove the inequality Sx ki (2. du((x) Σm= (em; (AC K.2 (em))m))p
(n - 1) K3 ||++ +mzn Om (see 31 in Appendix D.2 for details), and the desired inequality follows by combining Theorem and Theorem
3 (i.e., and )
Remark 2. Algorithm D with the Nystöm approximation can be decomposed into two parts: (a) Nystöm approximation by truncated singular value decomposition (SVD) (the first n 1 eigenvectors from an E-point sample), (b) Recombination from an N N-point empirical measure. The complexity of (a) is O (nl2), and it can also be approximated by randomized SVD in O + E2 log n) (20). The cost of part (b) is OnEN + nē' (og(N)))), where neN stems from the evaluation of ki,digg for all N sampling points. If we do not impose the inequality constraint regarding di,ii which still works well in practice, the cost of part (b) becomes O + n² (llllooo//) by using the trick
I EM U y:) - Ei=, k(Z,y.), y.), where Un-1 = [ui...., Un un-1] E ( is truncation of the matrix that appears in the Nyström approxiamtion 1o117 So the overall complexity is O(nEN + nl2 + n² log(N/no while an approximate algorithm (randomized SVD, without the inequality constraint) runs in OeN + l log n + ne log(N/n))).
N

---

2.4 Kernel Quadrature Using Expectations of Test Functions
Algorithm and the bound 5 can be generally applicable once we obtain low-rank approximation ku as we have seen in Section 2.2 and 2.3 However, since by construction we start by reducing the empirical measure given by NN. it mevitable have the (21/NN) term in the error estimate and performance. We can avoid this limitation by exploiting additional knowledge of expectations.
Let kn and h be positive definite kernels with k Ko k1. Let up
(f1 ((1:..... - be the vector of test functions that spans Hko- When we know the expectations of them, i.e., Jx (0(x) du(x)), we can actually construct a convex quadrature Qn = x;), satisfying
wapp(a;) = -fx (p(x) ()())))(())) du(x),
Wiski Wiski (Tx, an) S J k (1.
Σ
(13)
with positive probability by an algorithm based on random convex hulls (Appendix D
Algorithm
For this Qu. we have the following theoretical guarantee (see Theorem c in Appendix B Theorem 5. If convex quadrature On satisfies the condition
3 then have
wce (Ini Hk: n1)² < f
kn (z, dp(x).
ku is given the Mercer/Nyström approximations, we immediately have the following guarantees; they correspond to Mercer and Nyström in Table [1
See also Theorem 9 and 11 for details.
If ko y) = Em= Om Cmem (r) (v) is given by the Mercer approximation, we have
wce( Qni (oc
mm
for a convex quadrature On satisfying
13
Let Ko k-2 be given by the Nyström approximation 10 with Z being an 4-point
independent sample from je (with E > n). Then, for a convex quadrature Qn satisfying (3)
with probability least 1 - (with respect to Z) and kmax := suPzex k(a,,2), we have
wee Hk. 12)2 I no Σ 0m _om]
8(n 1)kmax
[+ / 2/1og
1
m>n
Numerical Experiments
In this section, we compare our methods with several existing methods. In all the experiments, we used the setting where we can compute fx k(x, y) dp(y) for E and [fxxx k(x, y) du(xx) du(y) since then we can evaluate the worst-case error of quadrature formulas explicitly. Indeed, if
quadrature formula Qn is given by points X = =(x;)² (2x) and weights ur = (wx) then we have
wee H4) = w 1(X, X) - 2 ww k(X. v)] + Ey. [k(y. v')l
(14)
for independent y.
11 under Sx / //( du(x) oc, which is well-known formula for An essential remark shown in Huszár and Duvenaud 25] is that the 49 with covariance kernel & given observation points (xx) (automatically)
the worst-case error 19
58
Bayesian quadrature estimates the integral as == (xz(();) with (wr) minimizing the above expression. Once given points (a:))= and additional knowledge of expectations,
can compute the optimal weights (w,) by solving convex quadratic programming (CQP), either without any restrictions or with the condition that (w:)r=1 1 is convex. Although the former can be solved by matrix inversion, we have used the optimizer Gurob for both CQPs avoid numerical instability. For the recombination part, we have modified the Python library by Cosentino et al. ||I implementing the algorithm of [61]
Our theoretical bounds are close to optimal in classic examples and we see that the algorithm even outperforms the theory in practice especially Section 3.T We also execute a measure reduction of large discrete measure in terms of Gaussian RKHS and our methods shows a fast convergence rate in two ML datasets in Section 3.2
Version 9.1.2, https: //www. gurobi com/
All done MacBook Pro, CPU: 2.4 GHz Quad-Cor Intel Core i5, RAM: 8 GB 2133 MHz LPDDR3.
7

---

0.-
logro(wce)²
-19-
log1o(wce)
-
N. + emp
---
N. emp
emp
-6-
-22
M. cmp
M. emp
M.
M.
-25-
iid Bayes
iid Bayes
Uniform Grid
-10-
Uniform Grld
Thinning
Thinning
-35-
Thin +
-12
Thin +
9.6
0.8
10
12
14
16
1.8
2.0
0.6
08
10
12
14
166
18
2.8
log1on
log1on
(a) d 1,==1
(b) d = 1.r
1092o(wcej'
0.0
Pogio(wcel)
N. 1 emp
N. emp
opt
Carlo
-1.0-
lld Bayes
lid Bayes
Halton opt
Halton 1 opt Thinning
Thinning
Thin opt
Thin opt
1.4
1.8
log1on
log10f
(d) 3,
Figure 1: Periodic Sobolev spaces with kernel k.dd; The average of log (wee(Qni Hy, 4) over 20 trials is plotted for each method of obtaining Qn The shaded regions show their standard deviation. The worst computational time per one trial was 57 seconds of Thin opt in (d, r, = (3, 3, 128). where Thinning was 56 seconds and N. emp [+ opt] was 22 seconds.
-1
10g1o (WCe)2
10g1o(wce)>
-4
+ emp
N. emp
N. emp
N. emp opt
Rd Bayes
idd Bayes
Herding
Herding
-7
Herd opt
Herd opt Thinning
Thinning
opt
1.
10g10n
log10f7
(a) 3D Road Network data
(b) Power Plant
Figure 2: Measure reduction in Gaussian RKHS with two ML datasets: The average of log (wce(Qu,, Qn: H&s (1)2) over 20 trials is plotted for each method of obtaining Qn- The shaded regions show their standard deviation. The worst computational time per one trial 14 seconds of Thinning [ opt] in Power Plant data with n = 128, where N. emp [ opt] was 6.3 seconds.

---

3.1 Periodic Sobolev Spaces with Uniform Measure
For positive integer r, consider the Sobolev space of functions on [0. 1| endowed with the norm II = So f(z) (2T)) fo (x)) dæ, where f and its derivatives f()),.. .... f(r) are periodic
(i.e, f(0) = f(1) and so forth). This function space can be identifies as the RKHS of the kernel
(x, 1+ (2x)2 B2r(1x -yl)
(2r)!
for z, y E 10, 1]. where B2, is the 2r-th Bernoulli polynomial 66 3..
If we let μ be the uniform mea- sure on [0, 1]. the normalized eigenfunctions (of the integral operator) are 1, G. v (2rm 2/2 sinmm for m ., and the corresponding eigenvalues are and m -2r (both for Cm and 9m). Although the rectangle formula f + n n- Ei=1 f(i/n) (a.k.a. Uniform 69 4vi in the sense of worst-case error, this RKHS is commonly used for testing the efficiency of general kernel quadrature methods [3] 7 28
and 5m (-)
Grid below) is known to be optimal for this kernel
We also consider its multivariate extension on [0, 1] i.e. the RKHS given by the product kernel k:d y) = II (x;, yk) for x =
= (y1...., ,9a) E [o, 1]d.
We carried out the experiment for rd,r) (1, (1, (2, 1), (3,3). For each (d, we compared the following algorithms for nt-oiint quadrature rules with n E {4, 16, 32, 64, 128).
N. + emp, N. + emp + opt: We used the functions (1
given by the Nys- as test functions (f...... :P in Algorithm
tröm approximation (Il) with The set was given as an ((==)(nnpponn) independent sample from μ.. We used N n samples from μ. In '+ opt' we additionally optimized the convex weights using
M. + emp, M. emp + opt (d 1): We used the first functions of the sequence of eigenfune- tions 1, C, 51, C22 $2,... as test functions .1.. in Algorithm We used N = n² samples from J. In '+ opt' additionally optimized the convex weights using
Monte Carlo, iid Bayes: With n-point independent sample (I:) from μı, we used uniform weights 1/n in Monte Carlo and the weights optimized using (14 in iid Bayes.
Uniform Grid (d = 1): We used the rectangle formula f + n E f(t/n). This is known to
be optimal [not just up to constant, but exactly: 69 47 and thus equivalent to the Bayesian quadrature on the uniform grid, i.e., the weights are already optimized.
Halton, Halton + opt (d > 2): For an n-point sequence given by the Halton sequence with Owen
scrambling 21 500 the uniform weights w; = 1/n is adopted in Halton and the weights are additionally optimized using [4] in Halton + opt.
Thinning, Thin + opt: Given an N-point independent sample (y: (w:) with N = n22 from μr, an
n-point subset (x;) taken from KT++ algorithm (kernel thinning [14 15 combined with Compress++ algorithm (56] with the oversampling parameter g = min 4, log, n). im- plemented with GoodPoints package: https://ggt is adopted in Thinning. In "+ opt' we additionally optimized the convex weights using (14).
com/ mi croooft /goo//oooonnsss
The results are given in Figure I In d = 1, the optimal rate given by Uniform Grid is known to be (n 2r). As the uniform sampling is equal to the optimized distribution of Bach 3 in this case, iid Bayes also achieves this rate up to log factors. Although our theoretical guarantee for M. + emp is Onn--2r + N11) with N = n² (Corollary 2). in the case (d,r) (1, 1), we can observe that in the experiment is better than iid Bayes and close to the optimal error of Uniform Grid, but slightly worse than Thinning. Moreover, N. + emp, which does not use the information of spectral decomposition, is remarkably almost as accurate as M. emp in = 1. Furthermore, if we additionally use the knowledge of expectations. which iid Bayes is already doing, M./. + emp + opt become surprisingly accurate even with N = n²? They are worse than Thinn + opt but well outperform it when r 3. Nonlineality in the graph of these methods when (d, r, n) = (1, 128) should be from numerical accuracy of the CQP solver (see also Section E.
when r
The accuracy of N. emp + opt becomes more remarkable in multivariate cases. It behaves almost the same as Halton opt in = 2 and clearly beats it in = Also, the sudden jump of our methods around 30 in (d, r) (3,3) seems to be caused by the jump of eigenvalues. Indeed, for the integral operator given by A3 with uniform measure, the eigenspace of the largest eigenvalue

---

1 is of dimension 27, and the next largest eigenvalue is 1/64. Again in the latter case, N. + emp + opt outperforms Thin opt, and these results suggest that our method works better when there is strong spectral decay, as is explicitly incorporated in our algorithm.
Note also that we can compare Figure T with Belhadji al. 7 Figure 1] which includes some other methods such as DPPs, herding and sequential Bayesian quadrature, as we did experiments under almost the same setting. In particular, in the case rd, = (1,3) where the eigenvalue decay fast, we see that our method substantially outperforms the sequential Bayesian quadrature.
3.2 Measure Reduction in Machine Learning Datasets
We used two datasets from UCI Machine Learning Repository (https://archiin ics. uci eduu ml/ /l/datasets/). We set μ as the equally weighted measure over (a subset of) the data points X =x()),,. ..... X (d 3, 5, respectively). where each entry centered and normalized. We considered the Gaussian kemeel exp(-1æ y/l /222)) whose hyperparameter À is determined by median heuristics [7,, and compared the performance N. emp, N. emp + opt (with é 10n, N ==n²), Monte Carlo, iid Bayes, Thinning, Thin opt. We also added Herding, an equally weighted greedy algorithm with global optimization [TC] and weight optimization Herd + opt within conver quadrature given by [4 We conducted the experiment for n E {4, 16, 32, 64, 128}.
The first 3D Road Network Data Set [3I. The original dataset is 3-dimensional real vectors at 434874 points. To be able to compute the worst-case eгror (4 efficiently to evaluate each kernel quadrature, we used a random subset X' of size 43487 434874/10] (fixed throughout the experiment) and defined as the uniform measure on it. We determined À with the median heuristic by using a random subset of X' with size 10000 and used the same X' and À throughout the experiment. The second is Combined Cycle Power Plant Data Set 5-dimensional real vectors 9568 points. We set the whole data as X' and defined μ as the uniform measure on it. We determined À with median heuristics by using the whole X..
32 63 The original dataset is
Figure 2 shows the results. We can observe that in both experiments N. emp opt successfully exploits the fast spectral decay of Gaussian kernel and significantly outperforms other methods. Also, even without using the knowledge of any expectations, N. + emp (and Thinning) show a decent convergence rate comparable to Herding or iid Bayes, which actually use the additional information. See also the end of Section E.2 for the plot of wee(Q H4. p²) for another set of empirical data μ².
Concluding Remarks
We leveraged classical measure reduction tool, recombination, with spectral properties of kernels to construct kernel quadrature rules with positive weights. The resulting algorithms show strong benchmark performance despite their restriction convex weights. Our method has also recently been applied to Bayesian inference problems
Although our method is applicable to fairly general situations, the usage performance can be limited when it difficult or inefficient to directly sample from the target measure μ.. Hence, an interesting follow up questions, is how one could replace the i.i.d. samples with smarter sampling (DPP, importance sampling, etc) before the recombination is carried out. Further, our theoretical results do not fully explain the empirical superiority; especially the 1/vē term does not match the experiments and it is challenging future research question to reduce this theoretical gap. Nevertheless, we believe our method is the first generally applicable algorithm with a guarantee from the spectral decay.
Acknowledgments and Disclosure of Funding
The authors would like to thank Chris Oates and Toni Karvonen for helpful remarks and discussions. The authors are also grateful to anonymous reviewers for detailed and constructive discussions that improved the paper. Harald Oberhauser and Terry Lyons are supported by the DataSig Program (FP//022347//1, the Alan Turing Institute (FP//510129/1] the Oxford-Man Institute, and the Hong Kong Innovation and Technology Commission (InnoHK Project CIMDA).
10

---

References
[1] M. Adachi, Hayakawa, M. Jørgensen, H. Oberhauser, and M. A. Osborne. Fast Bayesian inference with batch Bayesian quadrature via kernel recombination. In Advances in Neural Information Processing Systems, 2022. doi: 10448550/arXx..2206.47777
[2] A. Anastasiou, A. Barp, F.-. Briol, B. Ebner, R. E. Gaunt, F. Ghaderinezhad, J. Gorham, A. Gretton, Ley, Q. Liu, L. Mackey, J. Oates, Reinert, and Swan. Stein's method meets statisties: A review of some recent developments. arXiv preprint arXiv:2105.03481 2021.
[3] Bach. On the equivalence between kernel quadrature rules and random feature expansions. The Journal of Machine Learning Research, 18(1):714-751, 2017.
[4] Bach, S. Lacoste-Julien, and G. Obozinski. On the equivalence between herding and conditional gradient algorithms. In International Conference on Machine Learning, pages 1355-1362 1362, 2012.
[5] C. Bayer and J. Teichmann. The proof of Tchakaloff's theorem. Proceedings of the American Mathematical Society,
2006.
[6] A. Belhadji. An analysis of Ermakov-Zolotukhin quadrature using kernels. In Advances in Neural Information Processing Systems, volume 34, 2021.
[7] A. Belhadji, R. Bardenet, and P. Chainais. Kernel quadrature with DPPs. In Advances in Neural Information Processing Systems, volume 32, pages 12907-12917, 2019.
[8] A. Belhadji, R. Bardenet, and P. Chainais. Kernel interpolation with continuous volume sampling. International Conference Machine Learning, pages 725-735. PMLR, 2020.
[9] F.X. Briol, C. Oates, M. Girolami, and M. A. Oshorne. Frank Wolfe Bayesian quadrature: Probabilistic integration with theoretical guarantees. Advances in Neural Information Processing Systems, 28:1 111166111777 1170, 2015.
[10] Y. Chen, M. Welling, and A. Smola. Super-samples from kernel herding. In Conference on Uncertainty in Artificial Intelligence, pages 109-116, 2010.
[11] Cosentino, H. Oberhauser, and A. Abate. A randomized algorithm to reduce the support of discrete measures. In Advances in Neural Information Processing Systems, volume 33, pages 15100-15110, 10, 2020.
[12] S. De Marchi, R. Schaback, and H. Wendland. Near-optimal data tata-indeeendeet point locations for radial basis function interpolation. Advances in Computational Mathematics, 23(3):317-330, 2005.
[13] P. Drineas, M. W. Mahoney, and N. Cristianini. On the Nyström method for approximating Gram matrix for improved kernel-based learning. The Journal of Machine Learning Research, 6(1 6(12):2153-2175, 2005.
[14] R. Dwivedi and L. Mackey. Kernel thinning. In Conference Learning Theory, pages 1753-1753. PMLR, 2021.
[15] R. Dwivedi and L. Mackey. Generalized kernel thinning. In International Conference on Learning Representations, 2022.
[16] G. E. Fasshauer and M. J. McCourt. Stable evaluation of Gaussian radial basis function interpolants. SIAM Journal on Scientific Computing, 34())A737-A7722, 2012. URL https //doi. 1137/110824799
[17] D. Garreau, W. Jitkrittum, and M. Kanagawa. Large sample analysis of the median heuristic. arXiv preprint arXiv: 1707.00760766 2017.
[18] J.-y. Gotoh, A. Takeda, and K. Tono. DC formulations and algorithms for sparse optimization problems. Mathematical Programming, 169(1):141-17, 2018.

---

[19] A. Gretton, K. Borgwardt, M. Rasch, B. Schölkopf, and A. Smola. A kernel method for the poo-sample-problere Advances neural information processing systems. 19, 2006.
[20] N. Halko, P.-G. Martinsson, and A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217-288, 2011.
[21] J. H. Halton. On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. Numerische Mathematik, 2(1):84-90, 1960.
[22] S. Hayakawa. Monte Carlo cubature construction. Japan Journal of Industrial and Applied Mathematics, 38:561-577, 2021.
[23] S. Hayakawa, T. Lyons, and H. Oberhauser. Estimating the probability that a given vector is in the convex hull of a random sample. arXiv preprint vrXiv:2101.04250 2021.
[24] J. B. Hough, M. Krishnapur, Y. Peres, B. Virág, al. Determinantal processes and independence. Probability surveys, 3:206-229, 2006.
[25] Huszár and D. Duvenaud. Opeimally herding is Bayesian quadrature. In Conference on Uncertainty in Artificial Intelligence, pages 377-386, 2012.
[26] R. Jin, T. Yang, M. Mahdavi, Y.-F. Li, and Z.-H. Zhou. Improved bounds for the Nyström method with application to kernel classification. IEEE Transactions on Information Theory, 59 (0)::999-6699, 2013.
[27] M. Kanagawa, B. K. Sriperumbudur, and K. Fukumizu. Convergence guarantees for kemel- based quadrature rules in misspecified settings. Advances in Neural Information Processing Systems, 29: 3304,,344222 2016.
[28] M. Kanagawa, B. K. Sriperumbudur, and K. Fukumizu. Convergence analysis of determin- istic besnel-based quadrature rules in misspecified settings. Foundations of Computational Mathematics, 155, 2020.
[29] Karvonen, C. J. Oates, and S. Särkkä. A Bayes- cubature method. In Advances in Neural Information Processing Systems, volume 31, pages 5886-5897, 2018.
[30] T. Karvonen, C. Oates, and M. Girolami. Integration in reproducing kernel hilbert spaces of Gaussian kernels. Mathematics of Computation, 90(331):2009-22) 2233, 2021.
[31] M. Kaul, B. Yang, and C. S. Jensen. Building accurate 3d spatial networks to enable next generation intelligent transportation systems. In 2013 IEEE 14th International Conference Mobile Data Management, volume 1, pages 137-146. IEEE, 2013.
[32] H. Kaya, P. Tüfekci, and F. S. Gürgen. Local and global learning methods for predicting power of combined gas & steam turbine. In Proceedings of the International Conference on Emerging Trends in Computer and Electronics Engineering, pages 13-18, 2012.
[33] V. Koltchinskii and E. Giné. Random matrix approximation of spectra of integral operators. Bernoulli, 6(1):113-167, 2000.
[34] S. Kumar, M. Mohri, and A. Talwalkar. Sampling methods for the Nyström method, The Journal of Machine Learning Research, 13(1):911-1066, 2012.
[35] A. Kyrillidis, S. Becker, V. Cevher, and C. Koch. Sparse projections onto the simplex. In International Conference on Machine Learning, pages 235-243. PMLR, 2013.
[36] Larkin. Optimal approximation in Hilbert spaces with reproducing kernel functions. Mathe- matics of Computation, 24(1 (4(112):911992) 1970.
[37] P. Li, S. S. Rangapuram, and M. Slawski. Methods for sparse and low-rank recovery under simplex constraints. Statistica Sinica, 30(2):557-577, 2020.
[38] C. Litterer and Lyons. High order recombination and an application to cubature on Wiener space. The Annals of Applied Probability,
2012.
12

---

[39] F. Liu, X. Huang, Y. Chen, and J. A. Suykens. Random features for kernel approximation: A survey on algorithms, theory, and beyond. arXiv preprint arXiv:2004.11154. 2020.
[40] Q. Liu and Lee. Black-box importance sampling. In Artificial Intelligence and Statistics, pages 952-961. PMLR, 2017.
[41] J. Lu, G. Cheng, and H. Liu. Nonparametric heterogeneity testing for massive data. arXiv preprint arXiv: 0016001002200 2016.
[42] T. Lyons and N. Victoir. Cubature on Wiener space. Proceedings of the Royal Society of London Series A, 460:169-198, 2004.
[43] A. Maalouf, I. Jubran, and D. Feldman. Fast and accurate least-mean-squaes solvers. In
H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché Buc. E. Fox, and R. Gameet, editors, Advances in Neural Information Processing Systems, volume 32, pages 8305-8316, 2019.
[44] H. Q. Minh. Some properties of Gaussian reproducing kernel Hilbert spaces and their impli- cations for function approximation and learning theory. Constructive Approximation, 32(2): 307-338, 2010.
[45] H. Q. Minh, P. Niyogi, and Y. Yao. Mercer's theorem, feature maps, and smoothing. In International Conference Computational Learning Theory, pages 154-168. Springer, 2006.
[46] K. Muandet, K. Fukumizu, B. Sriperumbudur, and B. Schölkopf. Kernel mean embedding of distributions: review and beyond. Foundations and Trends in Machine Learning, 10(1-2): 1-141, 2017.
[47] E. Novak. Deterministic and stochastic error bounds in numerical analysis. Springer, 1988.
[48] C. Oates, M. Girolami, and N. Chopin. Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 79:695-715 -788,,7711711 2017.
[49] A. O' O'Hagan. Bayes Hermite quadrature. Journal of Statistical Planning and Inference, 29(3): 245-260, 1991.
[50] A. B. Owen. A randomized Halton algorithm in arXiv preprint arXiv:1706.02800 2017.
[51] V. Pan. On the complexity of a pivot step of the revised simplex algorithm. Computers & Mathematics with Applications, 11(11):1127 1140, 1985.
[52] M. Pilanci, L. El Ghaoui, and V. Chandrasekaran. Recovery of sparse probability measures via convex programming. In Advances in Neural Information Processing Systems, volume 25, pages 2420-2428, 2012.
[53] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, volume 20, pages 1177-1184, 2007.
[54] A. Sard. Best approximate integration formulas; best approximation formulas. American Journal of Mathematics, 71(1):80-91, 1949.
[55] R. Shamir. The efficiency of the simplex method: A survey. Management Science, 33(3): 301-334, 1987.
[56] A. Shetty, R. Dwivedi, and L. Mackey, Distribution compression in line-r time. In International Conference on Learning Representations, 2022.
[57] L. F. South, T. Karvonen, C. Nemeth, M. Girolami, C. Oates, et al. Semi-exact control functionals from sard's method. arXiv preprint arXiv:2002.00033, 2020.
[58] B. K. Sriperumbudur, A. Gretton, K. Fukumizu, B. Schölkopf, and G. R. Lanckriet. Hilbert space embeddings and metrics on probability measures. The Journal of Machine Learning Research, 1 1::117-1561,, 2010.
[59] Steinwart and C. Scovel. Mercer's theorem on general domains: the interaction between measures, kernels, and RKHSs. Constructive Approximation, 35(3):363-417) 2012.
13

---

[60] V. Tchakaloff. Formules de cubature mécanique à coefficients non négatifs. Bulletin des
Sciences Mathématiques, 81:123-134, 1957.
[61] M. Tchernychova. Carathéodory cubature measures. PhD thesis, University of Oxford, 2015.
[62] A. Tompkins and Ramos. Fourier feature approximations for periodic kernels in time-series modelling. In Proceedings of the AAAI Conference Artificial Intelligence, volume 32, 2018.
[63] Tüfekci. Prediction of full load electrical power output of base load operated combined cycle power plant using machine learning methods. International Journal of Electrical Power Energy Systems, 60:126-1444 2014.
[64] Tumne,, J. Liu, and P. Rigollet. statistical perspective on coreset density estimation. In International Conference on Artificial Intelligence and Statistics, pages 2512-2520. PMLR. 2021.
[65] U. Wagner and E. Welzl. continuous analogue of the upper bound theorem. Discrete & Computational Geometry, 26(2):205-215 2001.
[66] G. Wahba. Spline Models for Observational Data. Society for Industrial and Applied Mathe- matics, 1990.
[67] J. G. Wendel. A problem in geometric probability. Mathematica Scandinavica, (1 11(1):109-111 1963.
[68] C. Williams and M. Seeger. Using the Nyström method speed up kernel machines. In Advances in Neural Information Processing Systems, volume 13, pages 661-667, 2000.
[69] A. A. Zhensykbaev. Monosplines of minimal norm and the best quadrature formulae. Russian Mathematical Surveys, 36(4): 1981.
14

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Section
(c) Did you discuss any potential negative societal impacts of your work? [No] This is an improvement on a specific numerical method and believe it does not have forseeable negative social impact.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section and Appendix.
(b) Did you include complete proofs of all theoretical results? [Yes] See Section 2 and Appendix.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [Yes] Included the supplementary material as Zip file.
(b) Did you specify all the training details (c.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 3
(c) Did you report error bars (e.g.. with respect the random seed after running experi- ments multiple times)? [Yes] See Figures in Section g
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes] See the footnote 3 in page for the resource. We also included the computational time in each Figure (for the largest case) well the supplementary material (for all cases).
4. If you using existing assets (e.g., code, data, models) curating/releasing assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] See Section and the README file in the supplementary material.
(b) Did you mention the license of the assets? [Yes] See the README file the supple- mentary material.
(c) Did you include any assets either in the supplemental material URL? [Yes] See the README file in the supplementary material.
(d) Did you discuss whether and how consent was obtained from people whose data you're rou're using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
15


Reviews:
[{"review_id": "uHM8oC9Eb_m", "reviewer": "Reviewer_8cVE", "summary": "The authors investigate the properties of kernel quadrature rules with convex weights. The authors proposed new algorithms that take as input i.i.d. samples, a suitable kernel, and provide a small integration error in the worst-case over the associated reproducing kernel Hilbert space. Two strategies are proposed based on kernel approximations via Mercer decompositions and Nystrom approximations along with recombination algorithms for point measures.\n\n", "questions": "\n\nI have a few major concerns:\n\n1. For their Nystrom strategies (which are generically feasible for implementation), it is unclear if the mentioned bound can provide a better than Monte Carlo rate since that needs n / sqrt(l) < 1/ N which in turn means l > n^2 * N --- which is necessarily larger than N. On the other hand, the Mercer strategy in this work has a key limitation as it requires the knowledge of eigenfunctions, which is unknown for almost any generic (k, P) pair (other than a few exceptional pairs). Can the authors provide a detailed discussion about this aspect---as to the achieved rate in terms of the total number of queries---and provide concrete examples where their methods provide an improvement over prior strategies? In particular, further unpacking of Table 1 is needed, and it is unclear to me if the proposed strategy (E.g., see Table 2 https://arxiv.org/pdf/2105.05842.pdf].)\n\n2.  This work is missing relevant references on recent new developments on Kernel thinning (KT)---which would have a tick on all C, M, E, in Table 1 of this work. Moreover, KT provides an error of order 1/n^2 + 1/N (in terms of the notation of this work) up to logarithmic factors for the WCE with n point output, and N point input (typically referred to as MMD in the literature). See [Table 2 https://arxiv.org/pdf/2105.05842.pdf], [Table 3 https://arxiv.org/pdf/2110.01593.pdf] for error rates, and [Example 6 https://arxiv.org/abs/2111.07941] for the new variant that provides near-linear runtime (N log^3 N), which is significantly faster than the strategies mentioned. \n\n3. [Sec 1.2 https://arxiv.org/pdf/2105.05842.pdf] (KT paper) would also be a useful reference for some other works like P-greedy, black box importance sampling, etc. which are not referenced here and are relevant for the work.\n\n4. For experiments, I believe KT should be added. Moreover, like the strategies in this work, even KT and herding output should be combined with 'opt' to provide a fairer comparison to their best-performing variant (N.+emp.opt]. \n\n5. At least one of the algorithms from Remark 1 should be discussed with some brief intuition for the readers to become aware of the high-level ideas used.\n\n6. What is the probability in Thm 3 / Corr 4 taken over (is it the draw of Z)?\n\nMinor comments:\n\n- l 38/39: A reference is needed.\n", "limitations": "See questions.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nStrengths of the work include (a) perhaps a novel combination of two concepts: Spectral properties of kernels, and recombination results about point measured, and (b) the ability to simultaneously analyze strategies involving Mercer decomposition, and Nystrom approximations.\n\nLimitations are (a) lack of clear discussion on when the proposed methods improve over prior strategies, while a summary is provided in Table 1, it is unclear if the methods indeed provide an improvement, (b) missing prior works, and (c) lack of stronger baselines that can be used in experiments. (see questions for more details),\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "lXO6x6hAijt", "reviewer": "Reviewer_vQx7", "summary": "This problem focuses on the quadrature problem, more specifically on bounding the maximum mean discrepancy between a quadrature approximation mu_Q= sum w_i delta_{x_i} and a target distribution mu, where the (xi) is a n-subset of an available bigger sample DN of i.i.d. samples of mu.\nThey focus on convex rules i.e. when the weights are positive and sum to 1.\n\nThe main idea is to consider a basis of n functions phi_1.. phi_n approximating RKHS functions, then sample N points (y1…yN), and select a n-subset (x1..xn) reweighted by (w1..wn). There are thus two sources of error (1) approximating k/the RKHS by a n-finite dimensional one, e.g. through Nystrom or Mercer approximations, (2) the empirical measure supported on y1…yN.\n\nThey provide several results, in particular th 1, which becomes more explicit in Cor 2 through Mercer approximation and Cor 4 through Nystrom approximation. The authors evaluate empirically their algorithm on small dimensional datasets (with either Nystrom or Mercer approximations) and compare it to the Monte Carlo rate , uniform grid and iid Bayes (weights minimizing equation 14) and demonstrate it achieves faster rates in practice than Monte Carlo.\n", "questions": "- Is it possible to clarify the dependence in number of samples in the quadrature n and dimension d in the bounds of Cor 1 and 2?\n- can you report the slopes for various dimensions and infer a dependence dimension experimentally?\n\n\nMinor comments\nl55: “hybrid approaches” -> hybrid between what and what?\nl69: a subset (xi)_{i=1}^n (n instead of N)\nl71: (iii) in practice, constructing N>> n samples can be a challenge as well (eg MCMC methods are expensive)\nl138: Nystom ", "limitations": "The limitations of the theoretical and empirical results could be discussed much more extensively, see comments above. ", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\nThe paper is quite well written, and well referenced regarding kernel quadrature rules. It tackles an important problem, quadrature rules, which enables to approximate integrals accurately, with applications in Bayesian inference for instance.  The proposed algorithm has complexity O( nN + n 3 log(N/n) \u0001) which is competitive with herding in high dimensions, where the latter become intractable because of the global optimization subroutine. \n\nWeaknesses\n\nMaybe I am missing something, but I think the bound in Theorem 1 does not advocate for faster rates than Monte Carlo.\n\nIn theorem 1, the expected squared MMD (where the expectation is taken wrt the available sample DN) is bounded by a constant c_1(n) (later on it will depend on n) + c_2/N. \n\nClassical bounds on (squared) MMD (e.g. “A kernel two sample test”, Gretton et al.) between the empirical measure on DN (without quadrature weights) and the target mu is already of order c_2/N. Hence, from this bound it does not seem that the quadrature provides a better approximation of the target distribution mu, since if I understood well, Qn  (of size n \\le N) can be of the same size than DN (size N) a priori.\n\nThe bound of Theorem 1 starts to be interesting if n is always much smaller than N and the term int(k-k0)dmu is much smaller than 1/n.\nThe first condition is typically satisfied (the authors claim l98 they will take N ~ n^2 ) later, so c_2/N becomes of order c_2/n^2, which is faster than Monte Carlo (c/n) . However, for the second term c_1(n)  it is less clear. Corollary 2 (Mercer approximation) yields a term c_1(n), indeed decreasing with n,  but it is not clear how fast c_1(n) decreases, faster than 1/n? Same question for Corollary 4 (Nystrom approximation) which is not very explicit in n. \n\nif it is not faster than 1/n, then the quadrature rule (supported on n points) does not have a better upper bound currently than any empirical measure supported on n points of DN.\n\nHowever, the experimental results show that the steep is lower (better) than the Monte Carlo rate, which encourages the method. \nIn the experiments, it would be interesting to report the steeps and understand the dependence of the rate with respect to n but also with respect to the dimension.\n \nFor these reasons, I tend to reject the paper, but I am ready to revise my opinion if the authors can clarify their results.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "Eaxwvp3MNXx", "reviewer": "Reviewer_bzNP", "summary": "The paper addresses a problem of kernel quadrature rules and proposes a new quadrature with positive weights that leverages recombination algorithms (stemming from Caratheodory’s Theorem) and either Mercer’s or Nystrom approximation of the kernel. Based on each version, there is also a derivation of corresponding convergence rates. The proposed methods are empirically verified on a classical setup for testing efficiency of general kernel quadrature methods.\n", "questions": "What are the consequences of treating a less general quadrature problem without weighting function $g$?\n\nCould authors elaborate more on the need of convex quadrature? And how the positiveness of weights affects practical performance?\n\nHow the choice of recombination algorithm affects the performance?\n", "limitations": "The potential negative societal impact is not applicable to this submission.\n\nThe paper addresses a less general case of kernel quadrature, i.e. there is no weights function $g$ in the problem formulation (also stated in the related literature paragraph).\n\nThe method implies an ability to sample datapoints from target measure for recombination technique, however this could be problematic (this limitation is also stated by authors).\n", "rating": 6, "confidence": 3, "soundness": 1, "presentation": 1, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Originality: \n\nThe paper addresses an established task of kernel quadratures. The method is based on well-known techniques of recombination and spectral kernel approximation combined in a unexplored way to obtain convex quadrature method. The paper emphasizes the major difference of their method being th convex weights of the quadrature points.\n\nQuality: \n\nI’m not entirely convinced by the very brief motivation for the convex weights, perhaps it would benefit the paper to elaborate the need for the convex quadrature in a more comprehensible way.\n\nIt seems that extensive supplementary material provides a substantial theoretical support for the proposed method, unfortunately, I was not able to check that in detail.\n\nThe Experimental Section lacks explicit comparison with the other practical methods apart from the worst and best case baselines on a known setting with Sobolev and Korobov spaces, e.g. there is no DPP-based method in the Figures.\n\nIt is also not immediately clear whether convexity of the weights has a crucial role in any performance advantages, thus the motivation for convex quadrature rules lack empirical support (at least in the experiments conducted and the conclusions drawn from their results).\n\nThe only method that is consistently bringing advantages across all Figures is +opt variation which assumes knowledge of expectations andadditional optimization over the weights, however Table 1 does not include this information which brings a little bit of a confusion. I wish that was a little bit more clearly stated in the paper.\n\nClarity:\n\nThere is a general issue with exposition, at least for someone less acquainted with the area. Abstract as well as Related Work would greatly benefit from a more extended exposition. As mentioned above, the motivation for the convex weights is also very short and could be made more clear with more explanations. Experimental Section doesn’t make it any clearer in this regard too, perhaps adding an ablation study coud offer more transparent understanding. The paper doesn’t discuss how their derived bounds compare against other methods, even under different assumptions, perhaps there should be matching setups allowing comparison.\n\nSignificance:\n\nIt is definitely an interesting method, and perhaps will raise a good amount of interest in the kernel community. However, I’m not sure whether the paper makes a good work at demonstrating the advantages of the proposed approach over existing methods. It first argues that convex quadrature is essential, but later fails to confirm this empirically, unfortunately, and it seems positive weights do not play a major role in theoretical analysis neither.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "CKquFpF2P-O", "reviewer": "Reviewer_pDdL", "summary": "This article studies a family of quadrature rules suited for functions that belongs to an RKHS. The proposed construction is based on a recombination algorithm that takes a discrete measure $\\nu_{N}$, that approximates the initial measure $\\mu$ (for example, a Monte Carlo approximation),  and outputs another discrete  measure $\\mu_{n}$. The weights of the quadrature are obtained by enforcing that the quadrature rule is exact on n functions $\\phi_{1}, \\dots, \\phi_{n}$ that are taken to be equal to $\\phi_{i} := k_{0}(x_{i},.)$ where $k_{0}$ is a ‘low rank’ kernel. \n\nThe contributions of the paper may be summarised as follows:\n- A generic result (Theorem 1) that gives the upper bound of the worst-case error (on the unit ball of the RKHS) for the proposed algorithm for an arbitrary 'low rank' kernel $k_0$\n- The instantiation of Theorem 1 to the case when $\\nu_{N}$ is the Monte Carlo approximation of $\\mu$ and $k_0$ is obtained from the Mercer decomposition of $k$\n- The instantiation of Theorem 1 to the case when $\\nu_{N}$ is the Monte Carlo approximation of $\\mu$ and $k_0$ is obtained through the Nyström approximation\n- Several numerical simulations that illustrate the theoretical rates. \n\nThe article is well written and the proven results were widely discussed and compared to the existing results in the literature. \n\n", "questions": "\nQuestions and suggestions:\n\n* In Section 3, it is not clear what you mean by the opt version: the optimization of (15) is done in the simplex or $R^{N}$?\n\n* How do you explain the fact that in Corollary 2 the boundedness assumption (condition (a) in Theorem 8) is not required?\n\n* It would be nice to add the N-th eigenvalue as a benchmark in the graphs\n\n* Theorem 8 and Theorem 10 were mentioned implicitly in Table 1. In my opinion, I believe that they deserve to get an explicit mention in the main paper. ", "limitations": "-", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\n\nStrengths:\n* The empirical versions of the algorithms are flexible and may be used in domains where the eigenfunctions of the integration operator are not tractable.\n* The proposed quadrature rules are convex, which is an important property in misspecified settings\n* The theoretical analysis is insightful and may be used for other applications \n\n\nWeaknesses:\n \n* The empirical versions of the quadrature (when $\\nu_N$ is Monte Carlo approximation of $\\mu$) are very useful in practice, yet they come with weak convergence rates: the second term in the r.h.s. of the bound (5) in Theorem 1 is $\\mathcal{O}(1/N)$, which is typically a slow rate in the kernel-based quadrature literature. ", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept