Fast Summation of Radial Kernels via QMC Slicing

We would like to thank the reviewer for the evaluation of our paper. Please find below the detailed answers to the comments and questions.

• Regarding Weakness 1: In our numerical experiments, the distance QMC designs perform best among the (computable) QMC rules. The only exception is when the (provably optimal) spherical designs are available, which are only computable in $d\le4$. This is now written clearly in Section 4.2. For high dimension $d$ and few directions $P<d$, we showed in Appendix H that orthogonal points and distance QMC designs coincide, but the latter yield better errors for more directions $P>d$.

• Regarding Weakness 2: One possible explanation for this gap could be that our theoretical bounds on the error of QMC designs is based on worst-case error estimates in the Sobolev space $H^s(\mathbb S^{d-1})$. However, our spherical function $g_x$ from eqt. (11) has a very special form. In particular, it is symmetric, i.e., $g_x(\xi)=g_x(-\xi)$, such that the integral over the sphere can be replaced by the integral over an arbitrary hemisphere. If we choose this hemisphere properly, then the function $g_x$ is infinitely differentiable on it for all considered kernels, since the non-smoothness in our examples only appears at the equator $\langle \xi, x\rangle=0$. We have added a discussion in Appendix K.4 and referenced it in the outlook and limitations in Section 5.

• Regarding Question 1: We added an application of computing MMD flows with CIFAR10 (dimension $d=3072$) in the new Appendix K.3, which was used for image generation in the literature. Also in this case, QMC slicing performs better than normal slicing with random directions. RFF-based methods are not applicable here, since the considered negative distance kernel is not positive definite.

• Regarding Question 2: The sigmoid kernel is out of the scope of this paper, as it is not a radial one. We included experiments for the thin plate spline kernel in Figure 8, observing a significant improvement of QMC slicing over slicing also here.

Many thanks for your review! We answer to each of your comments separately.

Adaptability of the theoretical error bounds (addresses weakness 1): Our error bounds apply to any kind of data in Euclidean space and to any positive definite radial kernel, as well as the Riesz kernel and the thin plate spline. We now mention this in Section 1, Contributions.

Comparison to RFF (addresses weakness 2 and question 4): We numerically compare to RFF and related methods (like ORF and QMC-RFF) in our paper (Fig. 1, 2, 5, 6, 7, 8, 9, 10) and demonstrate in all cases a significant advantage of QMC slicing. From a theoretic side, we have outlined in detail (Appendix G), how the slicing method and RFF are related. Some fundamental differences include

• RFFs rely on Bochner's theorem and are consequently only applicable to positive definite kernels. We presented several examples in the paper, where this assumption is violated. In particular, they include the negative distance kernel (Fig. 4, 7, 8, 9, 10), which is widely used in the energy distance (Szeleky 2002). We additionally added the thin plate spline kernel as an example in Figure 8 and added an application to MMD flows in Appendix K.3. In both cases, RFF are not applicable, since the kernel is not positive definite.

• RFF integrate over the measure from Bochner's theorem, while slicing integrates over the unit sphere. In particular, we can exploit QMC designs and quadrature rules on the sphere for QMC slicing, which is a well studied problem from numerical analysis. In contrast, methods combining QMC with RFF always rely on a transformation of the measure from Bochner's theorem to the uniform distribution on the unit cube, which is only possible in very restrictive examples. Due to this difference, our QMC slicing significantly outperforms RFF and related methods like orthogonal Fourier features and QMC RFF (comparisons in Fig. 1, 2, 5, 6, 7, 8, 9, 10).

• We added a theoretical analysis of the complexity of QMC-Fourier-Slicing and RFF for the Gauss kernel in the new Proposition 5 in Appendix J.

Readability (addresses weakness 4): We streamlined the formulation of Theorem 1 and made several minor improvements to increase the readability. Taking into account the page limit, it seems unfeasible to add longer explanations to the main text without removing other content.

Applications (addresses weakness 3 and question 3): We added an application of computing an MMD flow in the new Appendix K.3, where we show a clear advantage of QMC slicing. Additionally, we would like to highlight that, typically, sensitivity with respect to noise is a property of the specific application and not of how the kernel sums are computed.

Higher dimensions (addresses questions 1 and 6): We added an application to MMD flows which acts on the CIFAR10 dataset ($d=3072$) in Appendix K.3. Here, RFF is not applicable because the negative distance kernel is not positive definite, and a direct computation takes considerably more time. Additionally, note that our paper already contained examples in 100+ dimensions: Fig. 9 considers (Fashion)MNIST with $d=784$.

Other kernels (addresses question 2): In Fig. 8, we have added the thin plate spline kernel. Note that it is not positive definite such that RFF-based methods are not applicable. Also in this example, we observe a clear advantage of QMC slicing.

Thank you for your review. Please find below our answers to your comments.

• Regarding readability (weakness 1): We increased the readability at several places, e.g., in making the notation more consistent at some places (for example, we write $g$ instead of $f$ now in Definition 2, in accordance with the notation $g_x(\xi)=f(|(\xi,x)|)$ in (10)), and streamlining the formulation of Theorem 1. Changes are highlighted in blue.

• Regarding applications (weakness 2): We have added a specific use-case of the fast kernel summation for minimizing an MMD functional by simulating a gradient flow in the new Appendix K.3. This was used in several papers for generative modelling and other applications. We note that the exact computation of the kernel sums is computationally intractable here. Also in this case, we observe a significant advantage of the QMC slicing approach.

• Regarding the questions about the minimization of (10), which is eqt (9) in the revised version: We do not consider the computation of the QMC designs as part of the computation time in Fig. 2, since it only depends on the dimension, the number of projections and the Sobolev smoothness. In particular, it is independent of the data points and the kernel such that an a-priori estimation seems reasonable. For $P>d$, we minimize (10) numerically by the Adam optimizer. The computation time of this procedure indeed heavily depends on $P$ and $d$ and ranges from a couple of seconds (in $d=3$) up to one hour (for $P=5120$ and $d=784$), but again, this computation is completely independent of the problem itself. For $P\leq d$ (and $s=\frac{d}{2}$), any orthogonal system is a minimizer in (10) as outlined in Appendix H such that the computation can be done in fractions of a second.

Differences to (Hertrich, 2024)

(addresses comment 2 of the weaknesses part and question 2)

We want to highlight two significant contributions of our paper above the mentioned paper.

First, our error bounds for the (non-QMC) slicing are significantly stronger than in (Hertrich, 2024). We stress here that the issue was not about confirming the $\mathcal{O}(P^{-1})$ rate for the MSE (i.e. $\mathcal{O}(1/\sqrt{P})$ for the RMSE/absolute error) of standard slicing, but to clarify how the constant depends on the dimension of the data space.

• For positive definite kernels and Riesz kernels, our error bounds are dimension-independent. In contrast, (Hertrich, 2024) only proves a dimension-independent absolute error bound for the Gauss kernel and conjectures that this should be true for the Laplace and Matern kernel as well. For the Riesz kernel, the error bound in (Hertrich, 2024) depends on the dimension by $\mathcal O(\sqrt{d})$.

• For the thin plate spline kernel, there did not exist an error bound at all before.

• For kernels with analytic basis functions, for the Riesz kernel and for the thin plate spline, we exactly compute the variance, which implies that there do not exist any tighter estimations of the error. For the Gauss kernel, this yields an improvement of the absolute error bound given in (Hertrich, 2024), as we observe that this bound actually decays proportional to $\|x\|^{-2}$ for $x\to\infty$.

Second, we exploit QMC sequences to improve the error rate of $\mathcal O(1/\sqrt{P})$. We prove smoothness results to ensure the applicability and demonstrate by several numerical examples that QMC slicing, most notably the proposed distance QMC designs, significantly outperforms non-QMC slicing and RFF-based methods.

We added the point-by-point comparison of our error bounds to the ones in (Hertrich, 2024) in Appendix D.1.

Slicing in OT vs kernel summation

(addresses question 1)

In the optimal transport literature, the sliced Wasserstein distance appeared since it is faster to compute than the standard Wasserstein distance. However, it is important to note that the sliced Wasserstein distance is not the same as the Wasserstein distance and its properties fundamentally differ (different sampling complexity, not a geodesic space, the equivalence estimate between Wasserstein and sliced Wasserstein depends exponentially on the dimension). In contrast, in slicing of kernels, the one-dimensional kernels are constructed such that the sliced kernel summation is an unbiased and consistent estimator of the non-sliced one, and the oscillations are independent of the dimension in many cases, as we proved.

We added the sliced Wasserstein distance to the related work section in the introduction.

Complexity of QMC Slicing

(addresses question 3)

We describe in detail the computational complexity of the RFF and Fourier slicing in the new Appendix J, and comment on this in Section 3.2. Since both methods have different parameters, we also compared the complexity to achieve a certain accuracy level $\varepsilon$ for the case of the Gauss kernel in Proposition 5. There we see that the asymptotic complexity of the QMC Fourier slicing is indeed lower.

We would like to thank the reviewer for the evaluation of our paper. Please find below our detailed answers to the raised comments. Before that, we want to clarify one important misunderstanding regarding the used QMC sequences, which appeared in the summary of the review and in comment 3 of the weaknesses part.

• The proposed distance QMC designs are applicable for any dimension $d$, where the minimization problem (9) is tractable, which includes all examples that we considered. We now highlight the use of the distance QMC designs in the "Contributions" in Section 1. They are obtained by minimizing eqt (9) and have a guaranteed convergence rate of $\mathcal O(P^{s/{d-1}})$ for $s\in(\frac{d-1}{2},\frac{d+1}{2})$. However, while the advantage of this convergence rate degrades for large $d$, we numerically observe much better convergence rates. This gap between the theoretical guarantees and the numerical observations is discussed in the limitations and, following the comments of the other reviewers, we have added possible reasons and perspectives to overcome it in Appendix K.4.

• Whether Sobol sequences projected to the sphere are QMC designs in the sense of Definition 2 is not clear at all. They were initially built as low discrepancy sequences on $[0,1]^d$. Whether they remain QMC designs after the transformation onto the sphere is unclear (we clarified that in Section 4.1). However, in all our numerical experiments, the Sobol sequence performs better than iid uniformly distributed slices. At the same time, the Sobol sequence always have a larger error than distance QMC designs. Therefore, the distance QMC designs are used in Section 4.3.

• In the special case that we can compute spherical designs (which is indeed a hard problem for $d>4$, for $d=4$, they can be computed numerically, see Womersley, 2018), we even get an exponential error rate with a smooth kernel such as the Gaussian.

Regarding the motivation, application and high dimensions

(addresses comments 1 and 4 of the weaknesses part and questions 4 and 5)

We disagree with the reviewer that the complexity reduction of kernel summations in the "era of distributed computing is not a crucial obstacle". We added an example of computing an MMD gradient flow, which was used for generative modelling by several papers. In particular, this includes a sequential evaluation of many large kernel sums (similar to an online-learning setting). While for (QMC-)slicing the computation of one step takes between 0.2 and 0.3 seconds, the exact evaluation using the high-performance library PyKeOps takes one hour. Considering that we run this gradient flow for 50000 steps, this clearly demonstrates the massive need of fast kernel summation methods. Our evaluation also shows the significant advantage of using a QMC sequence instead of the standard slicing approach.

Regarding high dimensions: We do not explicitly motivate our method by high data dimensions $d$, but by the case where the number of data points ($M$ and $N$) is large. We clearly write in the introduction that we mainly target problems in dimension up to 100 or 200, and we strongly believe that enough machine learning problems fall into this dimension range to consider the problem as interesting. On the other side, we agree with the reviewer that high dimensions appear in several applications and therefore consider also some higher-dimensional examples. The added application on MMD gradient flows now works in more than 3000 dimensions. Finally, we stress that common frameworks for the brute-force computation of large kernel sums have strict limitations on the dimensions. For instance, the widely used package KeOps (Charlier et al, 2021, we compared with its Python interface PyKeOps) states in its documentation that the performance significantly drops for dimensions $d>100$ due to register spilling, see the paragraph limitations in https://www.kernel-operations.io/keops/introduction/why_using_keops.html (we also observe this in our experiments with CIFAR10).