Fast Summation of Radial Kernels via QMC Slicing

## Abstract

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.

## 1 Introduction

We consider fast algorithms for computing the kernel sums

$$ $s_{m}=\sum_{n=1}^{N}w_{n}K(x_{n},y_{m}),\quad\text{for all}\quad m=1,...,M,$ (1) $$

where $x_{n},y_{m}\in\mathbb{R}^{d}$ and $w_{n}\in\mathbb{R}$ for $n=1,...,N$, $m=1,...,M$ and $K\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ is a radial kernel, i.e.,
$K(x,y)=F(\|x-y\|)$ for the Euclidean norm $\|\cdot\|$ and some $F\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$.
This summation problem appears in most kernel methods including kernel density estimation (Parzen, ; Rosenblatt, ), classification via support vector machines (Steinwart & Christmann, ), dimensionality reduction with kernelized principal component analysis (Schölkopf & Smola, ; Shawe-Taylor & Cristianini, ), distance measures on the space of probability measures like maximum mean discrepancies or the energy distance (Gretton et al., ; Székely, ), corresponding gradient flows (Arbel et al., ; Galashov et al., ; Hagemann et al., ; Hertrich et al., ; Kolouri et al., ), and methods for Bayesian inference like Stein variational gradient descent (Liu & Wang, ).
Computing equation 1 exactly for all $m=1,...,M$ has complexity $\mathcal{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., ; Potts et al., ), which approximates the kernel by a truncated Fourier series and reformulates equation 1 using the fast Fourier transform on non-equispaced data (Beylkin, ; Dutt & Rokhlin, ). We provide a short overview in Appendix I.
This kind of methods usually provides a computational complexity of $\mathcal{O}(M+N+N_{\mathrm{ft}}\log(N_{\mathrm{ft}}))$, where $N_{\mathrm{ft}}$ is the number of relevant Fourier coefficients, and admits very fast error rates for $N_{\mathrm{ft}}\to\infty$ (even exponential if the kernel is sufficiently smooth). However, the number $N_{\mathrm{ft}}$ 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 () proposed random Fourier features (RFF). They represent a positive definite kernel via Bochner’s theorem (Bochner, ) as the Fourier transform of a non-negative measure. Sampling randomly from this measure at $D\in\mathbb{N}$ points (features) leads to an approximation algorithm with computational complexity $\mathcal{O}(D(N+M))$ independent of the dimension $d$. However, the error decays only with rate $\mathcal{O}(1/\sqrt{D})$, which can be limiting if a high accuracy is required.
Moreover, RFF are limited to positive definite kernels and do not apply to other kernels, like the negative distance kernel $K(x,y)=-\|x-y\|$, which has applications, e.g., within the energy distance (Székely, ) that is used for defining a distance on the space of probability measures.

A related approach is slicing (Hertrich, ), which represents the kernel sum equation 1 as an expectation of one-dimensional kernel sums of the randomly projected data points with a different kernel.
Discretizing the expectation by sampling at $P$ random projections, the kernel sums equation 1 can be approximated 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 $\mathcal{O}(P(N+M+N_{\mathrm{ft}}\log N_{\mathrm{ft}}))$, where the expected error can be bounded by $\mathcal{O}(1/\sqrt{P})$.
For positive definite kernels, a close relation between RFF and slicing was established by Rux et al. (), see the short overview in Appendix G.
One advantage of slicing is the applicability to kernels that are not positive definite.

A way to improve the $\mathcal{O}(1/\sqrt{P})$ 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., ). Note that there also exist QMC approaches for RFF (Avron et al., ). However, they depend on the restrictive assumption that the measure from Bochner’s theorem decouples over the dimension, which is true for the Gauss and $L^{1}$-Laplace(^1^11In literature, there exist two versions of the Laplace kernel $K(x,y)=\exp(-\alpha\|x-y\|_{1})$ and $K(x,y)=\exp(-\alpha\|x-y\|)$, which differ in the used norm. Since our analysis focuses on radial kernels, we will only consider the second version in the rest of this paper.) kernel, but false for most other common kernels like the Laplace, Matérn or negative distance kernel.

##### Contributions

Our contributions to the slicing approach for fast kernel summation are as follows:

- -
We derive bounds on the slicing error for various kernels in Theorem 1. In particular, 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 $\mathcal{O}(1/\sqrt{P})$. To ensure the applicability of the QMC error bounds, we prove certain smoothness results for the function which maps a direction $\xi\in\mathbb{S}^{d-1}$ to the corresponding one-dimensional kernel in Theorem 3. The improved error rates are outlined in Corollary 4.
- -
We conduct extensive numerical experiments on standard test datasets for several kernels and different popular QMC sequences, and demonstrate that our QMC-slicing approach significantly outperforms the non-QMC slicing method as well as (QMC-)RFF. While the advantage of QMC-slicing is most significant in dimensions $d\leq 100$, it also performs better in higher dimensions.

##### Outline

In Section 2, we 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 5. Additional proofs, plots and evaluations are contained in the appendix. The code for the numerical examples is provided online(^2^22available at [https://github.com/johertrich/fastsum_QMC_slicing](https://github.com/johertrich/fastsum_QMC_slicing)).

### Related Work

##### Low-Dimensional Kernel Summation

Fast summation algorithms have been extensively studied in the literature.
They include fast kernel summations based on (non-)equispaced fast Fourier transforms (Greengard et al., ; Kunis et al., ; Potts et al., ), fast multipole methods (Beatson & Newsam, ; Greengard & Rokhlin, ; Lee & Gray, ; Yarvin & Rokhlin, ), tree-based methods (March et al., ; ) or H- and mosaic-skeleton matrices (Hackbusch, ; Minden et al., ; Tyrtyshnikov, ).
For the Gauss kernel, the fast Gauss transform was proposed by Greengard & Strain () and improved by Yang et al. (; ).
More general fast kernel transforms were considered by Ryan et al. ().

##### QMC and Quadrature on Spheres

QMC designs on spheres were studied by Brauchart et al. (). Here, the quadrature points optimizing the worst-case error in certain Sobolev spaces are given by spherical $t$-designs, which integrate all polynomials up to degree $t$ on the sphere exactly (Delsarte et al., ; Bannai & Bannai, ).
The construction of spherical designs is highly non-trivial and intractable in high dimensions.
For $\mathbb{S}^{2}$ and $\mathbb{S}^{3}$, several examples were computed numerically by Gräf & Potts () and Womersley ().
Gräf et al. () related quadrature rules on the sphere with halftoning problems.
The specific problem of QMC rules for the three-dimensional sliced Wasserstein distance was considered in Nguyen et al. ().

##### Random Fourier Features

Random Fourier features (RFFs) were proposed by Rahimi & Recht () and were further analyzed in several papers Bach (); Hashemi et al. (); Li et al. ().
To improve the error rates, Avron et al. () proposed a 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. () and Munkhoeva et al. (). Yu et al. () proposed orthogonal Fourier features.
In a very recent preprint, Belhadji et al. () derive an explicit quadrature rule in the Fourier space for efficient summations of the Gauss kernel.

## 2 Slicing of Radial Kernels

Let $K\colon\mathbb{R}^{d}\times\mathbb{R}^{d}$ be a radial kernel of the form $K(x,y)=F(\|x-y\|)$ for some basis function $F\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$.
Throughout this paper, we will assume that $K$ has the form

$$ $K(x,y)=\mathbb{E}_{\xi\sim\mathcal{U}_{\mathbb{S}^{d-1}}}[\mathrm{k}(\langle \xi,x\rangle,\langle\xi,y)\rangle]$ $$

for some one-dimensional radial kernel $\mathrm{k}\colon\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ with basis function $f\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$, where we surpress 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(\|x\|)=\mathbb{E}_{\xi\sim\mathcal{U}_{\mathbb{S}^{d-1}}}[f(|\langle\xi,x \rangle|)].$ (2) $$

A pair $(F,f)$ of basis functions in $L^{\infty}_{\mathrm{loc}}(\mathbb{R}_{\geq 0})$ fulfills this relation if and only if $F$ is the *generalized Riemann–Liouville fractional integral* transform given by

$$ $F(t)=\frac{2\Gamma(\frac{d}{2})}{\sqrt{\pi}\Gamma(\frac{d-1}{2})}\int_{0}^{1}f (ts)(1-s^{2})^{\frac{d-3}{2}}\mathrm{d}s,$ (3) $$

for $2\leq d\in\mathbb{N}$, see (Hertrich, , Prop 2 and Rubin, ).
In order to find the one-dimensional basis function $f$ for a given $F$, we have to invert the transform equation 3.
This can be done explicitly if

- i)
$F$ is analytic on $\mathbb{R}_{\geq 0}$, i.e., there exists $\left(a_{n}\right)_{n\in\mathbb{N}}\subset\mathbb{R}$ such that $F(x)=\sum_{n=0}^{\infty}a_{n}x^{n}$ for all $x\geq 0$, or
- ii)
$F(\|\cdot\|)\colon\mathbb{R}^{d}\to\mathbb{R}$ is continuous, bounded and positive definite, i.e., for all $N\in\mathbb{N}$, all pairwise distinct $x_{j}\in\mathbb{R}^{d}$ and all $a_{j}\in\mathbb{R}$ for $j=1,...,N$ it holds that
$\sum_{j,k=1}^{N}a_{j}a_{k}F(\|x_{j}-x_{k}\|)\geq 0$,

see (Hertrich, , Thm 3 and Rux et al., , Cor 4.11).
We include a list of pairs $(F,f)$ fulfilling equation 2 in Appendix A. In particular, it includes the basis functions of Gauss, Laplace and Matérn kernels. Moreover, $f$ 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 equation 1 efficiently, we approximate $F(\|\cdot\|)$ by projections and the one-dimensional basis function $f$, i.e., we aim to find directions $\xi_{1},...,\xi_{P}\in\mathbb{S}^{d-1}$ such that

$$ $F(\|x\|)\approx\frac{1}{P}\sum_{p=1}^{P}f(|\langle\xi_{p},x\rangle|)\quad\text {for all}\quad x\in\mathbb{R}^{d}.$ (4) $$

Then, the kernel sums equation 1 read as

$$ $s_{m}=\sum_{n=1}^{N}w_{n}K(x_{n},y_{m})=\sum_{n=1}^{N}w_{n}F(\|x_{n}-y_{m}\|) \approx\frac{1}{P}\sum_{p=1}^{P}\sum_{n=1}^{N}w_{n}f(|\langle\xi_{p},x_{n}-y_{ m}\rangle|).$ (5) $$

For computing the one-dimensional sums $\sum_{n=1}^{N}w_{n}f(|\langle\xi_{p},x_{n}-y_{m}\rangle|)$ for all $m=1,...,M$, there exists algorithms with complexity $\mathcal{O}(M+N)$ or $\mathcal{O}((M+N)\log(M+N))$ in literature. These include fast summations based on non-equispaced Fourier transforms (Kunis et al., ; Potts et al., ), fast multipole methods (Greengard & Rokhlin, ) or sorting algorithms (Hertrich et al., ).
In particular, we can approximate the vector $s=(s_{1},...,s_{M})$ via equation 5 with a complexity of $\mathcal{O}(P(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 equation 4.
To this end, we assume that the directions $\xi_{1},...,\xi_{P}$ are iid samples from the uniform distribution on the sphere.
Then, we exactly compute the variance

$$ $\mathbb{V}_{d}[f](x):=\mathbb{E}_{\xi\sim\mathcal{U}_{\mathbb{S}^{d-1}}}\left[ \left(f(|\langle\xi,x\rangle|)-F(\|x\|)\right)^{2}\right],$ (6) $$

which bounds the mean squared error through the Bienaymé’s identity as

$$ $\mathbb{E}_{\xi_{1},...,\xi_{P}\sim\mathcal{U}_{\mathbb{S}^{d-1}}}\left[\left( \frac{1}{P}\sum_{p=1}^{P}f(|\langle\xi_{p},x\rangle|)-F(\|x\|)\right)^{2} \right]=\frac{\mathbb{V}_{d}[f](x)}{P}.$ $$

In particular, our results show relative error bounds of the negative distance kernel $K(x,y)=-\|x-y\|$ and the thin plate spline (except around $\|x\|=1$), which do not depend on the dimension $d$. The proof is given in Appendix C.

###### Theorem 1 .

Let $F\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$ and $f\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$ fulfill the slicing relation equation 2.

- i)
If $F(\|\cdot\|)$ is continuous and positive definite on $\mathbb{R}^{d}$, then $\mathbb{V}_{d}[f](x)\leq f(0)^{2}-F(\|x\|)^{2}$.
- ii)
For the generalized Riesz kernel $F(\|x\|)=-\|x\|^{r}$ with $r>0$, we have
$\displaystyle\mathbb{V}_{d}[f](x)$
$\displaystyle=\left(\frac{\sqrt{\pi}\Gamma(r+\frac{1}{2})}{\Gamma(\frac{r+1}{2
})^{2}}\frac{\Gamma(\frac{d+r}{2})^{2}}{\Gamma(\frac{d}{2})\Gamma(r+\frac{d}{2
})}-1\right)\left(F(\|x\|)\right)^{2}$
$\displaystyle\leq\left(\frac{\sqrt{\pi}\Gamma(r+\frac{1}{2})}{\Gamma(\frac{r+1
}{2})^{2}}-1\right)\left(F(\|x\|)\right)^{2}.$
(7)
- iii)
For the thin plate spline $F(\|x\|)=\|x\|^{2}\log(\|x\|)$, we have
$\mathbb{V}_{d}[f](x)=\left(2+\frac{c_{1}}{\log(\|x\|)}+\frac{c_{2}+\mathcal{O}
(d^{-1}\log d)}{\log(\|x\|)^{2}}\right)\left(1+\tfrac{2}{d}\right)\left(F(\|x
\|)\right)^{2},$
with the numerical constants
$\displaystyle c_{1}$
$\displaystyle:=-3\gamma-\log(8)+8\approx 4.189,\text{ and}$
(8)
$\displaystyle c_{2}$
$\displaystyle:=\frac{3}{4}\left(\psi^{(1)}\left(\tfrac{5}{2}\right)+\left(\psi
^{(0)}\left(\tfrac{5}{2}\right)+\log(2)\right)^{2}\right)+2(\gamma+\log(2)-2)^
{2}\approx 2.895.$
- iv)
For $F(\|x\|)=\sum_{n=0}^{\infty}a_{n}\|x\|^{n}$ and $d\geq 3$ odd, we have
$\mathbb{V}_{d}[f](x)=\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}\left(\Pi_{i=1}^{
\frac{d-1}{2}}\left(1+\frac{k(n-k)}{(2i-1)(n+2i-1)}\right)-1\right)a_{k}a_{n-k
}\right)\|x\|^{n}.$
In particular,
-
for the Laplace kernel $F(\|x\|)=\exp(-\alpha\|x\|)$, we have
$\mathbb{V}_{3}[f](x)=\frac{1}{4\alpha\|x\|}\left(1-(2(\alpha\|x\|)^{2}+2\alpha
\|x\|+1)F(\|x\|)^{2}\right),$
-
and for the Gauss kernel $F(\|x\|)=\exp(-\|x\|^{2}/(2\sigma^{2}))$, we have
$\mathbb{V}_{3}[f](x)=\frac{\sigma^{2}}{2\|x\|^{2}}\left(1-\left(\frac{\|x\|^{4
}}{2\sigma^{4}}+\frac{\|x\|^{2}}{\sigma^{2}}+1\right)F(\|x\|)^{2}\right).$

Note that some weaker error bounds for the Gaussian and Riesz kernel were also shown in Hertrich ().
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 obtain an absolute error bound i) independent of $\|x\|$. For the Riesz kernel, we have a relative error bound in ii). For the thin plate spline, iii) behaves like a relative bound for $\|x\|\to\infty$ and $\|x\|\to 0$, but as a constant around $\|x\|=1$, which is a zero of $F$. For the Laplace and Gauss kernel, the dependence on $\|x\|$ changes drastically between $\|x\|\to\infty$ and $\|x\|\to 0$. In fact, $\mathbb{V}_{3}[f](x)$ is monotonically increasing in $\|x\|$ with global upper bound $1/(4\alpha)$, and converges quadratically in $\|x\|$ to zero for $\|x\|\to 0$. For the case $d>3$, we conjecture a similar behavior, see Appendix D for the discussion.

###### Theorem 1 .

## 3 Quasi-Monte Carlo Slicing

For directions drawn independently from the uniform measure $\mathcal{U}_{\mathbb{S}^{d-1}}$ on the sphere $\mathbb{S}^{d-1}$, our experiments from the numerical part suggest that the rate $\mathcal{O}(1/\sqrt{P})$ from Theorem 1 is optimal. As a remedy, we employ quadrature and quasi-Monte Carlo designs on the sphere for improving these error rates. To this end, we first revisit the corresponding literature in Subsection 3.1. Afterwards, we apply these results for our slicing method in Subsection 3.2.

### 3.1 Quasi-Monte Carlo Methods on the Sphere

Let ${\bm{\xi}}^{P}=(\xi_{1}^{P},...,\xi_{P}^{P})\in(\mathbb{S}^{d-1})^{P}$ for $P\in\mathbb{N}$. In the following, we aim to construct ${\bm{\xi}}^{P}$ such that the worst case error in a certain Sobolev space is asymptotically optimal.
The definition of the Sobolev space $H^{s}(\mathbb{S}^{d-1})$ is given in Appendix B.

###### Definition 2 .

A sequence $({\bm{\xi}}^{P})_{P}$ with $P\to\infty$ is called a sequence of QMC designs for $H^{s}(\mathbb{S}^{d-1})$ if there exists some $c(s,d)>0$ independent of $P$ such that the worst case error

$$ $\sup_{\begin{subarray}{c}f\in H^{s}(\mathbb{S}^{d-1})\\ \|f\|_{H^{s}(\mathbb{S}^{d-1})}\leq 1\end{subarray}}\left|\frac{1}{|\mathbb{S} ^{d-1}|}\int_{\mathbb{S}^{d-1}}f(\xi)\mathrm{d}\xi-\frac{1}{P}\sum_{p=1}^{P}f( \xi^{P}_{p})\right|\leq\frac{c(s,d)}{P^{s/(d-1)}}\in\mathcal{O}(P^{-s/(d-1)}).$ (9) $$

It was proven by Hesse () that the error rate $\mathcal{O}(P^{-s/(d-1)})$ is optimal, see also Brauchart et al. ().
For $s>\frac{d-1}{2}$, the existence of sequences of QMC designs is ensured by so-called spherical designs.
More precisely, ${\bm{\xi}}^{P}$ is called a spherical $t$-design
if the quadrature at these points integrates all polynomials of degree $t\in\mathbb{N}$ exactly, i.e., if it holds

$$ $\frac{1}{|\mathbb{S}^{d-1}|}\int_{\mathbb{S}^{d-1}}\psi(\xi)\mathrm{d}\xi= \frac{1}{P}\sum_{p=1}^{P}\psi(\xi^{P}_{p}),\quad\text{for all polynomials $ \psi$ of degree }\leq t.$ $$

It can be shown that for any $t$ there exists a spherical $t$-design with $P=\mathcal{O}(t^{d-1})$ points, see Bondarenko et al. ().
It was proven in (Brauchart & Hesse, , Cor 3.6) that a sequence of such spherical $t$-designs is a sequence of QMC designs, see also (Brauchart et al., , Sect 1) for a summary.

Unfortunately, the construction of spherical $t$-designs in arbitrary dimension is numerically intractable. Instead, many QMC methods rely on low-discrepancy point sets. It was shown in (Brauchart et al., , Thm 14) that a sequence $\bm{\xi}^{P}$ that minimizes the sum of powers of Euclidean distances

$$ $\mathcal{E}({\bm{\xi}}^{P})\coloneqq-\sum_{p,q=1}^{P}\|\xi^{P}_{p}-\xi^{P}_{q} \|^{2s-d-1}$ (10) $$

is a QMC design for $H^{s}(\mathbb{S}^{d-1})$ for $s\in(\frac{d-1}{2},\frac{d+1}{2})$.
In the numerics, we will consider $s=\frac{d}{2}$, so that we get a QMC design for $H^{d/2}(\mathbb{S}^{d-1})$.
Note that, up to a constant, equation 10 coincides with the maximum mean discrepancy with the Riesz kernel $K(x,y)=-\|x-y\|^{2s-d-1}$ between the probability measures $\frac{1}{P}\sum_{i=1}^{P}\delta_{\xi^{P}_{p}}$ and $\mathcal{U}_{\mathbb{S}^{d-1}}$, which is also known as energy distance Székely (). Furthermore, one can easily transform the QMC sequence into an unbiased estimator in equation 4, see Appendix E.

###### Definition 2 .

### 3.2 Smoothness of One-Dimensional Basis Functions

In order to apply the above theorems for our approximation equation 2, we need to ensure that for any $x\in\mathbb{R}^{d}$ the spherical function $g_{x}:\mathbb{S}^{d-1}\to\mathbb{R}$ given by

$$ $\xi\mapsto g_{x}(\xi)=f(|\langle\xi,x\rangle|)$ (11) $$

is sufficiently smooth on $\mathbb{S}^{d-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 $g_{x}$. In the special case $s=0$, this relates to the variance equation 6 by the formula
$\|g_{x}\|_{H^{0}\left(\mathbb{S}^{d-1}\right)}^{2}=|\mathbb{S}^{d-1}|(\mathbb{
V}_{d}[f](x)+F(\|x\|)^{2})$ for $(f,F)$ fulfilling equation 2.

###### Theorem 3 .

Let $x\in\mathbb{R}^{d}$ with $x\neq 0$.
For the Gauss, Riesz and Matérn kernel, the following smoothness results hold true:

- i)
For
$F(t)=\exp(-\frac{t^{2}}{2\sigma^{2}})$,
we have $g_{x}\in H^{s}(\mathbb{S}^{d-1})$ for all $s\geq 0$.
- ii)
For $F(t)=t^{r}$ with $t\geq 0$ and $r>-1$, we have $g_{x}\in H^{s}(\mathbb{S}^{d-1})$ if and only if $s<r+\frac{1}{2}$.
- iii)
For $F(t)=\frac{2^{1-\nu}}{\Gamma(\nu)}(\tfrac{\sqrt{2\nu}}{\beta}t)^{\nu}K_{\nu}(
\tfrac{\sqrt{2\nu}}{\beta}t)$, $t\geq 0$,
we have $g_{x}\in H^{s}(\mathbb{S}^{d-1})$ if $s<2\nu+\tfrac{1}{2}$.

Note that the theorem also includes the Laplace kernel, which is the Matérn kernel for $\nu=\frac{1}{2}$.
Combining 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\in\mathbb{N}$ and $s>\frac{d-1}{2}$. Then there exists a constant $c(s,d)$ and a sequence $({\bm{\xi}}^{P})_{P}$ with $P\to\infty$
such that for the Gauss and Matérn kernel with basis functions $F(t)=\exp(-\frac{t^{2}}{2\sigma^{2}})$ and $F(t)=\frac{2^{1-\nu}}{\Gamma(\nu)}(\tfrac{\sqrt{2\nu}}{\beta}x)^{\nu}K_{\nu}(
\tfrac{\sqrt{2\nu}}{\beta}t)$ with $\nu>\frac{2s-1}{4}$, respectively, it holds that

$$ $\sup_{x\in\mathbb{R}^{d}}\left|F(\|x\|)-\frac{1}{P}\sum_{p=1}^{P}f(|\langle\xi _{p},x\rangle|)\right|\leq\frac{c(s,d)}{P^{\frac{s}{d-1}}}.$ $$

For $s\in(\frac{d-1}{2},\frac{d+1}{2})$, such ${\bm{\xi}}^{P}$ are given by minimizers of equation 10.

###### Theorem 3 .

###### Corollary 4 .

## 4 Numerical Examples

In the following, we evaluate the kernel approximation with QMC slicing for several QMC sequences and compare our results with random Fourier feature-based (RFF, Rahimi & Recht, ) methods. We implement the comparison in Julia and Python, and provide the code online(^3^33available at [https://github.com/johertrich/fastsum_QMC_slicing](https://github.com/johertrich/fastsum_QMC_slicing)). In Subsection 4.1, we describe the used QMC sequences, RFF-methods and kernels. Afterwards, in Subsection 4.2, we numerically evaluate the approximation error in equation 4. Finally, we apply our approximation for fast kernel summations in Subsection 4.3.
We include additional numerical examples in Appendix J.

### 4.1 QMC Sequences and Kernels

##### QMC Sequences

Beside standard slicing where the projections ${\bm{\xi}}^{P}$ are drawn iid from the uniform distribution on $\mathbb{S}^{d-1}$, we consider the following sequences:

- -
Sobol Sequence: Two commonly used QMC sequences on $[0,1]^{d}$ are Sobol (Sobol’, ) and Halton (Halton, ) sequences. They can be transformed to QMC sequences for the normal distribution by applying the inverse cumulative density function along each dimension of the sequence, which was used for deriving a QMC method for random Fourier features Avron et al. (). To obtain a QMC sequence on the sphere, Nguyen et al. () proposed to project the QMC sequence for the multivariate normal distribution onto the sphere by the transformation $\xi=\theta/\|\theta\|$, see also Beltrán et al. (). To generate the original Sobol sequence on $[0,1]^{d}$, we use the implementation from SciPy (Virtanen et al., ) in Python and the Sobol.jl package in Julia. Our numerical experiments suggest that the Sobol sequence lead to slightly better results than the Halton sequence. Therefore, we omit the Halton sequence in our comparison.
- -
Orthogonal: Even though this is technically not a QMC sequence, we adapt the approach of orthogonal features (Yu et al., ) for slicing and generate directions $\xi$ as follows: We generate $\lceil\frac{P}{d}\rceil$ orthogonal matrices from the uniform distribution on $\mathrm{O}(d)$ (this can be done by taking the Q-factor of the QR decomposition applied on a matrix with standard normally distributed entries). Together, these matrices have $d\lceil\frac{P}{d}\rceil$ columns from which we choose ${\bm{\xi}}^{P}$ to be the first $P$ of those.
- -
Distance: In Section 3.1, we considered the QMC sequence ${\bm{\xi}}^{P}$ for $H^{\frac{d}{2}}(\mathbb{S}^{d-1})$ that is a minimizer of
$\mathcal{E}({\bm{\xi}}^{P})=-\sum_{p,q=1}^{P}\|\xi_{p}-\xi_{q}\|$, see equation 10.
In our application, we have the additional symmetry that $f(|\langle x,\xi\rangle|)=f(|\langle x,-\xi\rangle|)$. Therefore, we construct a QMC sequence ${\bm{\xi}}^{P}$ by minimizing the interaction energy functional $\mathcal{E}_{\mathrm{sym}}({\bm{\xi}}^{P})\coloneqq\mathcal{E}(({\bm{\xi}}^{P}
,-{\bm{\xi}}^{P}))$.
We do this numerically with the Adam optimizer (Kingma & Ba, ) and the PyKeops package (Charlier et al., ), which takes from a couple of seconds (for $d=3$) up to one hour (for $d=200$ and $P\approx 5000$) on an NVIDIA RTX 4090 GPU.
In Appendix H, we show that if $P\leq d$, the orthogonal points from above minimize $\mathcal{E}_{\mathrm{sym}}$, so this approach differs only if $P>d$.
- -
Spherical Design: For $d=3$, several spherical $t$-designs on the $\mathbb{S}^{2}$ were computed by Gräf & Potts () up to $t\leq 1000$ and $P\leq 1002000$ and are available online(^4^44[https://www-user.tu-chemnitz.de/~potts/workgroup/graef/quadrature/index.php.en](https://www-user.tu-chemnitz.de/~potts/workgroup/graef/quadrature/index.php.en)). Spherical $t$-designs for $\mathbb{S}^{3}$ were computed by Womersley () up to $t\leq 31$ and $P\leq 3642$.
Unfortunately, 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, ): see Appendix G for a description.
- -
Orthogonal Random Features (ORF, Yu et al., ): The directions of the RFF features are chosen in the same way as explained above for the orthogonal slicing.
- -
QMC-Random Fourier Features (Sobol RFF, Avron et al., ): For the Gauss kernel, we also compare with QMC random Fourier features, which are only applicable for kernels where the Fourier transform decouples as 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 $\nu=p+\frac{1}{2}$ for $p\in\{1,3\}$) and the negative distance kernel (Riesz kernel with $r=1$), see Table 2 in the appendix for the pairs $(f,F)$. The parameters $\sigma$, $\alpha$ and $\beta$ are chosen by the median rule (see, e.g., Garreau et al., for an overview and history). That is, we choose $\sigma=\beta=\frac{1}{\alpha}=\gamma m$, where $m$ is the median of all considered input norms $\|x\|$ of the basis functions $F$ and $\gamma$ is some scaling factor which we set to $\gamma\in\{\frac{1}{2},1,2\}$.

### 4.2 Numerical Evaluation of the Slicing Error

Figure: Figure 1: Loglog plot of the approximation error $\big{|}F(\|x\|)-\frac{1}{P}\sum_{p}f(|\langle\xi_{p},x\rangle|)\big{|}$ for approximating the function $F$ by slicing equation 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 ${\bm{\xi}}^{P}$ and $1000$ realizations of $x$. The kernel parameters are set by the median rule with scaling factor $\gamma=1$. We fit a regression line in the loglog plot for each method to estimate the convergence rate $r$, see also Table 1.
Refer to caption: x1.png

We examine the approximation error in equation 4 numerically. To this end, we draw a sample $x$ from $\mathcal{N}(0,0.1I)$ and evaluate the absolute error $\big{|}F(\|x\|)-\frac{1}{P}\sum_{p=1}^{P}f(|\langle\xi_{p},x\rangle|)\big{|}$. We average this error over $50$ realizations of ${\bm{\xi}}^{P}$ (whenever ${\bm{\xi}}^{P}$ is random) and $1000$ samples of $x$. The average results for the scale factor $\gamma=1$ and dimension $d\in\{3,10,50\}$ with the Gauss, Laplace and Matérn kernel are given in Figure 1.
We observe that all methods despite the spherical designs for the Gauss kernel converge with rate $\mathcal{O}(P^{-r})$ for some $r>0$.
To estimate the rate $r$ numerically, we fit a regression line in the loglog plot.
The resulting rates $r$ are given in Table 1. Further plots and tables are given in Appendix J.1 considering the negative distance kernel, higher dimensions and smaller/larger kernel widths.

Overall, the distance QMC points perform best in most examples. Only when the (provably optimal) spherical design is applicable, it outperforms the distance for smooth kernels and reaches machine precision already for $P\approx 250$.
In accordance with Definition 2, the benefits of QMC slicing are better for smooth kernels and in low dimensions. But also for $d=50$, a significant advantage of QMC slicing is visible.
In particular, we often observe much faster convergence rates than
the proven the worst case error rates $r=\frac{d}{2(d-1)}$ on $H^{s}(\mathbb{S}^{d-1})$ with the distance QMC sequence, see Section 3.
Furthermore, the slight advantage of the distance QMC points versus the spherical designs for the Laplace kernel is because the former ones are chosen specifically for symmetric functions.

### 4.3 Fast Kernel Summation

Figure: Figure 2: Loglog plot of the relative $L^{1}$ approximation error versus computation time for computing the kernel summations equation 1 with different kernels and methods. We use the Letters dataset ($M=N=20000$ points), MNIST (reduced to dimension $d=20$ via PCA, $M=N=60000$ points) and FashionMNIST (reduced to dimension $d=30$ via PCA, $M=N=60000$ points). We run each method $10$ times. The shaded area indicates the standard deviation of the error. For Fourier slicing, we use $P=5\cdot 2^{k}$ slices for $k=1,...,10$. In order to obtain similar computation times, we use $5\cdot 2^{k-1}$ slices for RFF-$10$ slicing and $D=2P$ features for RFF and ORF.
Refer to caption: x13.png

Finally, we test our kernel approximation for computing the whole kernel sums from equation 1. For computing the one-dimensional kernel sums $s_{m}=\sum_{n=1}^{N}w_{n}f(|\langle\xi_{p},x_{n}-y_{m}\rangle|)$, we use the following methods combined with either random or QMC points on the sphere:

- -
(QMC) Sorting-Slicing: For the negative distance kernel, we use the sorting algorithm from Hertrich et al. (), see also (Hertrich, , 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 summation can be found in the text book (Plonka et al., , Sec 7.5). Similar as in (Hertrich, , Sec 2.3), we do not evaluate the one-dimensional basis functions, but directly compute the Fourier transforms. We revisit the background on the one-dimensional fast Fourier summation and specify the used parameters in Appendix I.
- -
(QMC) RFF-$k$ Slicing: For positive definite kernels (Gauss, Laplace, Matérn), we use one-dimensional random Fourier features with $k$ features for the basis function $f$.
For iid or orthogonal slices and $k=1$, this approach is related to RFF and ORF, as outlined in Appendix G.

As QMC designs ${\bm{\xi}}^{P}$, we use the distance QMC points 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 points to obtain an unbiased estimator.
We evaluate the kernel sums on Letters dataset ($d=16$, Slate, ), MNIST (reduced to $d=20$ dimensions via PCA, LeCun et al., ) and FashionMNIST (reduced to $d=30$ dimensions via PCA, Xiao et al., ), where $(x_{1},...,x_{N})$ and $(y_{1},...,y_{M})$ constitute the whole dataset and the weights $(w_{1},...,w_{N})$ are set to $1$. In particular, we have $M=N=20000$ for the Letters dataset and $M=N=60000$ for MNIST and FashionMNIST. Then, we approximate the vector $s=(s_{1},...,s_{M})$ from equation 1 and report the absolute error $\|s-s_{\text{true}}\|_{1}$. We choose the kernel parameters by the median rule with scale factor $\gamma=1$ based on $1000$ example pairs $(x,y)$. We benchmark the computation times on a single thread of an AMD Ryzen Threadripper 7960X CPU and compare our results with RFF, ORF and (non-QMC) slicing. Since the QMC sequences ${\bm{\xi}}^{P}$ depend neither on the dataset nor on the kernel, we consider its construction not as a part of the computation time. For the Gauss kernel, we also compare with QMC (Sobol) RFF (Avron et al., ), which is not applicable for the other kernels. We use $P=5\cdot 2^{k}$ slices for $k=1,...,10$ in the slicing method. In order to obtain similar computation times, we use $5\cdot 2^{k-1}$ slices for RFF-$10$ slicing and $D=2P$ features for RFF and ORF.

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 higher dimensional datasets (including the full MNIST and FashionMNIST with $d=784$) and a GPU comparison are included in Appendix J.2.
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.

## 5 Conclusions

##### Summary and Outlook

We proposed a slicing approach to compute large kernel sums in $\mathcal{O}(P(N+M+N_{\mathrm{ft}}\log N_{\mathrm{ft}}))$ instead of the naïve $\mathcal{O}(NM)$ arithmetic operations. In the case of iid directions, we proved error bounds with rate $\mathcal{O}(1/\sqrt{P})$. 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\leq 100$.
In the future, we want to improve our theoretical analysis on QMC slicing 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 $\xi\to g_{x}(\xi)$ from Section 3.2 are always symmetric. From a practical side, we want to apply the slicing approach in some actual applications.

##### Limitations

In the numerical part, we observe significantly better error rates than we can prove theoretically.
Moreover, the computation of the QMC directions can be very costly and depends strongly on 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 a well-known effect for most QMC methods.

#### Acknowledgments

JH acknowledges funding by the EPSRC programme grant The Mathematics of Deep Learning with reference EP/V026259/1. TJ acknowledges funding by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy EXC-2046/1, Projekt-ID 390685689 (The Berlin Mathematics Research Center MATH+).
MQ acknowledges funding by the German Research Foundation (DFG) with reference STE 571/19-2 and project number 495365311, within the Austrian Science Fund (FWF) SFB 10.55776/F68 Tomography Across the Scales.
For open access purposes, the authors have applied a CC BY public copyright license to any author-accepted manuscript version arising from this submission.
We gratefully thank Gabriele Steidl for fruitful conversations.

## Appendix A Basis Function Pairs ( F , f ) 𝐹 𝑓 (F,f) ( italic_F , italic_f )

We include a list of pairs of basis functions $(F,f)$ that fulfill the slicing formula and the corresponding Fourier transforms in Table 2. The pairs are taken from (Hertrich, , Table 1). We use the convention that the Fourier transform of a function $g\in L_{1}(\mathbb{R}^{d})$ is defined by

$$ $\mathcal{F}_{d}[g](\omega)=\int_{\mathbb{R}^{d}}g(x)\mathrm{e}^{-2\pi\mathrm{i }\langle\omega,x\rangle}\mathrm{d}x.$ (12) $$

The basis functions from the table involve a couple of special functions, which are defined as follows:

- -
Gamma function: $\Gamma(z)=\int_{0}^{\infty}t^{z-1}\mathrm{e}^{-t}\mathrm{d}t$, $\mathrm{Re}(z)>0$,
- -
Modified Bessel function of first kind: $I_{\alpha}(x)=\sum_{m=0}^{\infty}\frac{1}{m!\Gamma(m+\alpha+1)}\left(\frac{x}{
2}\right)^{2m+\alpha},$
- -
Modified Bessel function of second kind:
$K_{\alpha}(x)=\frac{\pi}{2}\frac{I_{-\alpha(x)}-I_{\alpha}(x)}{\sin(\alpha\pi)},$
- -
Modified Struve function: $\mathbf{L}_{\alpha}(x)=\sum_{m=0}^{\infty}\frac{1}{\Gamma(m+\frac{3}{2})\Gamma
(m+\alpha+\frac{3}{2})}\left(\frac{x}{2}\right)^{2m+\alpha+1}$.

The formula for the Fourier transform of the Matérn kernel (and thus for the Laplace kernel with $\nu=1/2$) can be found in (Williams & Rasmussen, , (4.15)). Consequently, we recover the well-known results that the Fourier transforms of the Gauss, Laplace and Matérn kernels are the Fourier transforms of the Gauss, Cauchy and Student-$t$ (with $2\nu$ degrees of freedom) distribution.
To compute the Fourier transforms $\mathcal{F}_{1}^{-1}[f(|\cdot|)]$, we apply (Rux et al., , Prop 3.1), see also (Hertrich, , 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. () and (Wendland, , Sect 8.3) for details.
For the sliced Laplace $f(x)=\exp({-\alpha x})$, we have by (Hertrich, , (3)) and (Gradshteyn & Ryzhik, , 3.387.5)

$$ $F(t)=\int_{0}^{1}\exp({-\alpha st})(1-s^{2})^{\frac{d-3}{2}}\mathrm{d}s=\frac{ \sqrt{\pi}2^{\frac{d-4}{2}}\Gamma(\tfrac{d-1}{2})}{(\alpha t)^{\frac{d-2}{2}}} \left(I_{\frac{d-2}{2}}(-\alpha t)+\mathbf{L}_{\frac{d-2}{2}}(-\alpha t)\right).$ (13) $$

**Table 2: Basis functions $F$ for different kernels $K(x,y)=F(\|x-y\|)$ and corresponding basis functions $f$ from $\mathrm{k}(x,y)=f(|x-y|)$. We added the inverse Fourier transforms $\mathcal{F}_{d}^{-1}[F(\|\cdot\|)]$ and $\mathcal{F}_{1}^{-1}[f(|\cdot|)]$ to the table.**
| Kernel | $F(x)$ | $\mathcal{F}_{d}^{-1}[F(\|\cdot\|)](\|\omega\|)$ | $f(x)$ | $\mathcal{F}_{1}^{-1}[f](|\omega|)$ |
| --- | --- | --- | --- | --- |
| Gauss | $\exp(-\frac{x^{2}}{2\sigma^{2}})$ | $(2\pi\sigma^{2})^{d/2}\exp(-2\pi^{2}\sigma^{2}\omega^{2})$ | ${{}_{1}}F_{1}(\tfrac{d}{2};\tfrac{1}{2};\tfrac{-x^{2}}{2\sigma^{2}})$ | $\frac{\pi\sigma\exp(-2\pi^{2}\sigma^{2}\omega^{2})(2\pi^{2}\sigma^{2}\omega^{2<br>})^{(d-1)/2}}{2\Gamma(\frac{d}{2})}$ |
| Laplace | $\exp(-\alpha x)$ | $\frac{\Gamma(\frac{d+1}{2})2^{d}\pi^{\frac{d-1}{2}}}{\alpha^{d}}(1+\frac{4\pi^<br>{2}\omega^{2}}{\alpha^{2}})^{-\frac{d+1}{2}}$ | $\sum_{n=0}^{\infty}\frac{(-1)^{n}\alpha^{n}\sqrt{\pi}\Gamma(\frac{n+d}{2})}{n!<br>\Gamma(\frac{d}{2})\Gamma(\frac{n+1}{2})}x^{n}$ | $\frac{\Gamma(\frac{d+1}{2})2^{d}\pi^{\frac{d-1}{2}}|\omega|^{d-1}}{\Gamma(<br>\frac{d}{2})\alpha^{d}}(1+\frac{4\pi^{2}\omega^{2}}{\alpha^{2}})^{-\frac{d+1}{<br>2}}$ |
| Sliced Laplace | See equation 13 | $\frac{\Gamma(\frac{d}{2})}{\pi^{d/2}|\omega|^{d-1}}\frac{2\alpha}{\alpha^{2}+4<br>\pi^{2}\omega^{2}}$ | $\exp(-\alpha x)$ | $\frac{2\alpha}{\alpha^{2}+4\pi^{2}\omega^{2}}$ |
| Matérn | $\frac{2^{1-\nu}}{\Gamma(\nu)}(\tfrac{\sqrt{2\nu}}{\beta}x)^{\nu}K_{\nu}(\tfrac<br>{\sqrt{2\nu}}{\beta}x)$ | $\frac{\Gamma(\frac{2\nu+d}{2})2^{\frac{d}{2}}\pi^{\frac{d}{2}}\beta^{d}}{<br>\Gamma(\nu)\nu^{\frac{d}{2}}}(1+\frac{2\pi^{2}\beta^{2}\omega^{2}}{\nu})^{-<br>\frac{2\nu+d}{2}}$ | (Hertrich, , Appendix C) | $\frac{\Gamma(\frac{2\nu+d}{2})2^{\frac{d}{2}}\pi^{d}\beta^{d}|\omega|^{d-1}}{<br>\Gamma(\frac{d}{2})\Gamma(\nu)\nu^{\frac{d}{2}}}(1+\frac{2\pi^{2}\beta^{2}<br>\omega^{2}}{\nu})^{-\frac{2\nu+d}{2}}$ |
| Riesz for $r\in(0,2)$ | $-x^{r}$ | not a function | $-\frac{\sqrt{\pi}\Gamma(\frac{d+r}{2})}{\Gamma(\frac{d}{2})\Gamma(\frac{r+1}{2<br>})}x^{r}$ | not a function |
| Thin Plate Spline | $x^{2}\log(x)$ | not a function | $dx^{2}\log(x)+C_{d}x^{2}$, with<br>$C_{d}=\frac{d}{2}\big{(}H_{\frac{d}{2}}-2+\log(4)\big{)}$ | not a function |

## Appendix B Spherical Sobolev Spaces

We briefly introduce spherical harmonics and Sobolev spaces following Atkinson & Han ().
We denote by
$|\mathbb{S}^{d-1}|=\frac{2\pi^{d/2}}{\Gamma(d/2)}$ the volume of the unit sphere $\mathbb{S}^{d-1}$.
Let $Y_{n}^{k}$, $k=1,\dots,N_{n,d}$ be an $L_{2}$-orthonormal basis of spherical harmonics of degree $n\in\mathbb{N}_{0}$, 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

$$ $N_{n,d}=\frac{(2n+d-2)(n+d-3)!}{n!(d-2)!}\simeq\frac{2}{(d-2)!}n^{d-2}\quad \text{for}\quad n\to\infty.$ (14) $$

The spherical harmonics $Y_{n}^{k}$, $n\in\mathbb{N}_{0}$, $k=1,\dots,N_{n,d}$ form an orthonormal basis of $L_{2}(\mathbb{S}^{d-1})$.

The spherical Sobolev space $H^{s}(\mathbb{S}^{d-1})$ for $s\in\mathbb{R}$ can be defined as the completion of $C^{\infty}(\mathbb{S}^{d-1})$ with respect to the norm

$$ $\|g\|_{H^{s}(\mathbb{S}^{d-1})}^{2}=\sum_{n=0}^{\infty}\sum_{k=1}^{N_{n,d}} \left(n+\tfrac{d-2}{2}\right)^{2s}|\langle g,Y_{n}^{k}\rangle_{L_{2}(\mathbb{S }^{d-1})}|^{2},$ (15) $$

where $\langle g,Y_{n}^{k}\rangle_{L_{2}(\mathbb{S}^{d-1})}=\int_{\mathbb{S}^{d-1}}g(
\xi)Y_{n}^{k}(\xi)\mathrm{d}\xi$.
Note that the factor $(n+\frac{d-2}{2})^{2d}$ can be replaced by another one with the same asymptotic behavior with respect to $n$, yielding an equivalent norm.
For $s=0$, we can identify $H^{s}(\mathbb{S}^{d-1})$ with $L_{2}(\mathbb{S}^{d-1})$.
The Sobolev spaces are nested in the sense that $H^{s}(\mathbb{S}^{d-1})\subset H^{t}(\mathbb{S}^{d-1})$ whenever $s>t$.
If $s>\frac{d-1}{2}$, each function in the Sobolev space $H^{s}(\mathbb{S}^{d-1})$ is continuous (more specifically, it has a continuous representative).
If $s$ is an integer, then $H^{s}(\mathbb{S}^{d-1})$ consists of all functions whose (distributional) derivatives up to order $s$ are square integrable, cf. (Dai & Xu, , Sect 4.5 and 4.8).
An alternative characterization of the Sobolev norm uses the Laplace–Beltrami operator $\Delta_{*}$, which consists of the spherical part of the Laplace, and is given by

$$ $\|g\|_{H^{s}(\mathbb{S}^{d-1})}=\left\|\left(-\Delta_{*}+(\tfrac{d-2}{2})^{2} \right)^{s/2}g\right\|_{L^{2}(\mathbb{S}^{d-1})}.$ $$

If $s$ is an even integer, the operator applied to $g$ is a usual differentiable operator, otherwise it is a pseudodifferential operator.

## Appendix C Proof of Theorem 1

i): Since $F$ is continuous and positive definite, so is $f$ by (Rux et al., , Corollary 4.11). Further, because $\mathbb{E}[f\left(|\langle\xi,x\rangle|\right)]=F(\|x\|)$ and due to the fact that positive definite functions are maximal in the origin, we deduce

$$ $\mathbb{V}_{d}[f]=\mathbb{E}_{\xi\sim\mathcal{U}_{\mathbb{S}^{d-1}}}[f(| \langle\xi,x\rangle|)^{2}]-F(\|x\|)^{2}\leq f(0)^{2}-F(\|x\|)^{2}.$ $$

ii): We write ${\bm{\xi}}=\left(\xi_{1},...,\xi_{d}\right)^{T}$, where $\xi_{d}\in[-1,1]$ denotes the $d$-th component of ${\bm{\xi}}$. We assume w.l.o.g. that $x=\lambda\bm{e}_{d}$ with $\lambda\neq 0$ and recall that we can decompose the unnormalized surface measure on $\mathbb{S}^{d-1}$ as

$$ ${\mathrm{d}}\mathbb{S}^{d-1}({\bm{\xi}})={\mathrm{d}}\mathbb{S}^{d-2}(\bm{\eta })(1-t^{2})^{\frac{d-3}{2}}{\mathrm{d}}t,$ $$

where ${\bm{\xi}}=\sqrt{1-t^{2}}\bm{\eta}+t\bm{e}_{d}$ and $\bm{\eta}\in\mathbb{S}^{d-2}\times\{0\}$, see (Atkinson & Han, , (1.16)). Consequently,

$$ $\displaystyle\mathbb{E}_{{\bm{\xi}}\sim\mathcal{U}_{\mathbb{S}^{d-1}}}\left[| \langle{\bm{\xi}},x\rangle|^{2r}\right]=\frac{\|x\|^{2r}}{|\mathbb{S}^{d-1}|} \int_{\mathbb{S}^{d-1}}|\xi_{d}|^{2r}{\mathrm{d}}\mathbb{S}^{d-1}({\bm{\xi}})$ $\displaystyle~{}=\frac{\|x\|^{2r}}{|\mathbb{S}^{d-1}|}\int_{\mathbb{S}^{d-2}} \int_{-1}^{1}|t|^{2r}(1-t^{2})^{\frac{d-3}{2}}{\mathrm{d}}\mathbb{S}^{d-2}(\bm {\eta}){\mathrm{d}}t=\frac{\|x\|^{2r}}{|\mathbb{S}^{d-1}|}\left|\mathbb{S}^{d- 2}\right|\int_{-1}^{1}(t^{2})^{r}(1-t^{2})^{\frac{d-3}{2}}{\mathrm{d}}t.$ $$

Further, with $B(\cdot,\cdot)$ denoting the Beta function and the substitution $u=t^{2}$, we have

$$ $\int_{-1}^{1}t^{2r}(1-t^{2})^{\frac{d-3}{2}}{\mathrm{d}}t=\int_{0}^{1}u^{r- \frac{1}{2}}(1-u)^{\frac{d-1}{2}-1}{\mathrm{d}}u=B\left(r+\frac{1}{2},\frac{d- 1}{2}\right)=\frac{\Gamma\left(r+\frac{1}{2}\right)\Gamma\left(\frac{d-1}{2} \right)}{\Gamma\left(r+\frac{d}{2}\right)},$ (16) $$

and

$$ $\frac{|\mathbb{S}^{d-2}|}{|\mathbb{S}^{d-1}|}=\frac{2\pi^{\frac{d-1}{2}}}{ \Gamma\left(\frac{d-1}{2}\right)}\frac{\Gamma\left(\frac{d}{2}\right)}{2\pi^{ \frac{d}{2}}}=\frac{\Gamma\left(\frac{d}{2}\right)}{\sqrt{\pi}\Gamma\left( \frac{d-1}{2}\right)}.$ $$

Combining both expressions yields

$$ $\mathbb{E}_{{\bm{\xi}}\sim\mathcal{U}_{\mathbb{S}^{d-1}}}\left[|\langle{\bm{ \xi}},x\rangle|^{2r}\right]=\|x\|^{2r}\frac{\Gamma(r+\frac{1}{2})}{\Gamma(r+ \frac{d}{2})}\frac{\Gamma(\frac{d}{2})}{\sqrt{\pi}},$ (17) $$

and hence

$$ $\frac{\pi\Gamma\left(\frac{d+r}{2}\right)^{2}}{\Gamma\left(\frac{d}{2}\right)^ {2}\Gamma\left(\frac{r+1}{2}\right)^{2}}\mathbb{E}_{{\bm{\xi}}\sim\mathcal{U}_ {\mathbb{S}^{d-1}}}\left[|\langle{\bm{\xi}},x\rangle|^{2r}\right]=\frac{\|x\|^ {2r}\sqrt{\pi}\Gamma(r+\frac{1}{2})}{\Gamma(\frac{r+1}{2})^{2}}\frac{\Gamma( \frac{d+r}{2})^{2}}{\Gamma(\frac{d}{2})\Gamma(r+\frac{d}{2})}.$ $$

The asymptotic expansion $\Gamma(z+a)/\Gamma(z)\sim z^{a}$ for $z\to\infty$, see (, [5.11.12](https://dlmf.nist.gov/5.11)), implies that

$$ $\lim_{d\to\infty}\frac{\Gamma(\frac{d+r}{2})^{2}}{\Gamma(\frac{d}{2})\Gamma(r+ \frac{d}{2})}=1.$ $$

On the other hand, $\frac{\Gamma(\frac{d+r}{2})^{2}}{\Gamma(\frac{d}{2})\Gamma(r+\frac{d}{2})}$ is increasing with respect to $d$ because the identity $\Gamma(z+1)=z\Gamma(z)$ yields that

$$ $\frac{\Gamma(\frac{d+2+r}{2})^{2}}{\Gamma(\frac{d+2}{2})\Gamma(r+\frac{d+2}{2} )}\bigg{/}\frac{\Gamma(\frac{d+r}{2})^{2}}{\Gamma(\frac{d}{2})\Gamma(r+\frac{d }{2})}=\frac{\left(\frac{d+r}{2}\right)^{2}}{(\frac{d}{2})\left(r+\frac{d}{2} \right)}=\frac{\left(d+r\right)^{2}}{2dr+d^{2}}\geq 1.$ $$

Hence, we see that

$$ $\frac{\Gamma(\frac{d+r}{2})^{2}}{\Gamma(\frac{d}{2})\Gamma(r+\frac{d}{2})}\leq 1,$ (18) $$

and thus finally

$$ $\displaystyle\mathbb{V}_{d}[f]\leq\mathbb{E}_{{\bm{\xi}}\sim\mathcal{U}_{ \mathbb{S}^{d-1}}}\left[f(|\langle\xi,x\rangle|)^{2}\right]$ $\displaystyle=\frac{\pi\Gamma\left(\frac{d+r}{2}\right)^{2}}{\Gamma\left(\frac {d}{2}\right)^{2}\Gamma\left(\frac{r+1}{2}\right)^{2}}\mathbb{E}_{{\bm{\xi}} \sim\mathcal{U}_{\mathbb{S}^{d-1}}}\left[|\langle{\bm{\xi}},x\rangle|^{2r}\right]$ $\displaystyle\leq\frac{\sqrt{\pi}\,\Gamma(r+\frac{1}{2})}{\Gamma(\frac{r+1}{2} )^{2}}\|x\|^{2r},$ $$

which proves equation 7.

iii): We move on to the thin plate spline kernel with $f(t)=d|t|^{2}\log(|t|)+C_{d}|t|^{2}.$ As above, w.l.o.g. we assume that $x=s\bm{e}_{d}$ with $s\geq 0$. Then

$$ $\displaystyle\mathbb{E}\left[f(|\langle\xi,x\rangle|)^{2}\right]$ $\displaystyle={}$ $\displaystyle\frac{|\mathbb{S}^{d-2}|}{|\mathbb{S}^{d-1}|}2\int_{0}^{1}\left( ds^{2}\xi_{d}^{2}\log(s\xi_{d})-C_{d}s^{2}\xi_{d}^{2}\right)^{2}(1-\xi_{d}^{2} )^{\frac{d-3}{2}}{\mathrm{d}}\xi_{d}$ $\displaystyle={}$ $\displaystyle\frac{\Gamma\left(\frac{d}{2}\right)}{\sqrt{\pi}\Gamma\left(\frac {d-1}{2}\right)}\frac{d^{2}}{2}\int_{0}^{1}s^{4}\xi_{d}^{4}\left(2\log\left( \xi_{d}\right)+2\log(s)+H_{\frac{d}{2}}-2+\log(4)\right)^{2}(1-\xi_{d}^{2})^{ \frac{d-3}{2}}{\mathrm{d}}\xi_{d}$ $\displaystyle={}$ $\displaystyle\frac{d^{2}\Gamma\left(\frac{d}{2}\right)}{2\sqrt{\pi}\Gamma\left (\frac{d-1}{2}\right)}s^{4}\left(4\int_{0}^{1}\xi_{d}^{4}\log^{2}(\xi_{d})(1- \xi_{d}^{2})^{\frac{d-3}{2}}{\mathrm{d}}\xi_{d}\right.$ (19) $\displaystyle~{}\qquad\qquad\left.+4(2\log(s)+H_{\frac{d}{2}}-2+\log(4))\int_{ 0}^{1}\xi_{d}^{4}\log(\xi_{d})(1-\xi_{d}^{2})^{\frac{d-3}{2}}{\mathrm{d}}\xi_{ d}\right.$ (20) $\displaystyle~{}\qquad\qquad\left.+\left(2\log(s)+H_{\frac{d}{2}}-2+\log(4) \right)^{2}\int_{0}^{1}\xi_{d}^{4}(1-\xi_{d}^{2})^{\frac{d-3}{2}}{\mathrm{d}} \xi_{d}\right).$ (21) $$

We analyze the terms separately.
Denote by $\psi^{(n)}$ the $n$-th derivative of the digamma function and by $\gamma\approx 0.57$ the Euler–Mascheroni constant.
In the following, we use the asymptotics

$$ $H_{x}=\log(x)+\gamma+\mathcal{O}(x^{-1}),\quad\psi^{(0)}(x)=\log(x)+\mathcal{O }(x^{-1}),\quad\psi^{(1)}(x)=\mathcal{O}\left(x^{-1}\right).$ $$

By (Gradshteyn & Ryzhik, , 4.261.21), we have

$$ $\displaystyle\int_{0}^{1}\xi_{d}^{4}\log^{2}(\xi_{d})(1-\xi_{d}^{2})^{\frac{d- 3}{2}}{\mathrm{d}}\xi_{d}$ $\displaystyle\qquad=\frac{3\sqrt{\pi}\,\Gamma\left(\frac{d-1}{2}\right)}{32 \Gamma\left(\frac{d+4}{2}\right)}\left(\psi^{(1)}(\tfrac{5}{2})-\psi^{(1)}( \tfrac{d+4}{2})+\left(\psi^{(0)}(\tfrac{5}{2})-\psi^{(0)}(\tfrac{d+4}{2}) \right)^{2}\right),$ $\displaystyle\qquad=\frac{3\sqrt{\pi}\,\Gamma\left(\frac{d-1}{2}\right)}{32 \Gamma\left(\frac{d+4}{2}\right)}\left(\log(d)^{2}-2\left(\psi^{(0)}\left( \tfrac{5}{2}\right)+\log(2)+\mathcal{O}(\tfrac{1}{d})\right)\log(d)\right.$ $\displaystyle\qquad\qquad\left.+\psi^{(1)}\left(\tfrac{5}{2}\right)+\left(\psi ^{(0)}\left(\tfrac{5}{2}\right)+\log(2)\right)^{2}+\mathcal{O}\left(d^{-1} \right)\right),$ $$

and

$$ $\displaystyle\int_{0}^{1}\xi_{d}^{4}\log(\xi_{d})(1-\xi_{d}^{2})^{\frac{d-3}{2 }}{\mathrm{d}}\xi_{d}$ $\displaystyle=-\frac{\sqrt{\pi}\Gamma\left(\frac{d-1}{2}\right)}{16\Gamma\left (\frac{d+4}{2}\right)}\left(-8+3H_{1+\frac{d}{2}}+\log(64)\right)$ $\displaystyle=-\frac{\sqrt{\pi}\Gamma\left(\frac{d-1}{2}\right)}{16\Gamma\left (\frac{d+4}{2}\right)}\big{(}3\log(d)+3\gamma-8+\log(8)+\mathcal{O}(d^{-1}) \big{)}.$ $$

Recall from equation 16 that

$$ $\displaystyle\int_{0}^{1}\xi_{d}^{4}(1-\xi_{d}^{2})^{\frac{d-3}{2}}{\mathrm{d} }\xi_{d}=\frac{3\sqrt{\pi}\Gamma\left(\frac{d-1}{2}\right)}{8\Gamma\left(\frac {d+4}{2}\right)}.$ $$

Furthermore, using that $\log(4)=2\log(2)$, we obtain

$$ $\displaystyle 2\log(s)+H_{\frac{d}{2}}-2+\log(4)={}$ $\displaystyle 2\log(s)+\log(d)+\gamma+\log(2)-2+\mathcal{O}(d^{-1}),$ $\displaystyle\left(2\log(s)+H_{\frac{d}{2}}-2+\log(4)\right)^{2}={}$ $\displaystyle\big{(}2\log(s)+\log(d)+\gamma+\log(2)-2+\mathcal{O}(d^{-1})\big{ )}^{2}$ $\displaystyle={}$ $\displaystyle 4\log(s)^{2}+\log(d)^{2}+4\log(s)\log(d)$ $\displaystyle+2(\gamma+\log(2)-2)\log(d)+4\log(s)(\gamma+\log(2)-2)$ $\displaystyle+(\gamma+\log(2)-2)^{2}+\mathcal{O}(d^{-1}).$ $$

Plugging this into equation 19 yields

$$ $\displaystyle 4\frac{\Gamma\left(\frac{d+4}{2}\right)}{\sqrt{\pi}\Gamma\left( \frac{d-1}{2}\right)}\frac{\sqrt{\pi}\Gamma\left(\frac{d-1}{2}\right)}{\Gamma \left(\frac{d}{2}\right)d^{2}s^{4}}\mathbb{E}\left[f(|\langle\xi,x\rangle|)^{2 }\right]=\frac{1+2/d}{s^{4}}\mathbb{E}\left[f(|\langle\xi,x\rangle|)^{2}\right]$ $\displaystyle={}$ $\displaystyle\frac{3}{4}\left(\log(d)^{2}-2\left(\psi^{(0)}(\tfrac{5}{2})+\log (2)+\mathcal{O}(\tfrac{1}{d})\right)\log(d)+\psi^{(1)}(\tfrac{5}{2})+\left( \psi^{(0)}(\tfrac{5}{2})+\log(2)\right)^{2}+\mathcal{O}(\tfrac{1}{d})\right)$ $\displaystyle-\frac{1}{8}(2\log(s)+\log(d)+\gamma+\log(2)-2+\mathcal{O}(\tfrac {1}{d}))\left(3\log(d)+3\gamma-8+\log(8)+\mathcal{O}(\tfrac{1}{d})\right)$ $\displaystyle+\frac{3}{4}\left(4\log(s)^{2}+\log(d)^{2}+4\log(s)\log(d)+2( \gamma+\log(2)-2)\log(d)\right.$ $\displaystyle\qquad\left.+4\log(s)(\gamma+\log(2)-2)+(\gamma+\log(2)-2)^{2}+ \mathcal{O}\left(d^{-1}\right)\right)$ $\displaystyle={}$ $\displaystyle\log(d)^{2}\left(\frac{3}{4}-\frac{3}{2}+\frac{3}{4}\right)$ $\displaystyle+\log(d)\left(-\frac{3}{2}\left(\psi^{(0)}(\tfrac{5}{2})+\log(2)+ \mathcal{O}(\tfrac{1}{d})\right)-\frac{1}{2}\left(3\gamma-8+\log(8)\right)\right.$ $\displaystyle\qquad\qquad\left.-\frac{3}{2}\left((2\log(s)+\gamma+\log(2)-2 \right)+\frac{3}{4}\left(4\log(s)+2(\gamma+\log(2)-2)\right)\right)$ $\displaystyle+{3}\log(s)^{2}-\log(s)\left(3\gamma-8+\log(8)\right)$ $\displaystyle+\frac{3}{4}\left(\psi^{(1)}(\tfrac{5}{2})+\left(\psi^{(0)}( \tfrac{5}{2})+\log(2)\right)^{2}\right)+2(\gamma+\log(2)-2)^{2}+\mathcal{O} \left(d^{-1}\right)$ $\displaystyle={}$ $\displaystyle\log(d)\left(-\frac{3}{2}\psi^{(0)}(\tfrac{5}{2})-\frac{3}{2} \gamma-{3}\log(2)++\mathcal{O}(\tfrac{1}{d})\right)+{3}\log(s)^{2}-\log(s)(3 \gamma+\log(8)-8)$ $\displaystyle+\frac{3}{4}\left(\psi^{(1)}(\tfrac{5}{2})+\left(\psi^{(0)}( \tfrac{5}{2})+\log(2)\right)^{2}\right)+2(\gamma+\log(2)-2)^{2}+\mathcal{O}(d^ {-1}).$ $$

With equation 8 and the identity $\psi^{(0)}(\tfrac{5}{2})=-2\log(2)-\gamma+\frac{8}{3}$, we obtain

$$ $\frac{1+2/d}{s^{4}}\mathbb{E}\left[f(|\langle\xi,x\rangle|)^{2}\right]\\ =3\log(s)^{2}+c_{1}\log(s)+c_{2}+\mathcal{O}\left(d^{-1}\log d\right).$ $$

iv):
By (Hertrich, , Thm 3), the transformed function has the form

$$ $f(s)=\sum_{n=0}^{\infty}b_{n}x^{n}\quad\text{with}\quad b_{n}=\frac{\sqrt{\pi} \Gamma\left(\frac{n+d}{2}\right)}{\Gamma\left(\frac{d}{2}\right)\Gamma\left( \frac{n+1}{2}\right)}a_{n}.$ $$

Moreover,

$$ $\displaystyle f(s)^{2}=\sum_{n=0}^{\infty}c_{n}s^{n}:=\sum_{n=0}^{\infty}\left (\sum_{k=0}^{n}b_{k}b_{n-k}\right)s^{n}$ $$

and, by equation 17,

$$ $\mathbb{E}_{\xi\sim\mathcal{U}_{\mathbb{S}^{d-1}}}\left[f(|\langle\xi,x\rangle |)^{2}\right]=\sum_{n=0}^{\infty}\frac{\Gamma\left(\frac{d}{2}\right)\Gamma \left(\frac{n+1}{2}\right)}{\sqrt{\pi}\Gamma\left(\frac{n+d}{2}\right)}c_{n}\| x\|^{n}.$ $$

Using that $d$ is odd and applying the identities $\Gamma(z+1)=\Gamma(z)z$ and $\Gamma(\tfrac{1}{2})=\sqrt{\pi}$, we have

$$ $\displaystyle\tilde{c}_{k,n,d}$ $\displaystyle\coloneqq\frac{\Gamma\left(\frac{d}{2}\right)\Gamma\left(\frac{n+ 1}{2}\right)}{\sqrt{\pi}\Gamma\left(\frac{n+d}{2}\right)}\frac{\sqrt{\pi} \Gamma\left(\frac{k+d}{2}\right)}{\Gamma\left(\frac{d}{2}\right)\Gamma\left( \frac{k+1}{2}\right)}\frac{\sqrt{\pi}\Gamma\left(\frac{n-k+d}{2}\right)}{ \Gamma\left(\frac{d}{2}\right)\Gamma\left(\frac{n-k+1}{2}\right)}$ $\displaystyle~{}=\frac{\Gamma\left(\frac{n+1}{2}\right)}{\Gamma\left(\frac{n+1 +d-1}{2}\right)}\frac{\Gamma\left(\frac{k+1+d-1}{2}\right)}{\Gamma\left(\frac{ k+1}{2}\right)}\frac{\sqrt{\pi}\Gamma\left(\frac{n-k+1+d-1}{2}\right)}{\Gamma \left(\frac{d}{2}\right)\Gamma\left(\frac{n-k+1}{2}\right)}$ $\displaystyle~{}=\prod_{j=1}^{\frac{d-1}{2}}\frac{(k+2j-1)(n-k+2j-1)}{(n+2j-1) (2j-1)}.$ $$

Hence, we obtain

$$ $\displaystyle\mathbb{V}_{d}[f](x)=\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}( \tilde{c}_{k,n,d}-1)a_{k}a_{n-k}\right)\|x\|^{n}.$ $$

We apply the identity $4xy=\left(x+y\right)^{2}-\left(x-y\right)^{2}$ to the numerator (with $x=k+2j-1$, $y=n-k+2j-1$) and denominator (with $x=n+2j-1$, $y=2j-1$), and obtain

$$ $\displaystyle\tilde{c}_{k,n,d}$ $\displaystyle~{}=\prod_{j=1}^{\frac{d-1}{2}}\frac{(k+2j-1)(n-k+2j-1)}{(n+2j-1) (2j-1)}=\prod_{j=1}^{\frac{d-1}{2}}\frac{(n+2(2j-1))^{2}-(n-2k)^{2}}{(n+2(2j-1 ))^{2}-n^{2}}$ $\displaystyle~{}=\prod_{j=1}^{\frac{d-1}{2}}\left(1+\frac{n^{2}-(n-2k)^{2}}{(n +2(2j-1))^{2}-n^{2}}\right)=\prod_{j=1}^{\frac{d-1}{2}}\left(1+\frac{k(n-k)}{( 2j-1)(n+2j-1)}\right).$ $$

For $F$ the Laplace kernel, assuming w.l.o.g. that $\alpha=1$, we have $a_{n}:=(-1)^{n}/n!$. Consequently,

$$ $\displaystyle\mathbb{V}_{3}[f](x)$ $\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{k(n-k)}{n+1}\frac{1}{k!(n -k)!}(-\|x\|)^{n}=\sum_{n=0}^{\infty}\frac{(-\|x\|)^{n}}{(n+1)!}\sum_{k=0}^{n} {{n}\choose{k}}(kn-k^{2})$ $\displaystyle=\sum_{n=0}^{\infty}\frac{(-\|x\|)^{n}}{(n+1)!}\left(n^{2}2^{n-1} -n(n+1)2^{n-2}\right)=\frac{1}{4}\sum_{n=0}^{\infty}\frac{(-2\|x\|)^{n}}{(n+1) !}n(n-1)$ $\displaystyle=\frac{1}{4}\sum_{n=2}^{\infty}\frac{(-2\|x\|)^{n}}{(n+1)!}(n+1-2 )n=\frac{1}{4}\sum_{n=2}^{\infty}\frac{(-2\|x\|)^{n}}{(n-1)!}-\frac{1}{2}\sum_ {n=2}^{\infty}\frac{(-2\|x\|)^{n}}{(n+1)!}n.$ (22) $$

For the first term, we obtain

$$ $\displaystyle\sum_{n=2}^{\infty}\frac{(-2\|x\|)^{n}}{(n-1)!}=-2\|x\|\sum_{n=1} ^{\infty}\frac{(-2\|x\|)^{n}}{n!}=-2\|x\|\left(\mathrm{e}^{-2\|x\|}-1\right),$ (23) $$

and for the second

$$ $\displaystyle\sum_{n=2}^{\infty}\frac{(-2\|x\|)^{n}}{(n+1)!}(n+1-1)=\sum_{n=2} ^{\infty}\frac{(-2\|x\|)^{n}}{n!}-\sum_{n=2}^{\infty}\frac{(-2\|x\|)^{n}}{(n+1 )!}$ $\displaystyle~{}=\mathrm{e}^{-2\|x\|}-1+2\|x\|+\frac{1}{2\|x\|}\left(\mathrm{e }^{-2\|x\|}-1+2\|x\|-2\|x\|^{2}\right).$ (24) $$

Plugging equation 23 and equation 24 into equation 22 finally yields

$$ $\displaystyle\mathbb{V}_{3}[f](x)$ $\displaystyle=\frac{-\|x\|}{2}(\mathrm{e}^{-2\|x\|}-1)-\frac{1}{2}\left( \mathrm{e}^{-2\|x\|}-1+2\|x\|+\frac{\mathrm{e}^{-2\|x\|}}{2\|x\|}-\frac{1}{2\| x\|}+1-\|x\|\right)$ $\displaystyle=\frac{\mathrm{e}^{-2\|x\|}}{4\|x\|}\left(\mathrm{e}^{2\|x\|}-(2 \|x\|^{2}+2\|x\|+1)\right)$ $\displaystyle=\frac{1}{4\|x\|}\left(\mathrm{e}^{2\|x\|}-(2\|x\|^{2}+2\|x\|+1) \right)F(\|x\|)^{2}.$ $$

For the Gauss kernel, we can follow exactly the same lines, after considering w.l.o.g. $\sigma^{2}=1/2$ and replacing $\|x\|$ with $\|x\|^{2}$.

## Appendix D The variance of Gauss and Laplace kernel for d > 3 𝑑 3 d>3 italic_d > 3

For $d>3$, $d$ odd, we evaluated $\mathbb{V}_{d}(f)$ for $d=5,...,15$, and make the following conjecture.
Let $T_{n}(f)=\sum_{k=0}^{n}\frac{f^{(k)}(0)}{k!}x^{k}$ denote the Taylor expansion of $f$ in zero up to order $n$.
We conjecture that, for $d\geq 3$, $d$ odd, and the Laplace kernel $F(x)=\mathrm{e}^{-\alpha x}$, it holds that

$$ $\displaystyle\mathbb{V}_{d}[f](x)=\frac{\mathrm{e}^{-2\alpha\|x\|}}{c_{d}( \alpha\|x\|)^{d-2}}\left(\sum_{i=0}^{\frac{d-3}{2}}(\alpha\|x\|)^{2i}(-1)^{ \frac{d-5}{2}-i}c_{i,d}T_{d-2i}\left(\mathrm{e}^{2\alpha\|x\|}\right)+\sum_{i= 0}^{\frac{d-5}{2}}b_{i,d}(-1)^{i}(\alpha\|x\|)^{d+i}\right),$ $$

where

$$ $\displaystyle c_{d}$ $\displaystyle=\left(12+4(d-5)\right)c_{d-2}\quad\mbox{for}\quad d\geq 5\quad \mbox{with}\quad c_{3}=4,$ $\displaystyle c_{0,d}$ $\displaystyle=(d-2)^{2}(d-4)c_{0,d-2}\quad\mbox{for}\quad d\geq 5\quad\mbox{ with}\quad c_{0,3}=1,$ $\displaystyle b_{\frac{d-5}{2},d}$ $\displaystyle=2b_{\frac{d-5}{2},d-2}\quad\mbox{for}\quad d\geq 7\quad\mbox{ with}\quad b_{0,5}=4.$ $$

For the remaining coefficients, we have not found a general simple rule. Our numerical computations indicate that they are positive and that
$\mathbb{V}_{d}[f](x)\|x\|$ is monotonically increasing in $\|x\|$, with upper bound $s_{d}/\alpha$ and quadratic convergence to 0 0 for $x\to 0$. The upper bound $s_{d}$ increases slowly in $d$, and our simulations hint that it might be bounded by a constant independently of $d$, see Figure 3. For the Gauss kernel, the variance has the same form, except that $\alpha^{-1}$ is replaced with $\sigma/\sqrt{2}$ and $\|x\|$ with $\|x\|^{2}$ at all occurrences.

Figure: Figure 3: The scaled variance of the Laplace kernel $f_{L}$ with $\alpha=1$ (in red) and the Gauss kernel $f_{G}$ with $\sigma=1/\sqrt{2}$ (in blue) for dimension $d=3$ (solid lines), $d=9$ (dashed lines) and $d=15$ (dotted lines). The variance is multiplied times $\|x\|$ (Laplace kernel) and $\|x\|^{2}$ (Gauss kernel). We observe that the scaled variance increases monotonically in all cases, seemingly bounded from above by a constant (black solid line).
Refer to caption: x25.png

## Appendix E Randomization of a QMC Sequence

Having a QMC sequence $(\bm{\xi}^{P})_{P}$, see Section 3.1, we can easily obtain an unbiased estimator in equation 4 by considering the QMC sequence $(A\bm{\xi}^{P})_{P}$ with a uniformly chosen orthogonal matrix $A\sim\mathcal{U}_{\mathrm{O}(d)}$
because, according to (Ragozin, , (2.3)), we have

$$ $\mathbb{E}_{A\sim\mathcal{U}_{\mathrm{O}(d)}}[\psi(A\xi^{P}_{p})]=\int_{ \mathrm{O}(d)}\psi(A\xi^{P}_{p})\mathrm{d}\mathcal{U}_{\mathrm{O}(d)}(A)=\int_ {\mathbb{S}^{d-1}}\psi(\eta)\mathrm{d}\mathcal{U}_{\mathbb{S}^{d-1}}(\eta),$ $$

where $\mathcal{U}_{\mathrm{O}(d)}$ is the uniform distribution on the set $\mathrm{O}(d)$ of orthogonal $d\times d$ matrices.
Furthermore, note that $\frac{1}{P}\sum_{i=1}^{P}f(\xi^{P}_{i})$ converges to $\int_{\mathbb{S}^{d-1}}f(\xi)\mathrm{d}\mathcal{U}_{\mathbb{S}^{d-1}}(\xi)$ for $P\to\infty$ and all continuous functions $f$ because $H^{s}(\mathbb{S}^{d-1})$ is dense in $C(\mathbb{S}^{d-1})$ and the weights are all one, cf. (Atkinson, , Sect 5.2).

## Appendix F Proof of Theorem 3

Gauss kernel: A sufficient criterion for the function $g_{x}$ being in the Sobolev space $H^{s}(\mathbb{S}^{d-1})$ is that $f(|\cdot|)$ is $s$ times continuously differentiable.
For $F(t)=\exp(-\frac{t^{2}}{2\sigma^{2}})$, the one-dimensional basis function is given by the confluent hypergeometric function $f(t)={{}_{1}}F_{1}(\tfrac{d}{2};\tfrac{1}{2};\tfrac{-t^{2}}{2\sigma^{2}})$. We know by (Hertrich, , Lem 6) that $f=\mathcal{F}_{1}(g)$ with $g(\omega)=\frac{d\pi\sigma\exp(-2\pi^{2}\sigma^{2}\omega^{2})(2\pi^{2}\sigma^{
2}\omega^{2})^{(d-1)/2}}{\sqrt{2}\Gamma(\frac{d+2}{2})}$, where $\mathcal{F}_{1}$ denotes the one-dimensional Fourier transform.
Then, $g_{k}(\omega)\coloneqq\omega^{k}g(\omega)$ is absolutely integrable for any $k\in\mathbb{N}$ such that the differentiation/multiplication formula for the Fourier transform yields that the $k$-th derivative of $f$ exists and is given by the continuous function $f^{(k)}=(2\pi\mathrm{i})^{k}\mathcal{F}_{1}(g_{k})$.

Riesz kernel:
With fixed $x\in\mathbb{R}^{d}\setminus\{0\}$, we have $g_{x}(\xi)=\frac{\sqrt{\pi}\Gamma(\frac{d+r}{2})}{\Gamma(\frac{d}{2})\Gamma(
\frac{r+1}{2})}|\langle\xi,x\rangle|^{r}$.
We set $\eta\coloneqq\frac{x}{\|x\|}\in\mathbb{S}^{d-1}$, then $x=\|x\|\eta$ and accordingly $g_{x}(\xi)=\|x\|^{r}g_{\eta}(\xi)$.
The Sobolev norm equation 15 reads

$$ $\|g_{x}\|_{H^{s}(\mathbb{S}^{d-1})}^{2}=\sum_{n=0}^{\infty}\sum_{k=1}^{N_{n,d} }\left(n+\tfrac{d-2}{2}\right)^{2s}\left|\|x\|^{r}\int_{\mathbb{S}^{d-1}}g_{ \eta}(\xi)\overline{Y_{n}^{k}(\xi)}\mathrm{d}\xi\right|^{2}.$ $$

The last integral can be computed with the help of the so-called $\lambda$-cosine transform.
By (Rubin, , (3.5) and (3.8)), the $\lambda$-cosine transform of a function $h\in L^{1}(\mathbb{S}^{d-1})$ for $\lambda>-1$ is defined by

$$ $\mathcal{C}^{\lambda}[h](\omega)=\frac{1}{|\mathbb{S}^{d-1}|}\frac{\sqrt{\pi} \Gamma(-\frac{\lambda}{2})}{\Gamma(\frac{d}{2})\Gamma(\frac{\lambda+1}{2})} \int_{\mathbb{S}^{d-1}}h(\theta)|\langle\omega,\theta\rangle|^{\lambda}\mathrm {d}\theta,$ $$

for $\omega\in\mathbb{S}^{d-1}$ and satisfies

$$ $\mathcal{C}^{\lambda}[Y_{n}^{k}](\omega)=Y_{n}^{k}(\omega)\begin{cases}(-1)^{ \frac{n}{2}}\frac{\Gamma(\frac{n-\lambda}{2})}{\Gamma(\frac{n+d+\lambda}{2})}, &n\text{ even},\\ 0,&n\text{ odd}.\end{cases}$ $$

Hence, we obtain that

$$ $\displaystyle\int_{\mathbb{S}^{d-1}}g_{\eta}(\xi)\overline{Y_{n}^{k}(\xi)} \mathrm{d}\xi=\frac{|\mathbb{S}^{d-1}|\Gamma(\frac{d+r}{2})}{\Gamma(-\frac{r}{ 2})}\mathcal{C}^{r}\left[\overline{Y_{n}^{k}}\right](\eta)=(-1)^{\frac{n}{2}} \frac{|\mathbb{S}^{d-1}|\Gamma(\frac{n-r}{2})\Gamma(\frac{d+r}{2})}{\Gamma( \frac{n+d+r}{2})\Gamma(-\frac{r}{2})}\,\overline{Y_{n}^{k}(\eta)}$ $$

if $n$ is even and 0 0 if $n$ is odd.
The addition formula for spherical harmonics,

$$ $\sum_{k=1}^{N_{n,d}}|Y_{n}^{k}(\eta)|^{2}=\frac{N_{n,d}}{|\mathbb{S}^{d-1}|},$ $$

see Atkinson & Han (), and the substitution $n=2m$ yield

$$ $\|g_{x}\|_{H^{s}(\mathbb{S}^{d-1})}^{2}=\sum_{m=0}^{\infty}\left(2m+\frac{d-2} {2}\right)^{2s}\|x\|^{2r}\left(\frac{\Gamma(m-\frac{r}{2})\Gamma(\frac{d+r}{2} )}{\Gamma(m+\frac{d+r}{2})\Gamma(-\frac{r}{2})}\right)^{2}|\mathbb{S}^{d-1}|{N _{2m,d}}.$ $$

The asymptotic relations $N_{n,d}\sim n^{d-2}$ from equation 14 and $\Gamma(z+a)/\Gamma(z)\sim z^{a}$ for $z\to\infty$, see (, [5.11.12](https://dlmf.nist.gov/5.11)), yield that the $m$-th summand behaves asymptotically like a multiple of $(2m)^{2s-2r-2}$,
therefore the series converges if and only if $s<r+\frac{1}{2}$.

Matérn kernel:
The one-dimensional basis function of the Matérn kernel is (Hertrich, , Appx C.1)

$$ $f(t)=\underbrace{{}_{1}F_{2}\left(\tfrac{d}{2};\tfrac{1}{2},1-\nu;\tfrac{\nu t ^{2}}{2\beta^{2}}\right)}_{\eqqcolon f_{1}(t)}-\underbrace{|t|^{2\nu}}_{ \eqqcolon f_{2}(t)}\underbrace{\frac{\Gamma(1-\nu)\Gamma(\nu+\tfrac{d}{2}(2\nu )^{\nu}}{\Gamma(\tfrac{d}{2})\Gamma(2\nu+1)\beta^{2\nu}}{}_{1}F_{2}\left(\nu+ \tfrac{d}{2};\nu+\tfrac{1}{2},\nu+1;\tfrac{\nu t^{2}}{2\beta^{2}}\right)}_{ \eqqcolon f_{3}(t)}.$ (25) $$

The above hypergeometric functions $f_{1}$ and $f_{3}$, whose parameters do not contain any nonpositive integers, are entire functions in $\mathbb{R}$.
Hence, the corresponding spherical functions $(g_{1})_{x}$ and $(g_{3})_{x}$ from equation 11 are in $C^{\infty}(\mathbb{S}^{d-1})$, because they can be extended to smooth functions defined on a neighborhood of $\mathbb{S}^{d-1}$.
Furthermore, since $f_{2}$ is a multiple of the one-dimensional basis function of the Riesz kernel, we know from above that $(g_{2})_{x}\in H^{s}(\mathbb{S}^{d-1})$ for $s<2\nu+\tfrac{1}{2}$.
Hence, by (Quellmalz, , Thm 5.2), we see that the product $(g_{2})_{x}(g_{3})_{x}$ is also in $H^{s}(\mathbb{S}^{d-1})$. Note that the referenced theorem is only formulated for in integer $s$, but this can be easily extended using an interpolation argument (Quellmalz, , Sect 5.3), because the multiplication with a smooth function constitutes a continuous operator both $H^{\lfloor s\rfloor}(\mathbb{S}^{d-1})\to H^{\lfloor s\rfloor}(\mathbb{S}^{d-1})$ and $H^{\lfloor s\rfloor+1}(\mathbb{S}^{d-1})\to H^{\lfloor s\rfloor+1}(\mathbb{S}^
{d-1})$.

## Appendix G Relation between Slicing and RFF

We denote by $\mathcal{M}_{+}(\mathbb{R}^{d})$ the space of finite positive Borel measures on $\mathbb{R}^{d}$.
Such measures can be identified with linear functional on the space $C_{0}(\mathbb{R}^{d})$ of continuous functions that vanish at infinity.
The Fourier transform of measures is a linear operator defined by

$$ $\mathcal{F}_{d}\colon\mathcal{M}_{+}(\mathbb{R}^{d})\to C_{0}(\mathbb{R}^{d}), \quad\mathcal{F}_{d}[\mu](x)=\int_{\mathbb{R}^{d}}\mathrm{e}^{-2\pi\mathrm{i} \langle x,v\rangle}\mathrm{d}\mu(v),$ $$

cf. (Plonka et al., , Sect 4.4).
By Bochner’s theorem, the Fourier transform is bijective from $\mathcal{M}_{+}(\mathbb{R}^{d})$ to the set of positive definite functions on $\mathbb{R}^{d}$, see ii) in Section 2 for the definition.
If $\mu$ is a probability measure, i.e., $\mu(\mathbb{R}^{d})=1$, we have $\mathcal{F}_{d}[\mu](0)=1$.

In the following, let $F\circ\|\cdot\|$ be a positive definite function on $\mathbb{R}^{d}$ with $F(0)=1$.

RFF:
Random Fourier features (RFF), see Rahimi & Recht (), use that, by Bochner’s theorem, $F(\|\cdot\|)$
is the Fourier transform of a probability measure $\mu\in\mathcal{M}_{+}(\mathbb{R}^{d})$, i.e.,

$$ $F(\|x\|)=\mathcal{F}_{d}[\mu](x)=\mathbb{E}_{v\sim\mu}\left[\mathrm{e}^{-2\pi \mathrm{i}\langle x,v\rangle}\right].$ (26) $$

Using that $F$ is real-valued and sampling $D\in\mathbb{N}$ instances $v_{p}$ iid from $\mu$, we obtain the approximation

$$ $F(\|x\|)=\mathbb{E}_{v\sim\mu}\left[\mathrm{\cos}(2\pi\langle x,v\rangle) \right]\approx\frac{1}{D}\sum_{p=1}^{D}\cos(2\pi\langle x,v_{p}\rangle).$ (27) $$

By (Altekrüger et al., , Prop C.2), the radial measure $\mu$ can be decomposed as follows. We define $\iota\colon\mathbb{R}^{d}\to[0,\infty),$ $\iota(x)=\|x\|$ and its pushforward measure $\tilde{\mu}\coloneqq\iota_{*}\mu=\mu\circ\iota^{-1}\in\mathcal{M}_{+}([0,
\infty))$,
then we have

$$ $\mu=T(\tilde{\mu}\otimes\mathcal{U}_{\mathbb{S}^{d-1}}),\quad\text{where}\quad T (r,\xi)=r\xi.$ $$

Hence, sampling from $\mu$ can be realized by $v_{p}=r_{p}\xi_{p}$, where $\xi_{p}\sim\mathcal{U}_{\mathbb{S}^{d-1}}$ and $r_{p}\sim\tilde{\mu}$. Then equation 27 becomes

$$ $F(\|x\|)\approx\frac{1}{D}\sum_{p=1}^{D}\cos(2\pi r_{p}\langle x,\xi_{p} \rangle).$ (28) $$

Slicing:
The slicing approach uses the approximation equation 4, i.e.,

$$ $F(\|x\|)\approx\frac{1}{P}\sum_{p=1}^{P}f(|\langle x,\xi_{p}\rangle|),$ (29) $$

where $\xi_{p}\sim\mathcal{U}_{\mathbb{S}^{d-1}}$.
By (Rux et al., , Cor 4.11), the function $f(|\cdot|)$ is positive definite and hence possesses a Fourier transform, which is a probability measure on $\mathbb{R}$ because $f(0)=F(0)=1$.
Applying RFF with $Q$ points to the one-dimensional function $f(|\cdot|)$, we obtain

$$ $F(\|x\|)\approx\frac{1}{PQ}\sum_{p=1}^{P}\sum_{q=1}^{Q}\cos(2\pi r_{p,q} \langle x,\xi_{p}\rangle),$ (30) $$

where $r_{p,q}\sim\mathcal{F}_{1}^{-1}[f(|\cdot|)]$.
According
to (Rux et al., , Cor 4.11),
we have

$$ $\mathcal{F}_{1}^{-1}[f(|\cdot|)]=\mathcal{A}_{1}^{*}\iota_{*}\mu=\mathcal{A}_{ 1}^{*}\tilde{\mu},$ $$

where $\mathcal{A}_{1}^{*}\colon\mathcal{M}_{+}([0,\infty))\to\mathcal{M}_{+}(\mathbb
{R})$ is the symmetrization operator that extends a measure on $[0,\infty)$ to an even measure on the real line and is defined for any $\nu\in\mathcal{M}_{+}([0,\infty))$ and $g\in C_{0}(\mathbb{R})$ by $\langle\mathcal{A}_{1}^{*}\nu,g\rangle=\langle\nu,g+g(-\cdot)\rangle$.
Since the right-hand side of equation 30 is independent of the sign of $r_{p,q}$, it stays the same when we sample $r_{p,q}$ from $\tilde{\mu}$ instead of $\mathcal{A}_{1}^{*}\tilde{\mu}$.
Therefore, if $Q=1$, we see that the right-hand side of equation 30 is the same as the right-hand side of equation 28.
In particular, RFF can be viewed as a special case of slicing.

## Appendix H Relation between Distance QMC Points and Orthogonal Points

Consider the QMC design $\bm{\xi}^{P}$ that is a minimizer equation 10 for $s=\frac{d}{2}$.
To improve the error, Womersley () suggested symmetric QMC designs, meaning that for every point $\xi$ it contains also the antipodal point $-\xi$.
However, since the function equation 11 we want to integrate is symmetric, i.e., $g_{x}(\xi)=g_{x}(-\xi)$, we can discard one of the antipodal points $\xi$ and $-\xi$ and get the same result.
Therefore, we minimize the functional

$$ $\mathcal{E}_{\mathrm{sym}}(\bm{\xi}^{P})=\mathcal{E}((\bm{\xi}^{P},-\bm{\xi}^{ P}))=-2\sum_{p,q=1}^{P}\left(\|\xi^{P}_{p}-\xi^{P}_{q}\|+\|\xi^{P}_{p}+\xi^{P} _{q}\|\right).$ (31) $$

Our numerical trials indicate that indeed the minimizers of $\mathcal{E}_{\mathrm{sym}}$ yield a smaller integration error than the minimizers of $\mathcal{E}$.
Using $\|\xi_{p}\|=1$, we see that

$$ $\displaystyle-\left(\|\xi^{P}_{p}-\xi^{P}_{q}\|+\|\xi^{P}_{p}+\xi^{P}_{q}\| \right)^{2}$ $\displaystyle=-\left(\sqrt{2-2\langle\xi^{P}_{p},\xi^{P}_{q}\rangle}+\sqrt{2+2 \langle\xi^{P}_{p},\xi^{P}_{q}\rangle}\right)^{2}$ $\displaystyle=-4-2\sqrt{4-4\langle\xi^{P}_{p},\xi^{P}_{q}\rangle^{2}},$ $$

attains its minimum if and only if $\langle\xi^{P}_{p},\xi^{P}_{q}\rangle=0$.
Hence, if $P\leq d$ and $\bm{\xi}^{P}$ is an orthonormal system in $\mathbb{R}^{d}$, then $\bm{\xi}^{P}$ a minimizer of $\mathcal{E}_{\mathrm{sym}}$.
However, this argumentation does not work if $P>d$, as we can only choose $d$ orthogonal vectors.

## Appendix I Background on One-Dimensional Fast Fourier Summation

In this section, we review literature about one-dimensional fast Fourier summation used in Section 4.3 and specify the parameters used in our numerical examples.
Fast Fourier summations were proposed in Kunis et al. (); Potts et al. () based on the non-equispaced fast Fourier transform (Beylkin, ; Dutt & Rokhlin, ), which is implemented in several libraries (Knopp et al., ; Keiner et al., ; Shih et al., ). In our numerical examples, we use the Julia library Knopp et al. (). Here, we follow a similar workflow as in Hertrich ().

Let $x_{1},...,x_{N}\in\mathbb{R}^{d}$, $y_{1},...,y_{M}\in\mathbb{R}^{d}$, $w_{1},...,w_{N}\in\mathbb{R}$ and $\mathrm{k}(x,y)=f(|x-y|)=g(x-y)$. We want to compute for $\xi\in\mathbb{S}^{d-1}$, $x_{n,\xi}\coloneqq\langle x_{n},\xi\rangle$ and $y_{m,\xi}\coloneqq\langle y_{m},\xi\rangle$ the one-dimensional kernel sums

$$ $t_{m}=\sum_{n=1}^{N}w_{n}\mathrm{k}(x_{n,\xi},y_{m,\xi})=\sum_{n=1}^{N}w_{n}g( x_{n,\xi}-y_{m,\xi}).$ $$

##### Step 1: Rescaling

For Step 2 and 3, we will need two properties. First, since we will use discrete (fast) Fourier transforms, we require that $x_{n,\xi},y_{m,\xi},x_{n,\xi}-y_{m,\xi}\in[-\frac{1}{2},\frac{1}{2})$.
For the important example of positive definite kernels, which decay to zero, we often can derive explicit formulas for the Fourier transform of $g$ via Bochner’s integral and Rux et al. (), see Table 2 for some examples. In order to use this explicit formula for the fast Fourier summation, we will additionally require that $g(x)\approx 0$ for $|x|>\frac{1}{2}$.

Both properties can be achieved by rescaling the problem. More precisely, let $T<0.5$ and $c=\max\left\{\max_{n=1...,N}\|x_{n}\|,\max_{m=1,...,M}\|y_{m}\|\right\}$. For the case of decaying positive definite kernels, assume that $g(x)\approx 0$ for $|x|>g_{\max}$. Then, it holds that

$$ $t_{m}=\sum_{n=1}^{N}w_{n}g(x_{n,\xi}-y_{m,\xi})=\sum_{n=1}^{N}w_{n}\tilde{g}( \tilde{x}_{n,\xi}-\tilde{y}_{m,\xi}),$ $$

with

$$ $\tilde{g}(x):=g\left(\frac{x}{\tau}\right),\quad\tilde{x}_{n,\xi}:=\tau x_{n, \xi}=\langle\tau x_{n},\xi\rangle,\quad\tilde{y}_{n,\xi}:=\tau y_{n,\xi}= \langle\tau y_{n},\xi\rangle,$ $$

where the constant $\tau:=\min\left\{\frac{T}{c},\frac{1}{2g_{\max}}\right\}$ does not depend on $\xi$.
Then, by definition, the two properties from above are fulfilled. For the rest of the section, we will denote the rescaled points $\tilde{x}_{n,\xi}$, $\tilde{y}_{m,\xi}$ and the rescaled kernel function $\tilde{g}$ again by $x_{n,\xi}$, $y_{m,\xi}$ and $g$.

In the case of the Gauss, Laplace or Matérn kernels, the rescaled kernel $\tilde{\mathrm{k}}(x,y)=\tilde{g}(x-y)$ is again a Gauss, Laplace or Matérn kernel with the altered parameter $\tilde{\sigma}=\tau\sigma$, $\tilde{\alpha}=\alpha/\tau$ or $\tilde{\beta}=\tau\beta$. In our numerics, we set $g_{\max}=5m$ with $m=\sigma=\beta=\frac{1}{\alpha}$ for the Gauss, Laplace and Matérn kernel. 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.

##### Step 2: Computation of the Fourier Coefficients of the Kernel

In the next step, we expand $g$ into its Fourier series on $[-\frac{1}{2},\frac{1}{2})$ and truncate it by

$$ $g(x)=\sum_{k\in\mathbb{Z}}c_{k}(g)\mathrm{e}^{2\pi\mathrm{i}kx}\approx\sum_{k \in\mathcal{C}}c_{k}(g)\mathrm{e}^{2\pi\mathrm{i}kx},\quad c_{k}(g)=\int_{- \frac{1}{2}}^{\frac{1}{2}}g(x)\mathrm{e}^{-2\pi\mathrm{i}kx}\mathrm{d}x$ $$

for $\mathcal{C}=\{\frac{-N_{\mathrm{ft}}}{2},\dots,\frac{N_{\mathrm{ft}}}{2}-1\}$ with some even $N_{\mathrm{ft}}\in\mathbb{N}$.
To compute the Fourier coefficients $c_{k}(g)$, we employ that $g(x)\approx 0$ for $|x|>\frac{1}{2}$ such that Poisson’s summation formula (see, e.g., Plonka et al., , Thm 2.26) implies

$$ $c_{k}(g)\approx c_{k}\left(\sum_{k\in\mathbb{Z}}g(x+k)\right)=\mathcal{F}_{1}[ g](k),$ $$

where $\mathcal{F}_{1}[g](\omega)$ is the Fourier transform equation 12.

In our experiments, we choose $N_{\mathrm{ft}}=128$ for the Gauss, $N_{\mathrm{ft}}=512$ for the Matérn and $N_{\mathrm{ft}}=1024$ for the Laplace kernel.
Note that the coefficients $c_{k}(g)$ do not depend on the input points and need to be computed only once for different choices of $\xi$.
The function $\mathcal{F}_{1}[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

$$ $t_{m}\approx\sum_{n=1}^{N}\sum_{k\in\mathcal{C}}w_{n}c_{k}(g)\mathrm{e}^{2\pi \mathrm{i}k(y_{m,\xi}-x_{n,\xi})}=\sum_{k\in\mathcal{C}}c_{k}(g)\mathrm{e}^{2 \pi\mathrm{i}ky_{m,\xi}}\underbrace{\sum_{n=1}^{N}w_{n}\mathrm{e}^{-2\pi \mathrm{i}kx_{n,\xi}}}_{\eqqcolon\hat{w}_{k}}.$ $$

The computation of the second sum $\hat{w}_{k}$ is the adjoint discrete Fourier transform of the vector $w=(w_{1},...,w_{N})$ at the non-equispaced knots $(-x_{1,\xi},...,-x_{N,\xi})$. Afterward, the computation of the vector $t=(t_{1},...,t_{M})$ is the Fourier transform of the vector $(c_{k}(g)\hat{w}_{k})_{k\in\mathcal{C}}$ at the non-equispaced knots $(-y_{1,\xi},...,-y_{M,\xi})$. These Fourier transforms can be computed by the NFFT in time complexity $\mathcal{O}(N+N_{\mathrm{ft}}\log N_{\mathrm{ft}})$ and $\mathcal{O}(M+N_{\mathrm{ft}}\log N_{\mathrm{ft}})$ leading to an overall complexity of $\mathcal{O}(N+M+N_{\mathrm{ft}}\log N_{\mathrm{ft}})$.

## Appendix J Additional Numerical Results

### J.1 Additional Plots and Tables for Section 4.2

In the following, we redo the experiments from Figure 1 and Table 1 and vary some parameters. More precisely, we redo it 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 1 and Table 1 with the negative distance kernel. The results are given in Figure 4 and Table 3. We can see that the advantage of QMC slicing is not as large as for smooth kernels, which is expected considering the theoretical results from Section 3. In particular, the spherical function equation 11 is not in $H^{d/2}(\mathbb{S}^{d-1})$ if $d\geq 3/2$, so the assumptions of the bound equation 9 are not fulfilled. Nevertheless, QMC slicing is still significantly more 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 theorem does not apply.

##### Other Length Scales of the Kernels

We redo the experiment from Figure 1 and Table 1 with scale factors $s=\frac{1}{2}$ and $s=2$ of the kernel parameter. The results are given in Figure 5 and Table 4 for $s=\frac{1}{2}$ and in Figure 6 and Table 5 for $s=2$.
We observe that the advantage of QMC is more significant of for larger scale factors, which is expected since the function $\xi\mapsto f(|\langle\xi,x\rangle|)$ is more regular for larger $s$ than for smaller $s$.

##### Higher Dimensions

Finally, we do the same experiment as in Figure 1 and Table 1 for the higher dimension $d=200$. Here, we use the negative distance kernel, the Matérn kernel with $\nu=3+\frac{1}{2}$ and the Gauss kernel, where the parameters are chosen by the median rule with scale factor $\gamma=1$.
The results are given in Figure 7 and Table 6.
The advantage of QMC is less pronounced in such high dimensions, but still visible.

### J.2 Additional Results for Section 4.3

##### Negative Distance Kernel

We redo the experiment from Section 4.3 for
the negative distance kernel.
The results are given in Figure 8. We can see that QMC Fourier slicing outperforms standard slicing clearly in all cases.
Note that RFF based methods are not available for the negative distance kernel since it is not positive definite and Bochner’s theorem does not apply.

##### Higher Dimensions

We run the same experiment as in Section 4.3 on the MNIST and FashionMNIST 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

Finally, we want to demonstrate the advantage of our method in a GPU-comparison with a large number of data points.
As a test dataset we concatenate the MNIST and FashionMNIST in all eight orientations arising rotating and mirroring the images and reduce the dimension via PCA 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., ) as a baseline. The results are given in Figure 10. Even though the RFF-based methods parallelize a bit better on the GPU than the fast Fourier summations, the conclusions are mainly the same as 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 better than RFF-based methods.

Figure: Figure 4: Loglog plot of the approximation error in equation 4 versus the number $P$ of projections for the negative distance kernel. The results are averaged over $50$ realizations of ${\bm{\xi}}^{P}$ and $1000$ realizations of $x$. We fit a regression line in the loglog plot for each method to estimate the convergence rate, see also Table 3.
Refer to caption: x26.png

Figure: Figure 5: Loglog plot of the approximation error in equation 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 ${\bm{\xi}}^{P}$ and $1000$ realizations of $x$. The kernel parameters are set by the median rule with scale factor $\gamma=\frac{1}{2}$. We fit a regression line in the loglog plot for each method to estimate the convergence rate, see also Table 4.
Refer to caption: x29.png

Figure: Figure 6: Loglog plot of the approximation error in equation 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 ${\bm{\xi}}^{P}$ and $1000$ realizations of $x$. The kernel parameters are set by the median rule with scale factor $\gamma=2$. We fit a regression line in the loglog plot for each method to estimate the convergence rate, see also Table 5.
Refer to caption: x41.png

Figure: Figure 7: Loglog plot of the approximation error in equation 4 versus the number $P$ of projections for $d=200$. The results are averaged over $50$ realizations of ${\bm{\xi}}^{P}$ and $1000$ realizations of $x$. The kernel parameters are chosen by the median rule with scale factor $\gamma=1$. We fit a regression line in the loglog plot for each method to estimate the convergence rate, see also Table 6.
Refer to caption: x53.png

Figure: Figure 8: Loglog plot of the relative $L^{1}$ approximation error versus computation time for computing the kernel summations equation 1 with the negative distance kernel and different methods and datasets. MNIST and FashionMNIST are reduced to dimension $d=20$ and $d=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\cdot 2^{k}$ slices for $k=0,...,9$.
Refer to caption: x56.png

Figure: Figure 9: Loglog plot of the relative $L^{1}$ approximation error versus computation time for computing the kernel summations equation 1 with different kernels and methods on the MNIST and FashionMNIST dataset without dimension reduction. We run each method $10$ times. The shaded area indicates the standard deviation of the error. For the slicing method, we use $P=10\cdot 2^{k}$ slices for $k=0,...,9$. In order to obtain similar computation times, we run RFF and ORF with $D=2P$ features.
Refer to caption: x59.png

Figure: Figure 10: Loglog plot of the relative $L^{1}$ approximation error versus GPU computation time for computing the kernel summations equation 1 with different kernels on a 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, we use $P=10\cdot 2^{k}$ slices for $k=0,...,9$. In order to obtain similar computation times, we run RFF and ORF with $D=4P$ features. Since PyKeOps computes the exact kernel sum, the computation time is indicated by a vertical line.
Refer to caption: x69.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1° Can you comment on the interest of slicing in Optimal Transport versus in kernel summation ? 
2° It would be interesting to describe the differences point by point of this work and those of Hertrich (2024).
3° Could you clearly describe the analytical complexity in time of each QMC scheme ?
4° Can you find and experiment an example in Machine learning where the sliced kernels can make a difference ? 
5° is it possible to run the same tests on larger dimensional data ?

After rebuttal, I've increased my score (most of the concerns have been addressed except those concerning the experimental results.
soundness: 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.
- 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)).
- 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.
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.
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.

- Limited contribution given previous works:
The most important weakness of the paper is its incrementality with respect to the work of Hertrich (2024) well cited in the paper. 
The 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.
- No applicabe scale for the spherical deisgn: Unfortunately, the spherical design studied in depth here cannot be reasonably computed for  realistic values of d. 
So 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.

In 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.

Moreover 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.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: 1.How scalable is the QMC slicing method for very high-dimensional data (e.g., 100+ dimensions)?
2.Does the QMC slicing approach generalize to non-standard kernels often used in machine learning?
3.How do theoretical error bounds hold up on real-world datasets with noise and variability?
4.Can the authors compare the computational efficiency of QMC slicing to Random Fourier Features and other fast summation methods?
5.What are the practical trade-offs in terms of error convergence rate with QMC slicing?
6.Are there specific use cases or domains where this approach is expected to outperform current alternatives?
soundness: 2
strengths: 1.The QMC slicing method introduces an approach to fast kernel summation.
2.The methodology is well-documented, with explanations of error bounds and smoothness results.
3.QMC slicing has the potential to improve efficiency in kernel-based methods.
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.
weaknesses: 1.The adaptability of the theoretical error bounds to various types of data and kernels remains unclear.
2.While effective, the QMC slicing method lacks clear differentiation from existing fast summation methods, such as Random Fourier Features (RFF).
3.Experiments are focused on limited, synthetic datasets, raising concerns about the method's generalizability.
4.Some parts of the theory section are overly condensed, especially around error bounds and smoothness assumptions.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
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.
soundness: 3
strengths: - The theoretical analysis in the paper seemingly technically sound, though I couldn't follow the proofs.
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.
weaknesses: - Overall, the paper is difficult to follow those who are not familiar with the topic. Readability can be improved. 

- Empirical evaluation is only on purely kernel sum calculations. Showing benefit on actual computation of some learning model would have been convincing.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
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?
2. 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?
soundness: 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.
2. 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.
3. 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.
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.
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.
2. 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.

Meta-review: The paper investigates the problem of fast computation of large kernel sums via the method of slicing, which relies on random projections to one-dimensional subspaces and fast Fourier summation. In particular, the authors prove bounds for the slicing error and propose a quasi-Monte Carlo (QMC) approach for selecting the projections based on spherical quadrature rules.

Reviewers generally appreciate the proposed method, its theoretical analysis, and potential in improving the efficiency of kernel-based methods.

Decision: Accept (Poster)