Fast Summation of Radial Kernels via QMC Slicing
Authors: Johannes Hertrich
Tim Jahn
Michael Quellmalz
Venue: ICLR 2025 Poster (2025)
Keywords: fast kernel summation
slicing
quasi-Monte Carlo
non-equispaced fast Fourier transforms
random Fourier features
Primary area: other topics in machine learning (i.e., none of the above)

Abstract: The fast computation of large kernel sums is a challenging task, which arises as a subproblem in any kernel method. We approach the problem by slicing,  which relies on random projections to one-dimensional subspaces and fast Fourier summation. We prove bounds for the slicing error and propose a quasi-Monte Carlo (QMC) approach for selecting the projections based on spherical quadrature rules. Numerical examples demonstrate that our QMC-slicing approach significantly outperforms existing methods like (QMC-)random Fourier features, orthogonal Fourier features or non-QMC slicing  on standard test datasets.

Published as a conference paper at ICLR 2025
FAST SUMMATION OF RADIAL KERNELS VIA QMC SLICING
Johannes Hertrich', Tim Jahn", Michael Quellmalz"
1 Université Paris Dauphine-PSL. and UCL, johannes www.heihhhuuupphhhaeo...c pas.ipphhneese..............p..m psl.eu Technische Universität Berlin, {jahn, 2 tu-berlin.
ABSTRACT
The fast computation of large kernel sums is a challenging task, which arises as subproblem in any kernel method. We approach the problem by slicing, which relies on random projections to one-dimensional subspaces and fast Fourier sum- mation. We prove bounds for the slicing error and propose a quasi-Monte Carlo (QMC) approach for selecting the projections based on spherical quadrature rules. Numerical examples demonstrate that our QMC-slicing approach significantly outperforms existing methods like (QMC-)raadom Fourier features, orthogonal Fourier features or non-QMC slicing on standard test datasets.
INTRODUCTION
We consider fast algorithms for computing the kernel sums
8m W., K K Bm), for all 1,..., M,
(1)
where In, 3m E RRd and w., E TR for = Nm
M and K: Rd Rd is radial
kernel, i.e., K K(x, y) = F[[ - yll) for the Euclidean norm |-I|| and some F: R20 RR. This sum- mation problem appears in most kernel methods including kernel density estimation (Parzen 1962 Rosenblatt 1956). classification via support vector machines (Steinwart Christmann [2008), di- mensionality reduction with kernelized principal component analysis (Schölkopf & Smola 2202]
Shawe-Taylor & Cristianini] 2004), distance measures on the space of probability measures like maximum mean discrepancies or the energy distance (Gretton et al. 2006 Székely 2002), corre- sponding gradient flows (Arbel et al., [2099] Galashov et al., [004]] Hagemann et al. 2024] [Hertrich] al., 22224] Kolouri et al. (0022), and methods for Bayesian inference like Stein variational gradient descent (Liu Wang 2016). Computing (T) exactly for all m I,..,, M has complexity O(MN), which can be restricting if M and N are large.
In low dimensions, there is a rich literature on fast approximation algorithms, we include a (non- exhaustive) list in the "related work" section. One particular approach is the fast Fourier summation (Kunis et al.] 2006) Potts et al. 2004), which approximates the kernel by a truncated Fourier series and reformulates 1 using the fast Fourier transform non-equispaced data (Beylkin) 1995 Dutt & Rokhlin 1993). We provide a short overview in Appendix [1] This kind of methods usually provides computational complexity of O(M N + Nr. log( Nv. )). where Nt. is the number of relevant Fourier coefficients, and admits very fast error rates for N2 + co (even exponential if the kernel is sufficiently smooth). However, the number N₃ of relevant Fourier coefficients grows exponentially with the dimension d, such that this method is computationally intractable for dimensions larger than four.
As a remedy for higher dimensions, Rahimi & Recht 2007) proposed random Fourier features (RFF). They represent a positive definite kernel via Bochner's theorem (Bochner) 1933) as the Fourier transform of non-negative measure. Sampling randomly from this measure at D N points (features) leads to an approximation algorithm with computational complexity ((D(N+M) independent of the dimension However, the error decays only with rate OL/VVD), which can be limiting if high accuracy is required. Moreover, RFF are limited positive definite kerness and

---

Published as a conference paper at ICLR 2025
do not apply to other kerells, like the negative distance kernel [(z, y) = -11a y/l, which has applications, e.g., within the energy distance (Székely) 20022 that is used for defining a distance on the space of probability measures.
A related approach is slicing (Hertrich 2024). which represents the kernel sum T as an expectation of one-dimensional kernel sums of the randomly projected data points with different kernel. Dis- cretizing the expectation by sampling at random projections, the kernel sums (1) can be approx- imated by P one-dimensional kernel sums, which can be computed efficiently, e-g., by fast Fourier summation. Similarly as for RFF, this leads to a complexity of O(P(N + M + Nn log Nn.). where the expected error can be bounded by ( (11/P) /P). For positive definite kernels, close relation be- tween RFF and slicing was established by Rux et al. (2025). see the short overview Appendix G One advantage of slicing is the applicability to kernels that are not positive definite.
A way to improve the O(1/VP rate is to replace the uniformly chosen directions with specific sequences of points. This yields so-called quasi-Monte Carlo (QMC) algorithms on the sphere, see (Brauchart et al.] 2014). Note that there also exist QMC approaches for RFF (Avron et al.) 2016). However, they depend on the restrictive assumption that the measure from Bochner's 's theorem decouples over the dimension, which is true for the Gauss and L Lapplae keapppaaaeeeeelll... kernel, but false for most other common kernels like the Laplace, Matérn or negative distance kernel.
Contributions Our contributions for fast kernel summation in Rd via slicing follows:
We derive bounds on the slicing error for various kernels in Theorem iirremmmiincininng all positive definite radial kernels. Furthermore, we exactly compute the variance for the negative distance kernel, the thin plate spline, the Laplace kernel and the Gauss kernel.
We exploit QMC sequences on the sphere in order to improve the error rate OLL/VP). To ensure the applicability of the QMC error bounds, we prove certain smoothness results for the function which maps a direction & E g to the corresponding oie-dimensional kernel in Theorem The improved eтror rates are outlined in Corollary 4
We conduct extensive numerical experiments standard test datasets for several kernels and different QMC sequences, and demonstrate that our QMC slicing approach with the proposed distance QMC designs significantly outperforms the non-QMC slicing method as well as (QMC-)RFF. While the advantage of QMC slicing is most significant in dimensions d < 100, also performs better in higher dimensions.
Outline In Section 2 first revisit the slicing approach in detail and present our improved error bounds. Afterwards, in Section 3 we consider quadrature and QMC sequences on the sphere and prove the applicability of the approaches for slicing. We present our numerical results in Section 4 and draw conclusions in Section []] Additional proofs, plots and evaluations are contained in the appendix. The code for the numerical examples provided onlind"
RELATED WORK
Now-Dimensional Kernel Summation Fast summation algorithms have been extensively stud- ied in the literature. They include fast kernel summations based on (non-) jequippced fast Fourier transforms (Greengard et al. 2022, Kunis et al. 2006 Potts et al., 2004). J. fast multipole methods Beatson & Newsam 1992, Greengard & Rokhlin 1987; Lee & Gray) 2008, Yarvin & Rokhlin 1999). tree-based methods (March et al. 2015a b) or H- and mosaic-skeleton matrices (Hackbusch)
1999 Minden al., 2017 TTytyshnikov 1996p. For the Gauss kernel, the fast Gauss transform was
proposed by Greengard & Strain (1991) and improved by Yang et al. (2003 2004). More general fast kernel transforms were considered by Ryan et al. (2022).
QMC and Quadrature on Spheres QMC designs on spheres were studied by Brauchart et al. 2014). Here, the quadrature points optimizing the worst-case error in certain Sobolev spaces given by spherical (-designs, which integrate all polynomials up to degree t on the sphere exactly
In literature, there exist versions of the Laplace kernel K(x, y) exp(-ol|a- yil1) and K(x, y) exp(-al|æ v/l), which differ used norm. Since our analysis focuses on radial kernels, we will only consider the second version in the rest of this paper. "available at wthttps:/
2

---

Published as a conference paper at ICLR 2025
The construction of spherical designs is highly non-trivial and intractable in high dimensions. For S and S several examples were computed
Delsarte al. 1977| Bannai & Bannai 2009).
numerically by Gräf & Potts (201T) and Womeesley 2018). Gräf et al. (2012) related quadrature rules on the sphere with halftoning problems.
Sliced Wasserstein Distance The idea of slicing is also used in optimal transport. A sliced Wasser- stein distance was proposed Raaai et al. (2122. In contrast to the kernel summation problem, the sliced and non-sliced Wasserstein distance do coincide and have different properties. QMC rules for the three-dimensional sliced Wasserstein distance were considered by Nguyen et al. (2024).
Random Fourier Features Random Fourier features (RFFs) were proposed by Rahimi Recht (2077) and were further analyzed in several papers Bach 2017): Hashemi et al. (2023): Li et al.
2021). To improve the error rates, Avron al. (001)) proposed quasi-Monte Carlo approach for RFFs under the restrictive assumption that the measure from Bochner's theorem decouples over the dimensions. This approach was refined by Huang et al. (2024) and Munkhoeva al. (2018). Yu al. (2016) proposed orthogonal random features. In a very recent preprint, Belhadji et al. (2024) derive an explicit quadrature rule in the Fourier space for efficient summations the Gauss kernel.
SLICING OF RADIAL KERNELS
Let K: Rd x Rd be radial kernel of the form ((x,3) = F([ - y/l) for some basis function
F: Rzu + Throughout this paper, we will assume that K has the form
K (z, 3) Ee-Uga- [K( 2²) (S. v)
for some one-dimensional radial kernel k: R x IR + I with basis function f: R20 + R. where we suppress the dependence of f and F onto the dimension d. By inserting the basis functions of the kernels, this corresponds to the slicing relation
F([1==] = E (f(
(2)
A pair (F, f) of basis functions in Lice (R>o) fulfills this relation if and only if F is the generalized Riemann-Liouvile fractional integral transform given by
[(ls)( - 22 4=3 ds,
F(1) (
(3)
for 2 < d CENN N, see (Hertrich) 2024 Prop 2 and Rubin 2003). In order to find the one-dimensional basis function f' for given F. we have to invert the transform This can be done explicitly if
i) F' is analytic on Rzo, i.e., there exists (aw) C IR such that F(x) = > anæ" for all
x 2 0, ii) F(I-- II): I + R is continuous, bounded and positive definite, i.e., for all N E N, all pair- wise distinet 2, Rd and all aj R for 1,..., it holds that ΣjA-1 aga F( Ix 0,
see (Hertrich) 2024 Thm and Rux et al. 2025 Cor 4.11). We include a list of pairs (F. fulfilling ( in Appendix A In particular, includes the basis functions of Gauss, Laplace and Matérn kernels. Moreover, 1 can be computed for other important choices that fulfill neither i) nor ii), e.g., the thin-plate spline and the generalized Riesz kernel.
2.1 FAST KERNEL SUMMATION VIA SLICING
In order to compute the kernel sums ① efficiently, we approximate F(] II) by projections and the one-dimensional basis function f. i.e. we aim to find directions EI...., Ep gd-1 such that
F([
()))),z))) for all E Rd.
Then, the kernel sums I read as
8m
Wr, K (In, Um)
Σ
wn wn F(22
small)
Σ
Wn (1(Ep Bm))). (5)
p=1
3

---

Published as a conference paper at ICLR 2025
For computing the one-dimensional sums E Wn, onnfllle En ym))) for all I,..., M, there exists algorithms with complexity O(M + N) or O + N log(M + NN) in literature. These include fast summations based on non-equispaced Fourier transforms (Kunis et al. Poott et al. 2004), fast multipole methods (Grrengardd & Rokhlin 1987) or sorting algorithms (Hertrich et al.
2024 In partieular, we can approximate the vector (81, 8as) via (5) with complexity of PP (M + N)).
2.2
ERROR BOUNDS FOR UNIFORMLY DISTRIBUTED SLICES
To bound the error of the slicing procedure from the previous subsection, we consider error estimates for the approximation in (4). To this end, we assume that the directions E Ep are iid samples from the uniform distribution on the sphere. Then, we exactly compute the variance
Valf(()) == Eg [(f(KE)))) - F(I
(6)
which bounds the mean squared error through the Bienaymé's identity as
27
Valf(()) P
B Epl
P
x()
x)1)
In particular, our results show relative error bounds of the negative distance kernel K(x,y)
-||0 y/ and the thin plate spline (except around || = 1), which do not depend on the di-
mension d. The proof is given in Appendix
Theorem 1. Let F: R2o - and f: R200 + IR fulfill the slicing relation 2
If F(I I]) is continuous and positive definite on IRd, then Vals((>)= F(0)2 F(
ii) For the generalized Riesz kernel F([|=[]) -||=||11 with we have
(r())))(())) + 442)2
+ ())(()((())))))
VTI(r + 1)
valf(z))
1 F(1
F( (7)
1(
iii) For the thin plate spline F([[=1)) = ||0|| ||>l|||og log([1=]]). we have
+ C2 + O log d)
)
Valf((z) = 2 + log(11=]
(1+ 3) F([=+))?
(og((
with 01 4.189 and 02 = 2.895 given (8).
iv) For F(11=1)) En an ||ll=|11" and d 2 3 odd, we have
(d (d-1)/2
Vall(()) = Σ
II
(1
1+ (1+ (2i 1)(n [m[++ + 2i - 13 -1
k(n k)
anan -k ||
particular, for the Laplace kernel ([[[=]]) exp(-allaan we have
Valf((x)= भলحلo (1 (2(a[1=:1) + 20/|2=11 + 1)F(((=))))))
and for the Gauss kernel F[[|=]])= = exp( ||x=| (/axll/((20))) /22=))) we have
d² Valf((x) = "ाm (1 ||||| | + +1) F(I==1)))
2
Some weaker error bounds for the Gauss and Riesz kernel were also shown in Hertrich (2024), see Appendix Dopendix(D.Ijeeee for a detailed comparison. In all cases, the variance is independent of the dimension d. The dependence on x differs between the kernels: For positive definite kernels, which are always bounded, we have an absolute error bound i) independent of ||æ/|| For the Riesz kernel, we have relative error bound in ii). For the thin plate spline, iii) behaves like relative bound for ||>]]| + oo and ||>||| + 0, but as constant around || 1, which is zero of F. For the Laplace and Gauss kernel, the dependence on ||ж||| changes drastically between ||||| + 00 and ||ж||| + O. In fact, (x) is monotonically increasing in ||æ||| with global upper bound 1/ 1/(4a), and converges quadratically in |||||| to zero for ||ж||| + O. For the case > 3, we conjecture a similar behavior, see Appendix D.2 for the discussion.
4

---

Published as a conference paper at ICLR 2025
3 QUASI-MONTE CARLO SLICING
For directions drawn independently from the uniform measure Uga-1 on the sphere S our ex- periments from the numerical part suggest that the rate O(1/VP) from Theorem is optimal. As remedy, we employ quadrature and quasi-Monte Carlo designs the sphere for improving these error rates. To this end, we first revisit the corresponding literature in Subsection Afterwards, we apply these results for our slicing method in Subsection 3.2
3.1 QUASI-MONTE CARLO METHODS ON THE SPHERE
...., SP) E ( for P E N. In the following, we aim to construct & such that the Let Le 10 , worst case error in a certain Sobolev space is asymptotically optimal. The definition of the Sobolev space Hd H²(Sa-1) is given in Appendix B
Definition 2. A sequence (E²)p with P - oa is called a sequence of QMC designs for Ha if there exists some c(8, d) independent of P such that the worst case error
P
d(s,d ps/(d-1) E (((P-o/(a-1)).
sup gE (S4-1) (ل 1151
[dd-1 gd g(E)dE-
g
(8)
P
p=1
was proven by Hesse 20066 that the enror rate O(P->/(4-1) is optimal, see also Brauchart al. 2014 For الية the existence of sequences of QMC designs ensured by so-called spherical designs. More preciscly, &" is called spherical t-design if the quadrature at these points integrates all polynomials of degree E N exactly, i.e., if it holds
1 P Σ
(E)DE
(E
for all polynomials 1 of degree < t.
g4
It can be shown that for any t 'three exists a spherical (-design with P = O(d-11) points, see Bondarenko al. (2013). By (Brauchart & Hesse 2007, Cor 3.6), such sequence of spherical t-designs is a sequence of QMC designs. see also (Brauchart et al. 2014 Sect 1) for summary.
Unfortunately, the construction of spherical t-designs in arbitrary dimension tractable. Instead, many QMC methods rely on cow-discrepancy point sets. It was shown in (Brauchart al. 2014, Thm 14) that a sequence & that minimizes the sum of powers of Euclidean distances
numerically in-
E(EP)
Σ p-g=1
E 2-120-d-11
In the numerics, we will consider 2' so that we get QMC design for (dd2( which we call the distance QMC design. Note that, up to a constant, 99) coincides with the maximum mean discrepancy with the Riesz kernel K K(x,
QMC design for H= for E 1-1 d.d),
2 12)
|æ 2- between the probability measures 1 ΣP= dep and Uga-1, which is easily transform QMC
also known energy distance Székely (2002). Furthermore,
sequence into an unbiased estimator in see Appendix H
3.2 SMOOTHNESS OF
BASIS FUNCTIONS
In order apply the above theorems for our approximation
2 we need to ensure that for any
E IRd the spherical function 9x $d-1 IR given by
(E) = f((5. 1)
(10)
sufficiently smooth on gd-1, For some specific examples, this is verified in the next theorem, whose proof is given in Appendix F For part ii), we explicitly compute the Sobolev norm of ga-
Note that, in the special case O. this relates to the variance (6) by the formula ||gz (200)
(Val[()) + 2( for (J, F) fulfilling (2)
5

---

Published as a conference paper at ICLR 2025
Theorem 3. Let B with + O. For the Gauss, Riesz and Matérn kernel, the following smoothness results hold true:
12 i) For Ftt) exp(
nn ) we have Ir E H for all 2
for
ii For (t)) ts with t 2 0 and -1, we have 9u E H= and only if's < r
ii For F(1) = Nojj 2 1) v 22 t). 0, we have 9r E HH(SS-1) if < + 2*
Note that the theorem also includes the Laplace kernel, which the Matérn kernel for v/ = Com- bining this theorem with the results from the previous subsection leads to improved error bounds for the Gauss and Matérn kernel in the following corollary. For the Riesz kernel, the last theorem can be seen as a negative result that the theory from the previous subsection is not applicable.
Corollary 4. Let d E N and >> L Then there exists a constant c(s, d) and sequence (EP) with P > oo such that for the Gauss and Matérn kernel with basis functions F(1) = exp( 20z 122
22-" and (()) T(m) 2 r) Kul (2/20 1) with > 2 respectively, it holds that
B
1 P 1 E +)))
c(s.d) P
sup
F([+I))
ZER4
p=1
For ( (d=1 d+1), such E" are given by minimizers of (9) 2
In Appendix J we derive the complexity of the biciingbbased kernel summation with QMC se- quence and show in Proposition [ that it superior to the random Fourier feature approach for the Gauss kernel. The summation methods are described in Appendix Mand Appendix G respectively.
4 NUMERICAL EXAMPLES
In the following, we evaluate the kernel approximation with QMC slicing for several QMC se- quences and compare our results with random Fourier feature-based (RFF. Rahimi Recht 2007) methods. We implement the comparison in Julia and Python, and provide the code online Addi- tionally, we provide a general PyTorch implementation for fast kernel summations via slicing and non-equispaced fast Fourier transform In Subsection 4.1 we describe the used QMC sequences, RFF-methods and kernels. Afterwards, in Subsection 4.2 we numerically evaluate the approxima-
tion error in (4. Finally, we apply approximation for fast kernel summations in Subsection
We include additional numerical examples in Appendix K
4.1
QMC SEQUENCES AND KERNELS
QMC Sequences Beside standard slicing where the projections E" are drawn iid from the uniform distribution on gd-1 we consider the following sequences. We would like to emphasize that for the first two of them not clear, whether they are QMC designs in the sense of Definition even though they are sometimes called QMC sequences in the literature.
Sobol Sequence: Two commonly used QMC sequences on [0, 1] are Sobol (Sobol' 1967 and Halton (Halton 1960) sequences. They can be transformed QMC sequences for the normal distribution by applying the inverse cumulative density function along each dimension the sequence, which was used for deriving QMC method for random Fourier features Avron et al. (011)). To obtain a potential QMC sequence on the sphere, Nguyen (tt al. al. (2024) proposed to project the QMC sequence for the multivariate normal distribution onto the sphere by the transformation { = 0/||0]|, see also Beltrán et al. (2023). It is not known whether this constitutes a QMC sequence the sense of Definition 2 To generate the original Sobol sequence [0, 1]d, use the implementation from SciPy (Virtanen ct al. 2020 in Python and the Sobol. j1 package in Julia. Our numerical experiments
available at http:///w/aaaaaaaaooaooaaaaaa/aa available httpa//////gttth..ooo////etttttr herpppeeethh_hrp
QMC Sac_epiimmMmaaninnnn

---

Published as a conference paper at ICLR 2025
suggest that the Sobol sequence lead to slightly better results than the Halton sequence. Therefore, we omit the Halton sequence in our comparison.
Even though this is technically not QMC sequence, we adapt the approach of orthogonal features (Yu et al. 2016) for slicing and generate directions & as follows: We generate [1 orthogonal matrices from the uniform distribution on (dd) (this can be done by taking the Q-factor of the QR decomposition applied on matrix with standard normally distributed entries). Together. these matrices have dl / columns from which we choose
Orthogonal:
e" to be the first P of those.
Distance: In Section 3.1 we considered the distance QMC design E for (S which is a minimizer of E(EP) = -
|l - 1. see ( In our application, we Ep --Ep
construct the distance QMC designs &º by minimizing the functional (e ) == E((5". -E)))-
have the additional symmetry that ()) = f(((2. D) Therefore,
We do this numerically with the Adam optimizer (Kingma & Ba 2015) and the PyKeops package (Charlier et al. 202T). which takes from a couple of seconds (for d = 3) up to hour (for d 200 and P 5000) on an NVIDIA RTX 4090 GPU. In Appendix [H, we show that if P S. d, the orthogonal points from above minimize Esym. so this approach differs only if P >
Spherical Design: For d = 3, several spherical I-designs on the $ were computed by Gräf & Potts (001 up to S 1000 and P 1002000 and are available online Spherical t-designs for S* were computed by Womersley| (2018) up to S 31 and P < 3642. Unfor- tunately, the computation in higher dimensions appears to be intractable such that we only use the spherical designs for d = 3.
Compared Methods
We compare our results with the following methods:
Random Fourier Features (RFF, Rahimi & Recht, 2007 : see Appendix G for descrip-
tion.
Orthogonal Random Features (ORF, Yu et al. 2016):
The directions of the RFF features
are chosen in the same way as explained above for the orthogonal slicing.
QMC- RMC-Random Fourier Features (Sobol RFF, Avron et al., 2016): For the Gauss kernel, we also compare with QMC random Fourier features, which are only applicable for kernels where the Fourier transform decouples a product over the dimensions. For the kernels from Table 2 this is only true for the Gauss kernel. As a QMC sequence in [0, 1]d, we choose the Sobol sequence.
Kernels We use the Gauss, Laplace, Matérn (with = 3 for p E (1) the negative distance kernel (Riesz kernel with 1) and the thin plate spline kernel, see Table 2 in the appendix for the pairs F). The parameters a, a and 3 are chosen by the median rule (see, e.g., Garreau et al, 2017 for overview and history). That is, we choose all considered input norms ||>||| of the basis functions F and y is some scaling factor which we set to y 1,2).
3
ym, where m is the median of
4.2 NUMERICAL EVALUATION OF THE SLICING ERROR
We examine the approximation error in 4 numerically. To this end, we draw a sample x from (0, 0.1f) and evaluate the absolute error [[[[=]]) - Epe f(((5p. x) We average this error over 50 realizations of & (whenever & is random) and 1000 samples of The average results for the scale factor 7 = 1 and dimension d E {3, 10, 50} with the Gauss, Laplace and Matérn kernel are given in Figure [ We observe that all methods despite the spherical designs for the Gauss kernel converge with rate ((P-r) for some r > 0. To estimate the rate r numerically, we fit regression line in the loglog plot. The resulting rates r are given in Table Further plots and tables are given in Appendix K. considdeeg the negative distance kernel, higher dimensions and smaller/larger kernel widths.
Overall, the distance QMC designs perform best in most examples, except when the (provably op- timal) spherical designs are applicable, which are only computable in < 4 and outperform the
ht httpss//wwee nser. chennit xx.. poot www.///rrwgggggeeeooeeoo /tg.e///aurrrrreeeeerarr wetere///ddexxxpee wttpe/////deeeeeeeeee

---

Published as a conference paper at ICLR 2025
Gouss ww 6000
Gouss nemel with of -0.405 50
20
13
26-
13
12--
AP
-2-
Oaf
1858g
Sant
Scod RP?
13-0
Sendtitiing stoing
Sand sloing
23
Obcoree
20
Drn deir
Sunde
Rumber Proletorn
Natem cennel with p- 0.15, pm3 3.5.d d
Matorn komd p- 0.31 u-
with
20-2
10
134
off
Sctelsting orn siing
dicing
dingng
Sani dicin
Scnce slcing
11
sprotce drdde
Dittacce sidgg
10-0
Dittance ding
Matern kemil with 0.15 uw 3
Matern kemd with 0.31 ow
Mater karnal with 0.70 1.5 2.5.0
30
20.
20.
stong
Off
Eittace
Ornulees
Ore siiee
Lasiace barna o 6.74,
Lasiace kerne 3225
Laplace with gm
30"1
20-1
EXT O4R
Sobol D10
sling
diana
Gdoccaal coiiu
Difaace dorig
1.0
Figure 1: Loglog plot of the approximation error |F(1=l)) Σ, f(KSp. I) for approximating
the function F by slicing 4 versus the number P of projections (or the number D = P of features for RFF and ORF) for different kernels and dimensions (left d 3. middle d = 10, right d 50). The results are averaged over 50 realizations of &" and 1000 realizations of z. The kernel parameters set by the median rule with scaling factor 7 1. We fit regression line in the loglog plot for each method to estimate the convergence rate r, see also Table T
Table 1: Estimated convergence rates for the different methods. We estimate the rate by fitting a regression line in the loglog plot. Then, we obtain the estimated convergence rate p-" for some faster convergence. The resulting values of
r > Consequently, larger values of correspond r' are given in the below tables, the best values are highlighted in bold. The kernel parameters are the same as in Figure (median rule with scaling factor 7 1).
kernel with medion sallng y
Madie kemel with - and mesiaa rule scaline
REF 3e5-b
Slicing based
REE hood
Sixing besed
Dimension
RFF
Sobul
ORF
Slicing
Subol
Disumes
Dineesion
RFV
ORF
Slkcng
Sobel
Om
Distncc
spherical dedgi
R150
0.30
0.96
2,10
0.51
0.59
3.99
0.96
0.50
0.28
0.50
1.36
0.50
0.50
1150
10500
3.76
0.67
0.00
0.72
0.70
0.79
50
0.55
0.57
0.50
0.67
0.64
0.71
Mirica kemel with
and mediaa cale
semell with median
Ri-i-tesed
Slicing Peaed
RITi-bused
Sticing Nowed
Dimcnason
RFF
ORF
sphoniicl dloign
sphorical
2' 3
1U.0
0.31
0.51
0.56
0.58
2.11
2.24
3
BU2
0..0
0.50
0.88
0.52
1.31
1.28
10
(8.00
01.00
0.50
0.70
(150)
0.89
d 1D
BUG0
0.50
0.50
0.G]
(1.00
0.65
50
1U50
0.54
0.50
0.61
0.00
0.66
50
0.59
0:50
0.58
0.59
8.60
8

---

Published as a conference paper at ICLR 2025
Gauss kernel on Letters (dlm 16)
Gousd kernd on MMST (d 20)
Couss kernel on Fas HNIST (d-30))
PalMbve
CMC (Sckeb apd
CMe (sclit and
Oet (saavs and
Sarieg
099 Mi sariig
Oee siring
10
12"
10-
30
12'
angle trraded CPU
Time drglettraude
Matern tu Latters (di-
(o- 2.5] MNIST 1d
(15 FashionMSST 30)
12-1
Restov DDre
Reaatee Dore
Frasler sicino
sicine Owe Fourie Striig
20
Pouter dine
qpel Stireg
qeet siring
10
CPU
Time drgle-traaded
dng CPU
Matem 1.5) Letters
Motern (v MNIST
Matem (= 1.5) FashionMMST
32-1
10+1
10
Recatree 4 Spear
Sp
cep
ORF
dee 20
OMC NI shcing
Qee to skring
10-1
10-1
Time single whroaded CPU
Time teo
Thme 50 on cnroodee CPU
Laplace kemel on Lotters (d)
Lapace rerd on MNIST 201
Laploce kernel on Fashion (d-30)
ABF
RFF CeF
RFF
Touter Slicin
30
shiiing CMC Four Stidng
Fouter
OMC Fourne Sliiig
10
CMC Fourle Sliing
22
30
10"
10
16
Figure 2: Loglog plot of the relative L' approximation error versus computation time for computing the kernel summations (I) with different kernels and methods. We the Letters dataset (M N = 20000 points), MNIST (reduced to dimension = 20 via PCA, M = N = 60000 points) and FashionMNIST (reduced to dimension d 30 via PCA, M N 60000 points). We each method 10 times. The shaded area indicates the standard deviation of the error. For Fourier slicing, use P 5 2* slices for 1, .... 10. In order to obtain similar computation times, we use 5 2 slices for RFF-10 slicing and D = 2/' features for RFF and ORF.
distance QMC designs for smooth kernels as they reach machine precision already for P = 250. In accordance with Corollary the benefits of QMC slicing are better for smooth kernels and low dimensions. But also for d 50, significant advantage of QMC slicing is visible. In particular, we often observe much faster convergence rates than the proven worst case error rates r z(d11) on H (Sd-1) with the distance QMC designs, Section 3 Furthermore, the slight advantage of the distance QMC designs versus the spherical designs for the Laplace kernel is because the former ones are chosen specifically for symmetric functions.
4.3 FAST KERNEL SUMMATION
Finally, we test our kernel approximation for computing the whole kernel sums from ①.. For com- puting the one- mie-dimensoonll kernel sums Σm=1 w, ( -ym))) we use the following methods combined with either random QMC points the sphere:

---

Published as a conference paper at ICLR 2025
(QMC) Sorting-Slicing: For the negative distance kernel, the sorting algorithm from Hertrich et al. (2024), see also (Hertrich 2024, Sec 3.2).
(QMC) Fourier-Slicing: For the Gauss, Matérn and Laplace kernel, we use fast Fourier summation based on the non-equispaced fast Fourier transform (NFFT) for the one- dimensional kernel summation. A general overview on NFFTs and fast Fourier summa- tion can be found in the text book (Plonka et al., 2023 Sec 7.5). Similar as (Hertrich) 2024 Sec 2.3), we do not evaluate the one-dimensional basis functions, but directly com- pute the Fourier transforms. We revisit the background on the one-dimensional fast Fourier summation and specify the used parameters in Appendix
For positive definite kernels (Gauss, Laplace, Matérn), we use one mne-dimensional random Fourier features with & features for the basis function f. For iid orthogonal slices and k = 1, this approach is related to RFF and ORF, as outlined in Appendix
(QMC) RFF-K Slicing:
For E",, we use the distance QMC designs, since we have seen in the previous subsection that it performs best among the QMC rules. Moreover, we use the randomization from Appendix E] for the QMC design obtain unbiased estimator. We evaluate the kernel sums on Letters dataset (d = 16, Slate 1991), MNIST (reduced to = 20 dimensions via PCA, LeCun et al. 1998) and FashionMNIST (reduced to d 30 dimensions via PCA, Xiao et al. 2017), where (x11 -..g N.N' and (y1....,y/M) constitute the whole dataset and the weights (W1, .... WN] are set to 1. In particular, we have M- N 20000 for the Letters dataset and M N 60000 for MNIST and Fashion- MNIST. Then, we approximate the vector s = ($1, aM)) from T and report the absolute error ||s Ssaellll We choose the kernel parameters by the median rule with scale factor ? = 1 based on 1000 example pairs (z. y). We benchmark the computation times on single thread of an AMD Ryzen Threadripper 7960X CPU and compare our results with RFF, ORF and (non-QMe Since the QMC designs & depend neither on the dataset nor on the kernel, we consider its con- struction not as a part of the computation time. For the Gauss kernel, we also compare with QMC (Sobol) RFF (Avron et al. 2016), which is not applicable for the other kernels. We use P= 5 2* slices for k = 1.... 10 in the slicing method. In order to obtain similar computation times, we use 5 2 slices for RFF-10 slicing and D 2P features for RFF and ORF
slicing.
We visualize the approximation error (including standard deviations) in Figure 2 for the Gauss, Laplace and Matérn kernel. The results for the negative distance kernel, for the thin plate spline ker- nel, for higher dimensional datasets (including the full MNIST and FashionMNIST with d = 784) and GPU comparison are included in Appendix K.2 In addition, we apply the fast kernel summa- tion to computing MMD gradient flows in Appendix K.3 We can see that QMC Fourier-slicing has a significantly smaller error than the other methods. Moreover, it has the smallest standard deviation of the error.
CONCLUSIONS
Summary and Outlook We proposed slicing approach to compute large kernel sums in O(P(N + M + Nn log Nn)) instead of the naïve M) arithmetic operations. In the case of iid directions, we proved eтror bounds with rate O(1vP)). To improve this rate, we proposed a QMC approach based on spherical quadrature rules. We demonstrated by numerical methods that our QMC slicing approach outperforms existing methods, where the advantage is most significant for dimensions d < 100. In the future, we want improve our theoretical analysis on QMC slic- ing in order to match the convergence rate from the numerical section. One possible approach for that could be to study worst case errors for symmetric functions on the sphere, since the mappings + (x(E) from Section 3.2 are always symmetric, see Appendix K.4 for details. From practical side, we want to apply the slicing approach in some actual applications.
Limitations In the numerical part, observe significantly better error rates than we can prove theoretically, see Appendix K.4 for a discussion. Moreover, the computation of the QMC directions can be very costly and depends strongly the chosen method. For the spherical designs, it is even intractable for high dimensions. Finally, the advantage of QMC slicing becomes smaller for higher dimensions, which is well-known effect for most QMC methods.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
JH acknowledges funding by the German Research Foundation (DFG) within the Walter Ben- jamin Programme with project number 530824055 and by the EPSRC programme grant The Mathemattics of Deep Learning with reference EP/V026259/1. TJ acknowledges funding by the Deutsche Forschungssemmeinschatt under Germany's s Excellence Strategy EXC-2046/1, Projekt-II 390685689 (The Berlin Mathematics Research Center MATH+). MQ acknowledges funding by the Genman Research Foundation (DFG) with reference STE 571/19-2 and project number 495365311, within the Austrian Science Fund (FWF) SFB 10. 55766666/ Tomography Across the Scales. For open access purposes, the authors have applied CC BY public copyright license to any author- accepted manuscript version arising from this submission. We gratefully thank Gabriele Steidl for fruitful conversations.
REFERENCES
Fabian Altekrüger, Johannes Hertrich, and Gabriele Steidl. Neural Wasserstein gradient flows for discrepancies with Riesz kernels. In International Conference Muchine Learning, pp. 664- 690. PMLR, 2023.
Michael Arbel, Anna Korba, Adil Salim, and Arthur Gretton. Maximum mean discrepancy gradient flow. Advances in Neural Information Processing Systems, volume 32, 2019.
Kendall Atkinson. An Introduction to Numerical Analysis. Wiley, 1991. ISBN 978-0-471-66489-9)
Kendall Atkinson and Weimin Han. Spherical Harmonics and Approximattins the Unit Sphere: An Introduction. Springer, Heidelberg, 2012. doi: 10. 110009978-3-662-2598833
Haim Avron, Vikas Sindhwani, Jiyan Yang, and Michael W. Mahoney. Quasi-Monte Carlo feature maps for shift-invariant kernels. Journal of Machine Learning Research, 17(120):1-38, 2016.
Francis Bach. On the equivalence between kernel quadrature rules and random feature expansions. The Journal of Machine Learning Research, 18(1):714-751, 2017.
Eicchi Bannai and Etsuko Bannai. A survey on spherical designs and algebraic combinatorics on spheres. European Journal of Combinatorics, 30(6):1392-144 1425, 2009. doi: e0.1016//..jc.2008. 11.007.
Richard Beatson and Garry Newsam. Fast evaluation of radial basis functions: I. Computers & Mathematics with Applications, 24(12):7-19, 1992.
Ayoub Belhadji, Qianyu Julie Zhu, and Youssef Marzouk. On the design of scalable, high-precision spherical-radial Fourier features. arXiv preprint arXiv:2408.132. 2024.
Carlos Beltrán, Damir Ferizović, and Pedro R López-Gómez. Measure-preee preserving mappings from the unit cube to some symmetric spaces. arXiv preprint 2303.00405, 2023.
Gregory Beylkin. On the fast Fourier transform of functions with singularities. Applied and Com- putational Harmonic Analysis, 2(4):363-381, 1995.
Salomon Bochner. Monotone Funktionen. Stieltjessche Integrale und harmonische Analyse. Math- ematische Annalen, 108(1):378-)) -100, 1933.
Andriy Bondarenko, Danylo Radchenko, and Maryna Viazovska. Optimal asymptotic bounds for spherical designs. Annals of Mathematics, 178(2):443 452,, 2013. doi: 10.4007/amals.2133177.. 2.2.
Johann Brauchart and Kerstin Hesse. Numerical integration over spheres of arbitrary dimension. Constructive Approximation, 25(1):41-71, 2007.
Johann Brauchart, Edward Saff, Ian Sloan, and Robert Womersley. QMC designs: optimal order quasi Monte Carlo integration schemes on the sphere. Mathematics of Computation, 83(290): 2821-2851, 2014.
11

---

Published as a conference paper at ICLR 2025
Benjamin Charlier, Jean Feydy. Joan Alexis Glaunès, François- David Collin, and Ghislain Durif. Kernel operations on the GPU, with autodiff, without memory overflows. Journal of Machine Learning Research, 22(74))1-6, 2021.
Zonghao Chen, Aratrika Mustafi, Pierre Glaser, Anna Korba, Arthur Gretton, and Bharath Sriperumbudur. (de)-regularized maximum mean discrepancy gradient flow. arXiv preprint arXiv:2409.14900 2024.
Feng Dai and Yuan Xu. Approximation Theory and Harmonic Analysis on Spheres and Balls. Springer, New York, 2013. doi: 10. 11100/9978---461466660
Philippe Delsarte, J. M. Goethals, and Johan Jacob Seidel. Spherical codes and designs. Geometriae Dedicata, 6(3):363-388 1977. ISSN 1572-9168. doi: 10.10077/f03188704.
Alok Dutt and Vladimir Rokhlin. Fast Fourier transforms for nonequispaced data. SIAM Journal on Scientific Computing, 14(6):1368-113) 1993.
Alexandre Galashov, Valentin de Bortoli, and Arthur Gretton. Deep MMD gradient flow without adversarial training. arXiv preprint TrXiv:2405.067800 2024.
Damien Garreau, Wittawat Jitkrittum, and Motonobu Kanagawa. Large sample analysis of the median heuristic. arXiv preprint arXiv:1707.0720 2017.
Izrail S Gradshteyn and Iosif M Ryzhik. Table of Integrals, Series, and Products. Academic Press New York, seventh edition, 2007.
Manuel Gräf and Daniel Potts. On the computation of spherical designs by optimization approach based on fast spherical Fourier transforms. Numerische Mathematik, 1 19(())699-7224 2011.
Manuel Gräf, Daniel Potts, and Gabriele Steidl. Quadrature rules, discrepancies and their relations to halftoning the torus and the sphere. SIAM Journal Scientific Computing, 3465)22700-2791, 2012.
Leslie Greengard and Vladimir Rokhlin. A fast algorithm for particle simulations. Computational Physics, 73(2):325-348, 1987.
Journal of
Leslie Greengard and John Strain. The fast Gauss transform. SIAM Journal on Scientific and Statistical Computing, 12(1):79-99) 94, 1991.
Philip Greengard. Manas Rachh, and Alex Barnett. Equispaced Fourier representations for efficient Gaussian process regression from billion data points. arXiv preprint arXiv:2210.10210, 2022.
Arthur Gretton, Karsten Borgwardt, Malte Rasch, Bernhard Schölkopf, and Alex Smola. A kernel method for the two-sample-problem. Advances in Neural Information Processing Systems, 19, 2006.
Wolfgang Hackbusch. A sparse matrix arithmetic based H-Matrices, Part I: Introduction to H- Matrices. Computing, 62(2):89-108, 1999.
Paul Hagemann, Johannes Hertrich, Fabian Altekrüger, Robert Beinert, Jannis Chemseddine, and Gabriele Steidl. Posterior sampling based on gradient flows of the MMD with negative distance kernel. In International Conference on Learning Representations, 2024.
John Halton. On the efficiency of certain quasi-random sequences of points in evaluating multi- dimensional integrals. Numerische Mathematik, 2:84-90, 90, 1960.
Abolfazl Hashemi, Hayden Schaeffer, Robert Shi, Ufuk Topcu, Giang Tran, and Rachel Ward. Gen- cralization bounds for sparse random feature expansions. Applied and Computational Harmonic Analysis, 62:310-330, 2023.
Johannes Hertrich. Fast kemel summation in high dimensions via slicing and Fourier transforms. arXiv preprint arXiv:2401.0822 2024.
12

---

Published as a conference paper at ICLR 2025
Johannes Hertrich, Christian Wald, Fabian Altekrüger, and Paul Hagemann. Generative sliced MMD flows with Riesz kernels. In International Conference on Learning Representations, 2024.
Kerstin Hesse. A lower bound for the worst-case cubature error on spheres of arbitrary dimension. Numerische Mathematik, 103:413-4333 2006.
Zhen Huang, Jiajin Sun, and Yian Huang. Quasi-Monte Carlo features for kernel approximation. In International Conference Machine Learning, 2024.
Jens Keiner, Stefan Kumis, and Daniel Potts. Using NFFT3 software library for various noneq- uispaced fast Fourier transforms. ACM Transactions on Mathematical Software, 36:Article 19, 1-30, 2009.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. In Proceedings of the ICLR '15, 2015.
Tobias Knopp. Marija Boberg, and Mirco Grosser. NFFT.jl: Generic and fast julia implementation of the nonequidistant fast Fourier transform. SIAM Journal on Scientific Computing, 45(3):C179- C205, 2023.
Soheil Kolouri, Kimia Nadjahi, Shahin Shahrampour, and Umut Şimşekli. Generalized sliced prob- ability metrics. In IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 4513 2022.
Stefan Kunis, Daniel Potts, and Gabriele Steidl. Fast Gauss transforms with complex parameters using NFFTs. Journal of Numerical Mathematics, 14(4):295, 2006.
Yann LeCun, Leon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, (6(11):2278-233) 1998.
Dongryeol Lee and Alexander Gray. Fast high-dimensional kernel summations using the Monte Carlo multipole method. D. Koller, Schuurmans, Y. Bengio, and Bottou (eds.), Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008.
Zhu Li, Fran-Francois Ton, Dino Oglic, and Dino Sejdinovic. Towards unified analysis of random Fourier features. The Journal of Machine Learning Research, 22(1):4887-4937 2021.
Jen Ning Lim, Juan Kuntz, Samuel Power, and Adam Michael Johansen. Momentum particle max- imum likelihood. In Proceedings of the 41st International Conference Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 29816-29877 PMLR, 2024.
Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose Bayesian inference algorithm. Advances Neural Information Processing Systems, 29, 2016.
William March, Bo Xiao, and George Biros. ASKIT: Approximate skeletonization kemel- independent treecode in high dimensions. SIAM Journal on Scientific Computing, 37(2):A1089 A1110, 2015a.
William B March, Bo Xiao, D Yu Chenhan, and George Biros. An algebraic parallel treecode in ar- bitrary dimensions. In 2015 IEEE International Parallel and Distributed Processing Symposium, pp. 571-580. IEEE, 2015b.
Victor Minden, Anil Damle, Kenneth L Ho, and Lexing Ying. Fast spatial Gaussian process maxi- mum likelihood estimation via skeletonization factorizations. Multiscale Modeling & Simulation, 15(4):1584-161 16184-6611, 2017.
Marina Munkhoeva, Yermek Kapushev, Evgeny Burnaev, and Ivan Oseledets. Quadrature-bait features for kernel approximation. Advances in Neural Information Processing Systems, 31, 2018.
Khai Nguyen, Nicola Bariletto, and Nhat Ho. Quasi-Monte Carlo for 3d sliced Wasserstein. In International Conference on Learning Representations, 2024.
NIST. NIST Digital Library of Mathematical Functions. 1.2.1 of 2024-06-15. 2024.
httpp:://dlm .ist. gov/ Release
13

---

Published as a conference paper at ICLR 2025
Emanuel Parzen. On estimation of a probability density function and mode. The Annals of Mathe- matical Statistics, 33()):065-1076, 1962.
Gerlind Plonka, Daniel Potts, Gabriele Steidl, and Manfred Tasche. Numerical Fourier Analysis. Springer, 2 edition, 2023.
Daniel Potts and Gabriele Steidl. Fast summation at nonequispaced knots by NFFT. SIAM Journal on Scientific Computing, 24(6):2033-2037, 2003.
Daniel Potts, Gabriele Steidl, and Arthur Nieslony. Fast convolution with radial kernels at nonequi- spaced knots. Numerische Mathematik, 98:329-351, 2004.
Michael Quellmalz. The Funk-Radon transform for hyperplane sections through common point.
Analysis and Mathematical Physics, 10(38), 2020. doi: 10.1007/s13324- 020-0038--..
Julien Rabin, Gabriel Peyré, Julie Delon, and Mare Bemot. Wasserstein barycenter and its applica-
tion to texture mixing. In International Conference Scale Space and Variational Methods in Computer Vision, pp. 435-446. Springer, 2012.
David L. Ragozin. Rotation invariant measure algebras on Euclidean space. Indiana University Mathematics Journal, 23 (3(11):1139-1)) 1974.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in Neural Information Processing Systems, 20, 2007.
Murray Rosenblatt. Remarks on some nonparametric estimates of density function. The Annals of Mathematical Statistics, pp. 832-837, 1956.
Boris Rubin. Notes Radon transforms in integral geometry. Fractional Calculus and Applied Analysis, 6(1):25-72, 2003.
Boris Rubin. The X-cosine transforms with odd kernel and the hemispherical transform. Fractional Calculus and Applied Analysis, 17(3):765-80) 2014. doi: 10.2477/613540-014--0333
Nicolaj Rux, Michael Quellmalz, and Gabriele Steidl. Slicing radial functions: a dimension walk in the Fourier space. Sampling Theory, Signal Processing, and Data Analysis, 2025. doi:
10.100776443770-025-0000000
John P Ryan, Sebastian E Ament, Carla P Gomes, and Anil Damle. The fast kernel transform. In International Conference on Artificial Intelligence and Statistics, pp. 11669-11690. PMLR, 2022.
Bernhard Schölkopf and Alexander J Smola. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press, 2002.
John Shawe- Shawe-Taylor and Nello Cristianini. Kernel Methods for Pattern Analysis. Cambridge Univer- sity Press, 2004.
Yu-hsuan Shih, Garrett Wright, Joakim Andén, Johannes Blaschke, and Alex H Barnett. cuFIN- UFFT: a load-balanced GPU library for general-purpon nonuniform FFTs. In IEEE Inter- national Parallel and Distributed Processing Symposium Workshops (IPDPSW), 2021. doi:
10. 10.110//PPPP55591.2222..000
David Slate. Letter Recognition. http:://dooorr/110244422/////
UCI Machine Learning Repository, 1991. DOI:
yl Meerovich Sobol'. On the distribution of points in a cube and the approximate evaluation of integrals. Zhurnal Vychislitel' noiiteellnoniiiii Matematiki Matematicheskoi Fiziki, 7(4):784-802, 1967.
Ingo Steinwart and Andreas Christmann. Support Vector Machines. Springer Science & Business Media, 2008
Danica J Sutherland and Jeff Schneider. On the error of random Fourier features. In Proceedings of the Thirty-Fira Conference on Uncertainry in Artificial Intelligence, UAI' pp. 862-871. AUAI Press, 2015.
14

---

Published as a conference paper at ICLR 2025
Gabor Székely. E-statistics: The energy of statistical samples. University, 2002.
Techical Report, Bowling Green
Eugene Tyrtyshnikov. Mosaic-skeleton approximations. Calcolo, 33:47-57, 1996.
Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Couraaeaa,, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al. Scipy 1.0: funda- mental algorithms for scientific computing in Python. Nature methods, 17(3):261 2722, 2020.
Holger Wendland. Seattered Data Approximation. Cambridge University Press, 2004. doi: 10.1017/ (B0978051161753)
Christopher KI Williams and Carl Edward Rasmussen. Gaussian Processes for Machine Learning, volume 2. MIT Press, 2006.
Robert S Womersley. Efficient spherical designs with good geometric properties. Jn Dick, F. Kuo, and H. Wozniakowski (eds.), Contemporary Computational Mathematics A Celebration of the 80th Birthday of lan Sloan, pp. 1243-1285. Springer, 2018. doi: 10.10077978-333119-446--- 07.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-MNIST: a novel image dataset for bench- marking machine learning algorithms. arXiv preprint 1 2017.
Changjiang Yang, Ramani Duraiswami, Nail A Gumerov, and Larry S Davis. Improved fast Gauss transform and efficient kernel density estimation. In IEEE International Conference on Computer Vision, pp. 664-671 vol.1, 2003.
Changjiang Yang, Ramani Duraiswami, and Larry S Davis. Efficient kernel machines using the improved fast Gauss transform. Advances in Neural Information Processing Systems, 17, 2004.
Norman Yarvin and Vladimir Rokhlin. An improved fast multipole algorithm for potential fields on the line. SIAM Journal on Numerical Analysis, 36(2):629-66 -66626966666666 1999.
Felix Xinnan X Yu, Ananda Theertha Suresh, Krzysztof M Choromanski, Daniel N Holtmann-Rice, and Sanjiv Kumar. Orthogonal random features. Advances in Neural Information Processing Systems, 29, 2016.
15

---

Published as a conference paper at ICLR 2025
A BASIS FUNCTION PAIRS (F,
We include a list of pairs of basis functions (F, f) that fulfill the slicing formula and the correspond- ing Fourier transforms in Table 2 The pairs are taken from (Hertrich) 2024 Table 1). We use the convention that the Fourier transform of function E L, (IR)) is defined by
Falg](() = / g(x)e 2ni(w,,)
((x))==2m(u.m)))) da.
(11)
The basis functions from the table involve couple of special functions, which are defined fol- lows:
Gamma function: [(2) = So - dt, Re(z) 0,
Modified Bessel function of first kind: (x) -m=o mll (m+
2m
Modified Bessel function of second kind: (x( [an(s) (a)
sin(ow)
Modified Struve function: (x) = Em= Fim+ (m+++++++)))
2m+0+1
The formula for the Fourier transform of the Matérn kernel (and thus for the Laplace kernel with = 1/2) can be found in (Williams & Rasmussen 2006, (4.15)). Consequently, we recover the well-known results that the Fourier transtorms of the Gauss, Laplace and Matérn kernels are the Fourier transforms of the Gauss, Cauchy and Student-t (with 2v degrees of freedom) distribution. To compute the Fourier transforms FT- L [)], we apply (Rux et al. 2025 Prop 1), see also (Hertrich] 2024 Lem 6) for the Gauss kernel. For the Riesz and thin plate spline kernel, the Fourier transform does not exist in a classical sense, but only in a distributional one, see Rux et al. (2025) and (Wendland 2004 Sect 8.3) for details. For the sliced Laplace f(x) exp(-aaz -ax), we have by Herrtrich] 2024 (3)) and (Gradshteyn & Ryzhik 2007 3.387.5)
F(t) [ exp(-ast)( - 82 -) ds
T
(I (1aي2(-at) Lage (-ae(-at)) (12)
(at)
Table 2: Basis functions F for different kernels K K v) = F(I vll) and corresponding basis functions I from k(æ, v) [(x vl). We added the inverse Fourier transforms Fa [F II)] and F D] to the table.
Kemel
Flat
Mil
Gauss
expl- #
Laplace
exp(-out
Σ mop(-
Slicel Laplase
T7
Muzém
aire (at.))xKKK (Sgra))
Hertrich Appendix
Eigino
Riese for
mat"
nOt Function
not function
loafa) C with Ce (M 2 log(4)
Thin Plate Spline
log(s)
not & Function
not function
B SPHERICAL SOBOLEV SPACES
We briefly introduce spherical harmonics and Sobolev spaces following Atkinson & Han 2012) 2/2
We denote by | = T(d/2) the volume of the unit sphere gd-1, Let Y k = I,... Nn,a be
an Lotorrhonormal basis of spherical harmonics of degree E Np. i.e., harmonic polynomials in d variables that are homogeneous of degree n. Here, the dimension of the space of spherical harmonics of degree n is
(2n + d - 2) d 3)!
2
Nn
for n + oc.
(13)
n!(d - 2)
(d - 2)!
The spherical harmonics V No. h 1,..., Nn form orthonormal basis of (Sd
16

---

Published as a conference paper at ICLR 2025
The spherical Sobolev space HS(S=1) for E R can be defined as the completion of C (P(()))) with respect to the norm
Naad
(n+ 2
Kg M L2(84 1112.
(14)
==0
where (g. (2) L2184 1y = fga-1 (E)) Note that the factor (n + dz2)2d can be replaced by another with the same asymptotic behavior with respect to n, yielding equivalent norm. For s = 0. we can identify (Sd-11) with L2($4-1, The Sobolev spaces are nested in the sense that H=(S--1) C (S whenever $ t. If > d-1 each function in the Sobolev space 1*($d-1) is continuous (more specifically, it has a continuous representative). If s is an integer, then (I consists of all functions whose (distributional) derivatives up order are square integrable, cf. (Dai & Xu, 2013] Sect 4.5 and 4.8). An alternative characterization of the Sobolev norm uses the Laplace Beltrami operator A., which consists of the spherical part of the Laplace,
2
and is given by
( = (- + (
If is even integer, the operator applied to g is usual differentiable operator, otherwise it is pseudodifferential operator.
C PROOF OF THEOREM
i): Since F is continuous and positive definite, so is I' by (Rux et al. 2025 Corollary 4.11). Further, because Ff (E, 2)] F[[ and due to the fact that positive definite functions are maximal in the origin, we deduce
Valfl = Bg- [F((()) 2) F(I)=1)) z -
ii): We write & (Ei, ..., Ea)² where & E [-1, 1] denotes the d-th component of &. We assume w.l.o.g. that x xea with x + 0 and recall that we can decompose the unnormalized surface measure on gd-1 as
ds (E) d (1) dt,
where & /I tea and n gd-2 x {0], see (Atkinson & Han 2012 (I.16)). Consequently,
| 2r
Eg-ula-11 (1(5
Ed 2r asdd-1(())
3d-11 gd
| 2r
11/21(1
1) d
2 (())
|
2r
Sd Sd-21 [
(1) 2 dl.
3d-1 gd
Further, with B(.,)) denoting the Beta function and the substitution u = t2, have
[
2r Pr - - dt
//' -r-b (1- - u) 4 -u
du = r 2' 2
d
r (r +
r ( +
(15)
and
gd-2
2x -
(
P (
8d-1
I T(d-1) 2T
VRT (d-1)
Combining both expressions yields
T(r + 1)) T()))
E Ee
[[(5, 2) = ||x||
(16)
[(r + ) /T²
and hence
dtr r
|> / (r }
r d2rr
2 Be-ilde [|(2)))121]
(d)
2
(
r+1) 2
r( )2
[r++
17

---

Published as a conference paper at ICLR 2025
The asymptotic expansion + a))[(()) 2n for 2 oo, see NIST 5.11.12). implies that
r( (r))((++
lim d
d
r( (())
On the other hand, yields that
is increasing with respect to d because the identity I( 2I)2)
I( r( (( 4
I d )
(dtr 2
(d r)
2dr + d2 1.
T [r + )
[r + ()
Hence, we see that
r( =(
11
(r) 2
and thus finally
mI
(dtr) [
Va[f] & Eg
[f([(E) 2)1
Eg
Eg--laa- 2))
T (4)² 2
(++1 2
ViI(r + 2
||]||.r
)(
which proves 77.
iil): We move on to the thin plate spline kernel with f(1) = d/e/2 [og((()))) + Calds As above, w.l.o.g. we assume that x sea with 2 O. Then
E [((
(ds log(sEa) C2s (1 <2) 2 dEa
2
gd-
d²
( 8 sA5 (2log (Ea) + +-2 log(s) + Hg - -2 + log(4)))
2 (1 (1- d-3 dcd
V d-1
2
dr
/
4
4
E (E.) -3) d-3 dEd
(17)
2/
+ +4(2 (s) + H - 2 ++)g((4)) [' log(Ea)) -22) dea
dea
+ (2log(()) + H4 -2 llgg())) 2 / c3(1 -3-1-66 (2) dEa)
We analyze the terms separately. Denote by (b(()) the n-th derivative of the digamma function and by y 0.57 the Maleer-Mascheronn constant. In the following, we use the asymptotics
H (x) + 7 0(x ))),
(
log(æ) 00= (x) - O (z
By (Grddshteyn & Ryzhik
2007 4.261.21), we have
Ed (lgg)((())))
d-3 <2)
T
dea
3VFT
( - ( + (1)10) - (0 ( 2
32r d+4
3/
(log(d)² - (1,(0) (3) + +lg((2 + 0(H)) log(d)
32r (1+1)
+10(11 (8)
(0)))) (5) log(2)) *0 -0(1)))
18

---

Published as a conference paper at ICLR 2025
and L." (og((E)) (6-)) -23) -
VTl
d=1
ded
(-8+ 1H
1+4 + log(64))
161
(++)
2
VFl
($=1)
(3log(d) + - 8 log(8) O(d
16T (d+4)
Recall from 15 that
S. <1 22 A-3 dea
T 1
-
88
($111
Furthermore, using that log(4) () = 2 (2g(2), we obtain
2 2log(s) + Hg a - 2 + 1o) = 2 log(s) ()) + log( d) + + log(2) - 2 Old
(2log(s) - Ha - 2 + log(4)) = (2) + log(d) + 7 + l2) - 2 Odd
logg(s)² + log (d)² + 4 log(d)
+ 2(y + log(2) - 2) log(d) + 4 log(s) + log(2) - 2)
+ (7 + log(2) 2)² + O(dd)).
Plugging this into 17 yields
d+4
1=1
1/d E [[(({)=>)))]]
E [[((KS==)1))]
V π
d2s
s+
3
log(d)) - (//0) + log(2) + +(1)) log(d) + (5) + (3)+)) + +0())
4
1 (2log(s) + log(d) + 7 + log (2) - 2 0(5)) (3log(d) d) + 3 - 8 + log(8) + 0(z))
3
+
(4log(s)² + log(d) 2 + 4log(s log(d) - 2(n + (2g 2) log(d)
4
+4log( () + log(2) 2) + (og+ 2) 2 0
3 3 3
log(d)²
[ 4
+ log(d))
2 (1))(()) log(2) (
3
1
(3y - -8 + +log(8))
2
3
3
(()og(()) + + log(2) - 2) (4log(s) + 2(7 + +og(() - -2))
+ 3 log( s)² log(s) (37 8 + log(8))
3
1 (() + (2)0)))))) (3) + log(2)) ))
+
100())))))
4
log(d)
(-3
(B) 2010
3
3 2 -2llog(2) + +0(2)
+ +3 (og(s))))
+ log(8) -8)
3
(3)
(2)1)( + 1og(2)) 2(y + log(2) - 2)² + Odd
+
4
With the identity ()) -2 log(2) + $ and
= 37 - log(8)
+ log(2)) 69 := 2 (1) + (E)
3
+ log(2) - -2)2 2.895. + +
(18)
we obtain
+ 2/d
s5 E [f(Ks x) |> = {s) + c1 log(s) + + 0 log d)
19

---

Published as a conference paper at ICLR 2025
iv): By (Hertrich 2024
Thm 3), the transformed function has the form
f(s) = Σ
00
bn
with
bn
n
r
an.
7
n=
Moreover,
f(s)² =
sas
by bebn-
gr
n=
n=
h==0
and, by [6).
(1))
E -g [f(I(S x>1]]
Cn 1=|
n=0 VTl' (
Using that d is odd and applying the identities ['(z + 1) = [z))z and I( VT, we have
r ( r (1t)
k+d
Ckin.d-
VTI (ntd)
r
("=
I [ (+t)
I ((t++d-1)
|T
n=kt]+d-1
2
(
r
k+l 2
r (n-k+1)
(k+ 2j - ()(n))n(-- k -))
(n+ + 2j - 12)
j=1
Hence, we obtain
8
Valf(()) = Σ
Σ
[kknn.d - 1)a ||=||"-
n=0
\km0
We apply the identity 4zy (x + y)² (x - -y)² to the numerator (with k: 23 1, y
n2 k + 2j- - 1) and denominator (with x = n + 2j - 1, = 2j - 1), and obtain
(ks + 2 - 1) - k + 2j - 1)
(n + 2222 11)- (n (n 2k)²
ddd
1n+ 1)
2 - 1))² n²²
j=1
j=1
1=1
n² -² 22x
k(n - -k)
1+
(n + 2(2) - 1))2 n²
(2j - 1)(n + 2j 1)
==1
j=1
For F the Laplace kernel, assuming w.l.o.g. that a = 1, we have an := (-1)a/ /n!. Consequently,
n k(n - k) 1
x
| 1)
n
Va[s((>)
kl(n - k))! (-1|=|1)" n n+1
Σ
(km - k.2)
(n + 1)!
n=0
(-1==)>>
(1²22-11 nnnn + 1)2" 2)
1
(- 1)
n(n 1)
(n + 1)!
(n + 1)!
(-1||=])))
1
(- "
(-2
(n + 1 - 2)m
(19)
(n 1)!
7 n=2
1)!
(n + 1)!
n=
For the first term, we obtain
[-2
(n 1}! = -1|/=-/| Σ
211
-2 |||1||| (0 (0-211=11 -1) 1).
(20)
n!
2
and for the second
20
2 | 10
-||x=|>0*
x
(-211==))*
(n + 1 - 1) =
Σ
(n + 1}!
n!!
(n + 1)!
2=2
n=2
n=2
1
2||=||
1 2
||c|||| +
|||=||
1 + |||x=||
(21)
2||z|||
20

---

Published as a conference paper at ICLR 2025
Plugging 20 and
21
into
19 finally yields
|2|||
1
e |||||||
Vz(f((z)
(e
2||=||
2||=||| -1+ 2||2=||
+
+ 1 - ||=|1]
2
2|10|
2|/21
2|1
(22)
(2||=1||
2 + 2 2||2||| + 1))
4
|0:||
|||
(2 |0|
+ 2||3||| + +)) F F
4||ar||
For the Gauss kernel, we can follow exactly the same lines, after considering w.l.o.g. o² = 1/2 and replacing || with ||==||?
D DISCUSSION OF THEOREM
D.1 COMPARISON TO THE ERROR BOUNDS FROM HERTRICH (024
In the following, we give a short summary how Theorem improves the error bounds from Hertrich 2024
For positive definite kernels and Riesz kernels, the bound from Theorem | is dimension- independent. In contrast, Hertrich (2024) only proves dimension -indeeidenndeeeeeeent absolute error bound for the Gauss kernel (and conjectures that this also true for the Laplace and Matérn kernel). For the Riesz kernel, the error bound in Hertrich (2024) depends on the dimension by (Vvd).
Theorem bounds the error of thin plate spline kernel. For this kernel Hertrich (2024) does not provide a bound.
For kernels with analytic basis functions and for the Riesz kernel, we exactly compute the variance Wa[f]. In particular, Theorem provides an exact calculation for the mean square error, i.e., no tighter estimation possible. E.g., for the Gauss and Laplace kernel, this yields an improvement of the absolute error bounds given in Hertrich (2024). observe that this bounds actually decay to zero for || - DC.
D.2 THE VARIANCE OF GAUSS AND LAPLACE KERNEL FOR d >
For d > 3, d odd, evaluated Wa(f) for d 5, ..., 15, and make the following conjecture. Let TT(J) = Σn k= ((6)(() * denote the Taylor expansion of f in zero up order We conjecture it holds that
that, for d 2 3, odd, and the Laplace kernel F(x)
2a|=
(a])(a) ca(a{|[=]))03
Σ (-1)
Ta- (22
b b (allx)))))
where
Cd-2 Co.d (d 2) 4)co0 for 2 with
Ca (12+ + (12+4(d- 4(d-)))) cd-2 for d 5 with
C3
Co.3
1,
7 with b0.5 = 4. ba-6 2ba_b 4-2 ,d-2 for d 2
For the remaining coefficients, we have not found general simple rule. Our numerical computations indicate that they are positive and that [z)f(()) [x]]] | is monotonically increasing in ||]||, with upper for 0. The upper bound sa increases slowly in d,
bound sa/a and quadratic convergence and our simulations hint that it might be bounded by constant independently of d, see Figure 3 For the Gauss kernel, the variance has the same form, except that a -1 replaced with a o/v2 and ||>||| with ||= all occurrecces.
21

---

Published as a conference paper at ICLR 2025
0.35
0.25
(XI))
IIII"Vs(foXxx)
0.15
IIxIFVs(faXxx)
iI(I"\
0.05
y=1/3
2
4
The scaled variance of the Laplace kernel fz with 1 (in red) and the Gauss kernel fo with = 1/v2 (in blue) for dimension d = 3 (solid lines), d = 9 (dashed lines) and d = 15 (dotted lines). The variance is multiplied times ||=||| (Laplace kernel) and |>=|| (Gauss kernel). We observe that the scaled variance increases monotonically in all cases, scemingly bounded from above by constant (black solid line).
Figure 3:
E RANDOMIZATION OF A QMC SEQUENCE
Having QMC sequence (5") see Section 3.1 we can easily obtain an unbiased estimator in 44 by considering the QMC sequence (AEP), P with a uniformly chosen orthogonal matrix A (o(() because, according to (Ragozin 1974 (2.3)), we have
BA~Won (A55)) = [ Joran (dEE))dLoog JdLooa - ((n)dM4s-s1)))
=
where Mo(d) is the uniform distribution on the set ((d) of orthogonal x d matrices. Furthermore, note that $ f(EP) converges Jgs f(E)dl4g-+1() for P oo and all continuous func- tions f because 1Sdd dense in C(Sd-1) and the weights are all one, ef. (Atkinson 1991 Sect 5.2).
F PROOF OF THEOREM 3
Gauss kernel: A sufficient criterion for the function being in the Sobolev space (gd-))) is that S(I- is times continuously differentiable. For F(() exp( 12 ) the one-dimensional basis
function is given by the confluent hypergeometric function ftt) = Fi
_oz) We know by Herrtiic 2024 Lem 6) that f FF(g) with g(w) doo exp( - 2 R 22 W 13/2 where Fi
-t²
3:
/2 212)
denotes the one-dimensional Fourier transform. Then, (w) = * ((w) is absolutely integrable for any l E N such that the pnterrentiaiiommmutippiiaaa formula for the Fourier transform yields that the k-thh derivative of ff exists and is given by the continuous function ((k) (2mi)* Fi()).
Riesz kernel: With fixed x E Rdd {0}. we have (E)
<),=))". We set ==
|>ा E
gd-1 then || and accordingly (E) || (E). The Sobolev norm 14 reads
80 N
(n + 1=2)
20
||=|||
/
(E)))
The last integral can be computed with the help of the so-called A-cosine transform. By (Rubin 2014 (3.5) and (3.8),, the X-cosine transform of a function h E (Sd-11 for À -1 is defined by
(
/
c^(N](w)
(w,)
Ad,
T A+1
22

---

Published as a conference paper at ICLR 2025
for w E Sd-1 and satisfies
(-1) r(
c* (w)
even, odd.
0.
Hence, we obtain that
| (E))
| |
d+r T(
S
Y
(1) = (-1)
(,((n)
r(
((=td ((=td-r) (--3)
n is even and O n is odd. The addition formula for spherical harmonics,
No
Nn. 84-11
Jmm
see Atkinson & Han 2012), and the substitution
2m yield
d-2
2
2s
I(m -221((
2
l 1 (844 (84 1)
2m +
||a|||
2r
| -
2
T(m + 2 r
The asymptotic relations Nn.d d- from (3) and [(2 + a)/((=) 2 for - oc, see NIST 5.11.12), yield that the m-th summand behaves asymptotically like multiple of (2m) 220-2r22 27-2 therefore the series converges and only if < <
Matérn kernel: C.I)
The one-dimensional basis function of the Matérn kernel is (Hertrich 2024
Appx
T(1 () (2
d. ((t) = 1
2v
F2 v S. v
vt 1 2
2:2
2
[(2) 1)32
2
(2())
=ii(()
fs(()
The above hypergeometric functions and fa, whose parameters do contain any nonpositive integers, are entire functions in R. Hence, the corresponding spherical functions (91). and (ga)z from (1O) are in C (Sr-1) because they can be extended smooth functions defined on neigh- borhood of 1d-1 Furthermore, since fz is multiple of the one mie-iinesio basis function of the Riesz kernel, we know from above that (92)z E (Sd--11) for
+ Hence, by (Quellmalzy 2020 Thm 5.2), we see that the product (gs) is also in H*(S*-1), )). Note that the referenced theorem is only formulated for in integer but this be easily extended using interpolation argument (Quellmalz) 2020] Sect 5.3), because the multiplication with smooth function constitutes
a continuous operator both H|s] (Sd-1) + H|* (Sd-1) and H [x]+ (gd--11) H|*) (gd-1)).
G BACKGROUND ON RFF AND RELATION WITH SLICING
We denote by (Rd) the space of finite positive Borel measures on R.d. Such measures can be identified with linear functional on the space Co (R)) of continuous functions that vanish at infinity. The Fourier transform of measures is linear operator defined by
Fa: M+ (R + (Rd)), Falu((x)
e 2 dp(()))))dpp)))
cf. (Plonka et al. 2023 Sect 4.4). By Bochner's theorem, the Fourier transform is bijective from (RR)) to the set of positive definite functions on IRd, see ii) Section for the definition. If y probability measure, i.e., (( 1, we have (a) = 1.
the following, let F be a positive definite function on Rd with F(0) 1.
RFF: Random Fourier features (RFF), see Rahimi & Recht (200)), use that, by Bochner's theorem, F(II []) is the Fourier transform probability measure μ E M+ (Rd), i.e.
F([=[]) Fz)) Evnn (exp(2 v))
D Σ exp exp((2ti(z, up)),
(22)
23

---

Published as a conference paper at ICLR 2025
where sample D E N instances Up iid from μ.. Using that F is real-valued, we replace the exponential by the cosine in (22). By (Altekrüger et al. 2023| Prop C.2), the radial measure μ can
be decomposed as follows. We define RTT + [0, 00), (x) ||||| and its pushforward measure
ji = baf4 = po port E EM+ ((0, o0)), then we have
T(p @ Ugr-1), where [(r, = rG.
Hence, sampling from j be realized by Up rppp- where &p Uga-s and rp p.. Then 22 becomes
D Σ
F([|=]) ~
(2T Ep)).
(23)
D
p=1
RFF Summation: For approximating the kernel sum we insert the RFF (2 and obtain
m
Wn
2m -.n
DDD
W 2mi(kn.pp)
(24)
where U1. UD are iid samples of Since the inner sum over is independent of m and therefore has to be evaluated only D) times, the total computational complexity of computing (24) for all m 1, ,..., M is O(D(N + M)).
Slicing: The slicing approach uses the approximation i.e.,
F([ ~
Σ ((Kx,))))) P
1
where Sp ~ Lgga-1. By (Rux et al. 2025 Cor (.14).1) the function f(l-l) is positive definite and hence possesses Fourier transform, which probability measure on R because f(0) F(0) 1. Applying RFF with Q points to the one-dimensional function f(1 (), we obtain
F([
cos 2TT (x, Eph)s
(25)
PQ
where Tna FT [f( D).. According RRux al. 2025
Cor 4.11), have
T Aic. Ajū,
where Ai: ( 20)) (R) is the symmetrization operator that extends measure on [0, oo) to an even measure on the real line and is defined for any # ([0, oo)) and E (R)) by (Ajv, g) (v, + g)- Since the right-hand side of (25) independent of the sign of Tp. it stays the same when we sample pp.g from jú instead of A p. Therefore, if Q = 1, we see that the
right-hand side of (25) the same as the right-hand hand side of (23 In particular, RFF can be viewed
as a special case of slicing. The respective slicing summation given by (5) combined with the one-dimensional summation in Appendix below.
H RELATION BETWEEN DISTANCE QMC DESIGNS AND ORTHOGONAL POINTS
Consider the QMC design _ that is a minimizer p for ersley (0188 suggested symmetric QMC designs, meaning that for every point contains also the antipodal point -E. However, since the function (TO) we want to integrate is symmetric, i.e., (E) Iz (- we can discard one of the antipodal points and - and get the same result. Therefore, we minimize the functional
To improve the error, Wom-
(sym(E") =E()= -(()""-E')))
Σ (IE p-9=1
=
11 + +1(1))
24

---

Published as a conference paper at ICLR 2025
Our numerical trials indicate that indeed the minimizers of Esym yield smaller integration error than the minimizers of E. Using l1Epll = 1, we see that
-I1Ep - E + 1 <(11)² (V2 2- (g) + 2+ 2 2(55, -
= -44 - 2/ 4 55 E
attains its minimum if and only if <<p sp) ) = 0. Hence, P S d and & is an orthonormal system in Rd, then e minimizer of Eyym- However, this argumentation does not work if P > d, as we can only choose d orthogonal vectors.
I BACKGROUND ON ONE- DIMENSIONAL FAST FOURIER SUMMATION
In this section, review literature about one-dimensional fast Fourier summation used in Sec- tion 4.3 and specify the parameters used in our numerical examples. Fast Fourier summations were proposed in Kunis al. (2006); Potts al. (2004) based the non-equispaced fast Fourier trans- form (Beylkin 995 Dutt & Rokhlin 1993). which is implemented in several libraries (Knopp et al., 2009 Shih al. 2021). In our numerical examples, we use the Julia library Here, we follow similar workflow as in Hertrich (2024).
2023 Keiner al.
Knopp et al. (2023).
Let 21...., TN E Rdd yi.... VM Rd W1, Wi,..., uN and k(r,,) y) f([== - vl) = g(x y). We want to compute for & E ga-1, Zn.g == (En, E) and Fm. = (ym. {) the one-dimensiona kernel sums
N
N
tm
Σ w wonk(angg Vmm.)) n=1
Σ Wng(Fn. - ym.g)-
Step 1: Rescaling For Step and 3, we will need two properties. First, since we will use discrete (fast) Fourier transforms, we require that Zn,g: Mm..; Tn. Dmm,g E [- )). For the important example of positive definite kernels, which decay for the Fourier transform of g via Bochner's 's integral and Rux et al. (2025), see Table for some examples. In order to use this explicit formula for the fast Fourier summation, we will additionally require that g(x: Rs 0 for [x| > 2'
zero, we often can derive explicit formulas
Both properties can be achieved by rescaling the problem. More precisely, let T < 0.5 and c max.=I....N || maxm=1,...MM Il |l. For the case of decaying positive definite kernels, assume that g(x) Rs 0 for [ > Imax. Then, holds that
N
Im Σ
Wng(Ins ymm.s) = Σ ung(@nn "m.5),
with
② Fngg TInA E) Dn.5 == Tyn.g = (TDn,
((x) = g
where the constant properties from above are fulfilled. For the rest of the section, we will denote the rescaled points En,g Um.g and the rescaled kernel function ğ again by Zn,gg Vmn.g and g.
min T ) does depend on E. Then, by definition, the two
In the case of the Gauss, Laplace or Matérn kernels, the rescaled kernel k(x,, y) = G(z -y) again a Gauss, Laplace or Matérn kernel with the altered parameter = TO, â = a/r or = -3. In our numerics, we set gmax 5m with a Moreover, we set the threshold T to 0.3 for the Gauss kernel, to 0.2 for the Matérn kernel and to 0.1 for the Laplace kernel.
for the Gauss, Laplace and Matérn kernel.
Step 2:
Computation of the Fourier Coefficients of the Kernel In the next step, we expand into its Fourier series on [- 1. 1) and truncate it by
Nn./2-1
((()) /
( Σ Cs ( 2nike KEZ
Σ
ca 2aik
2nike
g(a)e -xikke dx
k= N
25

---

Published as a conference paper at ICLR 2025
with some even Na N. To compute the Fourier coefficients cz(().. employ that g(z) =s 0 for |x|| > 2/2 such that Poisson's 's summation formula (see, e.g., Plonka et al. 2023 Thm 2.28) implies
(g)
= Ck
Σ E(( (k), Jez
where [W] is the Fourier transform III.
In our experiments, we choose Nn. 128 for the Gauss, N 512 for the Matérn and Nn. = 1024 for the Laplace kernel. Note that the coefficients ei (g(g) do not depend on the input points and need to be computed only once for different choices of &. The function Fi [g] is analytically given for the Gauss, Laplace and Matérn kernel in Table 2
Step 3: Fast Fourier Summation Finally, we use this expansion to compute the kernel sums
N
Nn/2-1
Na/2-1
Im
Σ Nn./2
(n 2 Fnnee)
Σ
cx(9) (
WyC 2rikzn. (26)
1=
=10
The computation of the second sum dy is the adjoint discrete Fourier transform of the vector W = (Wi,..., WN) at the non-equispaced knots (-I1.g,..., EN.g). Afterward, the computation Na/2-1 of the vector (t,, ((,..),,)))) is the Fourier transform of the vector (cu(g)(n )== -n-N the non- equispaced knots (-y1,g,... yMgg). These Fourier transforms can computed by the NFFT in time complexity O(N + No. log Na) and O(M Nn log Nn) leading overall complexity of O(N + M + Nn log Na..
J COMPUTATIONAL COMPLEXITY
We consider the computational complexity of the random Fourier features (RFF) and the Fourier slicing summation. We denote the RFF summation 24 with D features by
N
D
D
N
3RFF m
Wn
2milgm -pn:0p)
Σ
2milymm.pp
Σ
w pp)
(27)
D
D
pul
p=1
Computing $RFF for all m 5m
M has complexity of O(D(N M)).
We consider the Fourier slicing as described in Appendix with the slight modification that instead of rescaling the problem, we take 27 perioddi Fourier series. To this end, we choose
R > max n 1=1.
||x,|||
max m=1.....M
yl ||-
(28)
Then, in the one-dimensional sum in (5). it holds that (x,, - ym { S. for for all & Sd-1, Therefore, we can replace f(( 1) (5) by any function g: without changing the sum. In particular, we choose as sufficiently smooth 27 -prriodi function for some Rt in order to achieve convergent Fourier series of g. With this, we insert the one-dimensional summation (66 with a 2T -perioicc function into the sliced kernel sum 5 and
that satisfies g(t) = f((t)) for all |t| R
obtain
Na/2-1 Σ
Wr.Ck -
(29)
P
N2.22
with the Fourier coefficients ea (g) = (27)-1 f (() rioo-) dt. As already noted in Ap- pendix we interchange the sums in order to achieve a fast summation
P
2
N Ep)/T c
JSS m
Σ
Σ
w Tike
m I.... M.
p=1
Jim
Na/2
n=1
Utilizing the NFFT for the sums over n and k, this has a complexity of
O(P(N M + No. log N₂.)
26

---

Published as a conference paper at ICLR 2025
with P the number of slices and Na the number of Fourier coefficients.
The asymptotic complexities with respect to M, N the two methods depend on different param- eters, most notably the number D of random features for RFF and the number P of points/slices in the QMC sequence for Fourier slicing, which need to be chosen depending on the desired accuracy. The following proposition compares their complexity achieve relative emror
E
|||s
RFF /RFFII [|,]] E O(e) and
| ||s
1/ OEE))
for 4 0,
|
1|
where ||mm|| = n= wel and || = max taken with respect to the random Up ~ p.
/mm) and, in case of RFF, the expectation
Proposition 5. Let F(t) = exp(- - be the Gauss kernel with a > 0 and sliced kernel f. Assume there exists R > satisfying (28) independently of N and M.
The computation of the RFF sum gRFF = [S ) see (2 with the parameter
D loge) for + O achieves the relative error Ellls ///ol//|ll|1 O(s) and the numerical complexity
O(N + 2)
For q E N arbitrary, let g (RRR)) be 27 2p-rrioo function with T R that satisfies g(t) for [t] S R. Furthermore, assume that the slicing error has rate >> O, i.e.
F([==)) 1 P 2)11 EO(P))
Σ S(SS p=1
sup sERd
Then the computation of the Fourier slicing sum gts parameters P e-/rr and Nu. e=1/91logel1 for E 4 O achieves the relative error |s 3FS //////////// E O(e) and the numerical complexity
(5 mm= see with the
O(e-Vi((--Vo 1/9 + M)).
The assumption that the slicing error is bounded with rate r is fulfilled with r- 2[d-1 2 for the distance QMC designs, which minimize 4 with d/2, as we proved in Corollary However, our numerical results suggest that might be even larger, ef. Table Because E N can be chosen arbitrarily and r > 1/2, the asymptotic complexity the slicing summation is lower than for RFF. Furthermore, note that the second part of the proposition holds for any kernel function F for which t + f((t) is in (R).
Proof. RFF: By (Sutherland & Schneider 2015 Prop 3), have for the Gauss kernel
D
Eug..... Ppap
sup F(II II)
2mi((.op)
p=1
24 1/2 + /d /2ggg
Eo (1 log(D) D)
ovD
D
where By the ball in Rdd of radius T > Hence, the error of the RFF summation bounded by
N
D
B|]|s JRFF
[10
<
Wml
max m= B.....,,M M
F([IEn - Fm (l)
2mi(ym - -xaop))
m=
p=1
1
2 kn.
log D)
<
max
F( 3lll)
D
D
For the desired accuracy 6, choose D -logge(e)). Then Ing(D) (EE2), so that we obtain relative enror OEE) with a complexity of
82 Inn( Ing(e)) 22 Ing(s)
D
log(e)
O (D(N M)) O ((N + M) log(s)))
27

---

Published as a conference paper at ICLR 2025
Slicing: For the Gauss kernel, the sliced kernel + f((tl) is analytic function, Table Therefore, can be constructed via two-point Taylor approximation, see Potts & Steidl (2003). We estimate the error
Nn/2-1 Σ Nn.
|s 8 ||w||2
max
F( Vmlll)
Cllg) 20 In,
Since, R > ||0, Vml |l for all n and by 28 we have
8
3
1000
sup ||
F([[=|)
Σ f(1 f(((Ep=2))) x\1)
|w/11
p
Na/2-1
1
Σ (
+ sup |
((((5p:))1) -
P p=1
k=-Na/2
Since g(t)
f((t) for all [t] S R. we have
|s
FFS
100
sup zER4
F([
1
F x)1
sup
g(t)
Nn/2-1
Σ cx(g)erikt/T N4/2
1
p=1
For the first term, we use our assumpution that
sup ZER4
F(
1 |
Σ ((((sp;>)) p=1
E (O(P-)
with some > O. For the second term, Bernstein's 's theorem for the convergence of Fourier series implies that there exists C, such that
Nn/2-1
sup tEk
g(t)
Σ
I Cy llog (Na),
N2
cf. (Plonka et al.
2023 Thm 1.39). Combining these estimates, we obtain
|s FFS
E EO(P + N logg(N)))
Choosing P ~ e 11/r and Nn 8-1///////0/// hoel/oo(e=-) for 1 O, we see that
P + N log(Nn) ~ =(q-~log( log(E log(log(E e O(s). log(=-
Hence, the complexity of the Fourier slicing to achieve relative error E is O(P(N + M + Na log Na)) = (/16- - N M)).
K ADDITIONAL NUMERICAL RESULTS
K.1 ADDITIONAL PLOTS AND TABLES FOR SECTION
In the following, redo the experiments from Figure and Table 1 and vary some parameters. More precisely, we redo for the negative distance kernel, choose other length scale parameters of the kernel and perform it in higher dimensions (d -200).
Negative Distance Kernel We do the same experiment as in Figure and Table with the negative distance kernel. The results are given Figure and and Table We
that the advantage of QMC slicing is not as large as for smooth kernels, which is expected considering the theoretical results from Section 3 particular, the spherical function (10) not in HSd2((dd1)) d 2 3/2, so the assumptions of the bound (8) are not fulfilled. Nevertheless, QMC slicing still significantly more
28

---

Published as a conference paper at ICLR 2025
accurate than non-QMC slicing. Note that RFF based methods are not available for the negative distance kernel since it is not positive definite and therefore Bochner's 's theorem does not apply.
Other Length Scales of the Kernels We redo the experiment from Figure ! and Table with scale factors 8 = 2 and s = of the kernel parameter. The results are given in Figure and Table
for & = and in Figure 6 and Table ( for & = 2. We observe that the advantage of QMC is more significant of for larger scale factors, which is expected since the function 5 f((s, 2) is more regular for larger s than for smaller s.
Higher Dimensions Finally, we do the same experiment in Figure D and Table | for the higher dimension d = 200. Here, we use the negative distance kernel, the Matérn kernel with y/ = 3 and the Gauss kernel, where the parameters are chosen by the median rule with scale factor ? 1. The advantage of QMC is less pronounced in such
The results are given in Figure 7 and Table 6 high dimensions, but still visible.
K.2 ADDITIONAL RESULTS FOR SECTION
Negative Distance Kernel We redo the experiment from Section 4.3 for the negative distance kernel and the thin-plate spline kernel. The results are given in Figure We can see that QMC Fourier slicing outperforms standard slicing clearly in all cases. Note that RFF based methods not available for these kernels since they are not positive definite and Bochner's 's theorem does not apply.
Higher Dimensions We run the same experiment as in Section 4.3 on the MNIST and FashionM- NIST dataset without dimension reduction and therefore d 784. The results are given in Figure 9
We can see that the advantage of QMC slicing is smaller than for the lower-dimensional examples but still clearly visible for some kernels. In accordance with the considerations of Appendix H the advantage comes in when P > d.
GPU Comparison We want to demonstrate the advantage of our method in a GPU-comparison with large number of data points. As test dataset we concatenate the MNIST and FashionMNIST in all eight orientations arising rotating and mirroring the images and reduce the dimension to d 30. The arising dataset has N M 960000 entries. Then, we compare RFF, ORF, QMC (Sobol) RFF, Slicing and QMC Slicing, where the QMC directions for slicing are chosen by minimizing the distance functional, see Section 4.1 This experiment is implemented in Python using PyTorch and we use brute-force kernel summation by the PyKeOps library (Charlier et al.) 2021 baseline. The results given in Figure TO Even though the RFF-based methods parallelize a bit better on the GPU than the fast Fourier summations, the conclusions are mainly the same for the CPU experiments. We can clearly see the advantage of QMC slicing over the comparisons. Particularly, for non-smooth kernels, slicing-based methods work much hetter than RFF-based methods.
PCA
K.3 MMD GRADIENT FLOWS
Finally, we use our fast summation method in a specific application. Here, we consider gradient flows of the maximum mean discrepancy (MMD), which have been considered several papers
for generative modeling and other applications, see, e.g., Arbel et al. (2019); Chen et al. (002)); Galashov al. (2024): (Hertrich al. (2024): [im et al. (2024).
Background For dataset y (yi...,, War) of target particles, we consider the discrete MMD functional Fy: (Rd)) N + defined by
N Σ
N.M Σ R.
M Σ
1
Gz) MMD K (x. v)
K(Z,,
K(I,, U,) + 2M f-
K Uz).
2.
MN
We note that (under suitable assumptions on the kernel), the MMD is metric on the space of probability distributions on Rd such that (æ()) is always non-negative and zero if and only if the particles and y coincide. Here, we interpret the discrete probability measure 1 DN 5 8z.-
Now, we minimize G, by simulating the gradient flow
a - (2²) (0) x (0)
29

---

Published as a conference paper at ICLR 2025
starting with some random initial samples x (0) (Rd)) using the explicit Euler discretization
x (k+1) - T (
(30)
where the gradient VGy is computed via backpropagation. Throughout the flow, the partieles will then move towards the target distribution
Experimental Setup We choose the negative distance kernel (xx, y) - yll and consider the CIFARIO dataset (M = 50000 and d = 3072) as target points y. For æ. we choose N = M 50000 samples. Then, we run the (discretized) gradient flow from (30), with initial particles x0()) drawn iid from standard normal distribution. For speeding up the convergence, we follow Hertrich iet al. (2024): Lim et al. (2024) and add momentum parameter m = 0.9. That is, we modify the equation (30
(k+1) = (Gy(x + mx x [k+1] = (k) TD (k+1)
with initial value (0) O. We run the MMD flow for 50000 steps with step size T 1, where we compute the function Gy by (Monte Carlo) Slicing and by QMC Slicing with P = 1000 projections.
Remark 6. We would like to point out that running the gradient flow with exact computation of Gy computationally intractable. Using (QMC-)Slicin)) (Slicing, one iteration (30) takes between 0.2 and 0.3 seconds on NVIDIA RTX 4090 GPU. On the other hand, the exact gradient evaluation via PyKeOps takes about one hour. Considering that we are running the flow for 50000 steps, this underlines the need of (QMC-)Slicingg
Results We plot the objective value of (, versus the computation time in Figure |1 We observe that the smaller error in the gradient evaluation by QMC Slicing significantly improves the convergence behavior.
K.4 ON THE GAP BETWEEN THEORETICAL GUARANTEES AND NUMERICAL RESULTS
In our numerical part, we observe significantly better error rates than we can prove theoretically. One possible explanation is the following. Our theoretical guarantees for QMC Slicing are based on worst-case errors in Sobolev spaces on the sphere and consequently rely on the smoothness How-
of the functions Oz from (10), which is not satisfied for some kernels in Theorem ever, these results are only worst-case error rates that do not account for the specific properties of O2- First, the function (E) depends only (xr thus having lower effective dimen- sion. Furthermore, by construction (x(() (-(-E) such that for all E (E) Bg-Rl 1 (E)). Moreover, 9z is infinitely often differentiable for on the hemisphere {s gd-1, ((,x) > 0} for all considered kernels of Appendix A such that tighter error bounds could apply. Consequently, exploring QMC designs on the hemisphere could be an interesting direction for further improving theoretical analysis.
Rd it holds that
30

---

Published as a conference paper at ICLR 2025
pegetive deoance cemel of
tegative delance temel
negetive disance
16
Ro.ae 20-1
Send stond
Senc alchn
201
Figure 4: Loglog plot of the approximation error in versus the number P of projections for the negative distance kernel. The results are averaged over 50 realizations of & and 1000 realizations of æ. We fit regression line in the loglog plot for each method to estimate the convergence rate, see also Table 3
Table 3: Estimated convergence rates for the different methods and the negative distance kernel. We estimate the rate by fitting regression line in the loglog plot. Then, obtain the estimated convergence rate O(P-r) for some > 0. Consequently, larger values of correspond to faster convergence. The resulting values of are given in the below tables, the best values are highlighted in bolt. See Figure for a visualization.
Negative distance kernel
Slicing-basee based
Dimension
Slicing
Sobol
Orth
Distance
spherical design
d
0.49
0.94
0.55
1.27
1.29
d == 10
0.50
0.72
0.50
0.71
d = 50
0.50
0.69
0.65
0.70
31

---

Published as a conference paper at ICLR 2025
Gouss with 00.000000000000000000000
Gouds reme) with of -0.485.
12
20%
13
12-
AV
-2-
56
1828g
Sant
Scod
13-1
Sendtatiing stoing
Sanol
Obcoree 1ikh
Drn decrs
12-4
Dra tie
Sunde eicecttio
Proeetor
Nater cernel with p- 0.15, pm 35.dm3 d
komd
d'
with
A11
Dff
Sptelstonn orn sidng
diiing
diding
Sani dicing
Scnce slcing
11m
sprotce drace
Dittacce sidig
10
Dittance cing
Matern kemil with 0.15 uw
Matern kemd with 0.30, uw di
Matern karnal with 0.70 5.5.0
10-1
30-1
20.1
20-4
stong
Ofr
Eittace
Ornulees
Ore
Lasiace kerna 326 20
Laplace
30.1
20-1
12
Sobul
diana
Gdoecal coiie
Didaace cocing
1.0
Figure 5: Loglog plot of the approximation eтror in (4) versus the number P of projections for different kernels and dimensions (left d = 3, middle d = 10. right d = 50). The results are averaged over realizations of & and 1000 realizations of æ. The kernel parameters are set by the median rule with scale factor 5- We fit a regression line in the loglog plot for each method to estimate the convergence rate, see also Table
Table 4: Estimated convergence rates for the different methods. We estimate the rate by fitting regression line in the loglog plot. Then, we obtain the estimated convergence rate O(P-") for some > O. Consequently, larger values of correspond faster convergence. The resulting values of r are given in the below tables, the best values are highlighted in bolt. The kernel parameters are the same as in Figure (scale factor 7
scaling
scaling Slicing Nesed
Ri9--aasdd
Sticing hased
RE-H-bosed
Dinession
Sobol
Sliciss
Sobol
Disunss
Dinession
Slicng
Sobol
Disunee
sphesical
3
11:00
11:00
0.7
2.00
0:40
0.06
0.50
2:11
4.01
0.50
0.55
0.90
0.50
0.27
0.50
1.30
d
0.40
0.30
0.50
0.74
0,50
1.12
99
0150
0.77
0.67
0.50
072
0.70
0.77
d 50
0.57
0.57
0.30
0.08
0.63
0.70
Munirr kemmll with
komal scaling
Laphacs kernel with Kornel wcaling
BFF-Iuacd
Slicieg-tacc
REF-Isacd
Sticin--back
Dinession
Ori
doiign
desgg
2 3
1150
(1.33
0.51
0.96
0.50
2.11
2.24
3
1150
01.30
6.40
11988
11:22
1.35
3.25 28
10
(.50
0.50
0.50
0.20
0.50
0.88
10
[Ug
0.50
0.50
0.64
(U:0
0.60
0.58
0.00
0,51
0.50
32

---

Published as a conference paper at ICLR 2025
Gouss semol wwt.://wwhhoo.omm. w00.040...000..ocooom0 6000
Gouds reme) with of +0.408. 50
12
26"
13
12-
AP
20.4
06
Oaf
1858g
Sant
Scod RP?
23
Sendtitiing stoing
Sanol sloing
Obcoree
20-7
Drn deirs
15-4
Sunde Aicectiio
Natem cernel with p 0.15, pm dm
komd
with
134
111
Dff
Scbelstong
diiing
diding
orn siing
Sani dicin
Sonde slcing
110
sprotce drace
Dittacce sidig
10
Dittane digg
Matern kemal with 0.15 ow
Matern kemd with 0.31 ow
Mater karnal with 0.71 2.5.0
30-1
20-1
20-4
stong
Off
Eittace
Ornulees
Ore
Rindor
6.55
Lasiace kerne a- 321.
Laplace with gm
30.1
20-1
EXT
Sobul
Wling
diana
Gdoccaa coiie
Didaace coring
1.0
Figure 6: Loglog plot of the approximation eтror in (4) versus the number P of projections for different kernels and dimensions (left d = 3, middle d = 10. right d = 50). The results are averaged over 50 realizations of & and 1000 realizations of æ. The kernel parameters are set by the median rule with scale factor 2. We fit a regression line in the loglog plot for each method to estimate the convergence rate, see also Table
Table 5: Estimated convergence rates for the different methods. We estimate the rate by fitting regression line in the loglog plot. Then, we obtain the estimated convergence rate O(P-") for some faster convergence. The resulting values of
> O. Consequently, larger values of correspond r are given in the below tables, the best values are highlighted in bolt. The kernel parameters are the same as in Figure 6 (scale factor 7 2).
walling
wollng Sicing--asol
REF-bbaed
Sliciag-hasll
REF-baaod
Dinession
Slicng
Sobol
Disunce
sphesical deign
8
0.50
1,00
0.61
0.61
0.07
0.90
2.08
0.50
0.87
0.54
2:11
4.02
1150
3.85
0:30
0.50
0.7
0.50
1:37
d
0U50
0.50
0.50
0.78
0,50
1.12
0.50
0.70
0.66
0,50
0.72
0.70
0.77
d 50
0.56
0.57
0.30
0.67
0.03
0.71
Matém kenael with
ana kennel aciling
Laplace kemnel with kemel scaling =
bused
Slicing Resed
Slecnp bosed
Dimension
RBY
Schol
Clrb
Tiitanco
design
Dimarziion
RFP
Slicing
Sobel
Clrbh
Piieanoo
dedgg
2 3
(U50
0.49
0.50
0.08)
0.54
2.100
2.24
S
0.50
0.51
0.00
(1.5%
1.201
1.28
to
1151
01:31
0.30
0.02
0.50
1.854
10
1151
01:00
0.50
0.02
(1.50
1..09
0.50
0.54
0.50
0.00
33

---

Published as a conference paper at ICLR 2025
distance kemel
Gauss keme with 1,995 d00
20
De4 siliin
wlogg
Diyy diiigg
D11h dliing
alling
Number Proection
Wirder refeeeeeee
Number witeeectees
Figure 7: Loglog plot of the approximation eπrоr in 4 versus the number P of projections for = 200. The results are averaged over 50 realizations & and 1000 realizations of x. The kernel parameters are chosen by the median rule with scale factor = 1. We fit regression line in the loglog plot for each method to estimate the convergence rate, see also Table
Table 6: Estimated convergence rates for the different methods and kernels for d 200. We estimate the rate by fitting a regression line in the loglog plot. Then, we obtain the estimated convergence O. Consequently, larger values of correspond to faster convergence.
rate O ((P")) for some The resulting values of r are given in the below tables, the best values are highlighted in bolt. The kernel parameters are the same in Figure 7 (scale factor 7 1).
Dimension d 200
RFs-b
baaee
Kernel
RFF
Sobol
ORF
Slicing
Sobol
Orth
Distance
Negative Distance
0.50
0.68
0.69
0.70
Matérn, v 3
0.60
0.59
0.50
0.67
0.67
0.68
Gauss
0.50
0.72
0.72
0.71
0.73
0.76
negative distance on Letters (d 36)
10
Scrtie Sicrn
Sarding Sliing
Serting Skiiig
in dngle th cru
Time 30 digieetiodee sivededddddd dku
Time 30 on thigg eiwddddd CPu
thin plate spline on Letters 16)
thin plate spiin MMST (d 20)
thin plate spline on MashionMNisT id
109
100
Pecatver Eran
20
OMC Fourer Sioing
OMC Fouree Stoner
OMC Fouror bkingg
1011
1001
16
singlo
Time wtgphutrucce
Figure 8: Loglog plot of the relative L approximation етгоr versus computation time for comput- ing the kerne summations ( with the negative distance and thin plate spline kernel and different methods and datasets. MNIST and FashionMNIST are reduced to dimension d = 20 and = 30 via PCA. We run each method 10 times. The shaded area indicates the standard deviation of the error. For the slicing method, we use P 10 2 slices for k 0,..,
34

---

Published as a conference paper at ICLR 2025
MNIST
32-4
ResEINe of
CPC FFF-20 slicing
CMC RFT 1FF-20 stiing
CMC REF-LO sliiin
Slicigg
dicing
cing
20
Time to single tnd CPU
Timee tegll CPU
Thmse 50 on dormoede CPU
Laplace bernel in MMST (d 7841
negalive distance on MNBT 1d 784)
slicin GMC Fournr Stiiing
10
Sertin
GMC Sortin Steie
MaShionMNIST
Gauss kernel on FashicreMNIST (dim 784)
Matern (urm 3.5) on FashiorMasT (s 784)
Mater (um 1:51 FashionMMST (d 784
CHE
12
RFF
Ri
22-4
QMC Fourer Siiing
GMC Fastersiiing
QMC Fouter siiiig
FashionMNiST
regative FachienMMST
10.
22-4
OMC so Saring sticoo
og Sicke
Time digglo www.lleetaddddd......oo... CPU
Time wtgle-traddedd CPU
Figure 9: Loglog plot of the relative L' approximation eπrо versus computation time for computing
the kernel sunmations (I) with different kernels and methods on the MNIST and FashionMNIST indicates the
dataset without dimension reduction. We each method 10 times. The shaded standard deviation the enror. For the slicing method, we use P 10 2 slices for k: = O..., 9. In order obtain similar computation times, we RFF and ORF with 2P features.
35

---

Published as a conference paper at ICLR 2025
10P1
Reoddve 12
QMC Fouter Siiing
QMC Surin
OeT Fooes siring
skioo
Faster cine
The NOEA Galema RXX 9030 GPJ
1hrano en Wili Galara RX 0030 GFJ
Time Do Galarrs XIX 1060 GPU
Leplac cene with 1 193033.4-13
negulve Docae Wirh 20- DE0000, c1-30
22-1
2314
TANAAWE Stee
RFF
Slicig
Serting PyRe0zi
Py
10
12?
Figure 10: Loglog plot of the relative L' approximation error versus GPU computation time for computing the kernel summations (I) with different kernels large dataset (M N 960000). We run each method 10 times. The shaded area indicates the standard deviation of the error. For the slicing method, times, we run RFF and ORF with D = 4P features. Since PyKeOps computes the exact kernel sum, the computation time indicated by a vertical line.
P 10 2 slices for k O,.... 9.
In order obtain similar computation
Objective value vs time for the MMD flow on CIFAR10
slicing
OMC-siking
10
objcctive value
4x10-
xx10
2x10-
30>
10"
computation time
Figure 11: Plot of the objective value F, (x) for the MMD gradient flow the CIFARIO dataset an NVIDIA GeForce RTX 4090 GPU. Note that computing this flow with exact kernel summations computationally intractable. The plot omits the first 2000 steps of the gradient flow improve the scaling the plot.
36


Reviews:
[{"review_id": "yME5dpnSPU", "reviewer": "Reviewer_WPrq", "summary": "In order to reduce the compute time of kernel summations in large dimension, this paper studies 1D slicing when computing kernels on R^d of the form (x,x’) -> F(\\| x – x’ \\|). The paper is directly built on a work of Hertrich (2024), that introduced the framework of the generalized Riemann-Liouville fractional integral transform to characterize the relationship between the sliced approximation and its target. In a first contribution, a novel error bound is proven with a characterization of the variance of the estimator when Slices are uniformly distributed. In a second contribution, the authors discuss quasi-Monte-Carlo designs (spherical, Sobol) to choose the slices and improve the error rates. The paper is completed by numerical simulation on summations of kernels on artificial data and on MNIST data. While QMC-slicing is interesting in terms of approximation error (convergence rate) their design is expensive and prohibitive for d>= 3.", "questions": "1° Can you comment on the interest of slicing in Optimal Transport versus in kernel summation ? \n2° It would be interesting to describe the differences point by point of this work and those of Hertrich (2024).\n3° Could you clearly describe the analytical complexity in time of each QMC scheme ?\n4° Can you find and experiment an example in Machine learning where the sliced kernels can make a difference ? \n5° is it possible to run the same tests on larger dimensional data ?\n\nAfter rebuttal, I've increased my score (most of the concerns have been addressed except those concerning the experimental results.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Although slicing has been thoroughly explored in reducing the complexity in time of optimal transport, slicing the computation of kernels is a novel topic, very rencelty introduced.\n- The novel error bound differs from Hertrich (2024, SIAM Journal on Mathematics of Data Science and arxiv) with a characterization of the vairance of the slicing estimator, which confirms the rate in O(1/sqrt(P)).\n- The most interesting and promising part from my point of view is related to the exploration of quadrature rules to construct the sequence of slices and not sampling it. The proposed scheme (approximation of a spherical design) comes with a bound on the approximation error.", "weaknesses": "- About the motivation:The number of data is usually considered as the most emblematic issue with kernels in Machine Learning and existing approximation schemes aim at reducing the compute time involved by operations with the Gram matrix as well as the complexity in memory. The authors motivate their work by the case when the dimension of data is large. Contrary to optimal transport problems that are defined in a variational way,  computing sums of kernels at the era of distributed computing is not a crucial obstacle.\n\n- Limited contribution given previous works:\nThe most important weakness of the paper is its incrementality with respect to the work of Hertrich (2024) well cited in the paper. \nThe revisited error bound in the case of the Monte-Carlo estimate does not seem sufficient to make a difference so the new component of the work deals with the proposition of the quasi-Monte-Carlo design for the slicing method.\n- No applicabe scale for the spherical deisgn: Unfortunately, the spherical design studied in depth here cannot be reasonably computed for  realistic values of d. \nSo finally either Monte-Carlo estimation or classic Sobol sequences have to be used in practise (dataset MNIST) which again limits dratsically the novelty of the paper.\n\nIn terms of experiments, I would not consider the dimension of MNIST data as a computational obstacle and would definitevely expect higher dimensional data for a conference like ICLR.\n\nMoreover it would have been interesting to apply this new approximated way to compute sums of kernels in a statistical kernel-based test or in an online learning algorithm of a kernel machine where this approximation can make sense. Studying then the convergence of the obtained estimator would have been entirely new.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "e3mRqHr5ox", "reviewer": "Reviewer_EyTQ", "summary": "This paper proposes a quasi-Monte Carlo (QMC) slicing approach for fast radial kernel summation, which is commonly required in kernel methods across machine learning. By using QMC sequences and spherical quadrature rules, the authors derive error bounds for various kernels, such as the Gauss, Laplace, and Matérn kernels, and show that QMC slicing significantly improves computation time and accuracy. Numerical experiments confirm that QMC slicing outperforms non-QMC slicing and other approximation methods, particularly in lower-dimensional cases, while maintaining favorable error rates.", "questions": "1.How scalable is the QMC slicing method for very high-dimensional data (e.g., 100+ dimensions)?\n2.Does the QMC slicing approach generalize to non-standard kernels often used in machine learning?\n3.How do theoretical error bounds hold up on real-world datasets with noise and variability?\n4.Can the authors compare the computational efficiency of QMC slicing to Random Fourier Features and other fast summation methods?\n5.What are the practical trade-offs in terms of error convergence rate with QMC slicing?\n6.Are there specific use cases or domains where this approach is expected to outperform current alternatives?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1.The QMC slicing method introduces an approach to fast kernel summation.\n2.The methodology is well-documented, with explanations of error bounds and smoothness results.\n3.QMC slicing has the potential to improve efficiency in kernel-based methods.", "weaknesses": "1.The adaptability of the theoretical error bounds to various types of data and kernels remains unclear.\n2.While effective, the QMC slicing method lacks clear differentiation from existing fast summation methods, such as Random Fourier Features (RFF).\n3.Experiments are focused on limited, synthetic datasets, raising concerns about the method's generalizability.\n4.Some parts of the theory section are overly condensed, especially around error bounds and smoothness assumptions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yizfwZwhzc", "reviewer": "Reviewer_T9eK", "summary": "The paper proposes a slicing based fast calculation of the sum of radial kernels. The authors analyses bounds on slicing error, and ensure a better error rate based on the smoothness result.", "questions": "- How is the optimization problem (10) is solved in practice? How large is this additional computational cost? Figure 2 contains it? In particular, when P and d is large, it seems expensive.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The theoretical analysis in the paper seemingly technically sound, though I couldn't follow the proofs.", "weaknesses": "- Overall, the paper is difficult to follow those who are not familiar with the topic. Readability can be improved. \n\n- Empirical evaluation is only on purely kernel sum calculations. Showing benefit on actual computation of some learning model would have been convincing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lXWAkDTGjY", "reviewer": "Reviewer_Edek", "summary": "This paper derived the error bounds of slicing, a fast approximation algorithm for computing kernel sums, for various kernels. The paper also incorporated quasi-Monte Carlo (QMC) sequences on the sphere into the slicing method to improve its error rate. The paper conducted experiments to demonstrate that the QMC-slicing approach outperforms the non-QMC slicing method and other baselines like Random Fourier Features (RFF) in terms of error rate.", "questions": "1. The authors mentioned that QMC-slicing exhibits a smaller advantage in high dimensions (d=200), but it still outperforms other baselines. What if the dimensionality is even higher, e.g. d=10,000? Does QMC-slicing still perform better than slicing and (QMC-)RFF?\n2. In this paper, only Gauss, Laplace, Matern and negative distance kernels are considered. How would QMC-slicing perform in the case of other kernels, e.g. sigmoid kernel?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The authors proposed the QMC-slicing method and demonstrated that it has better error rates than the original slicing method. The QMC-slicing method exhibits a better expected error bound than RFF while maintaining a similar time complexity; this demonstrates its effectiveness and potential as a fast kernel summation algorithm.\n2. The authors provided a comprehensive analysis of the slicing error for the negative distance kernel, the thin plate spline, the Laplace kernel and the Gauss kernel. The derived error bounds could be useful for future research on slicing and fast kernel summation.\n3. Overall, the paper is well-written and clear to understand. Background knowledge on fast kernel summation, slicing and QMC are well-explained. Details of derivation, e.g. the applicability of the QMC methods for slicing, have been sufficiently provided.", "weaknesses": "1. The authors considered several QMC sequences in the experiments, and observed difference in slicing error between these different sequences. It would be helpful if the authors could further analyze the results and discuss more on how the choice of QMC sequences affects slicing error, under what circumstances would some QMC sequences outperform others, etc.\n2. As mentioned in Section 5, the observed error rates of QMC-slicing are significantly better than the theory. It would be helpful if the authors discuss more about this phenomenon and the reason behind it.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)