Repulsive Monte Carlo on the sphere for the sliced Wasserstein distance
Authors: Vladimir Petrovic
Rémi Bardenet
Agnes Desolneux
Venue: Accepted by TMLR ()

Abstract: In this paper, we consider the problem of computing the integral of a function on the unit sphere, in any dimension, using Monte Carlo methods. Although the methods we present are general, our guiding thread is the sliced Wasserstein distance between two measures on $\mathbb{R}^d$, which is precisely an integral of the $d$-dimensional sphere. The sliced Wasserstein distance (SW) has gained momentum in machine learning either as a proxy to the less computationally tractable Wasserstein distance, or as a distance in its own right, due in particular to its built-in alleviation of the curse of dimensionality. There has been recent numerical benchmarks of quadratures for the sliced Wasserstein, and our viewpoint differs in that we concentrate on quadratures where the nodes are repulsive, i.e. negatively dependent.  Indeed, negative dependence can bring variance reduction when the quadrature is adapted to the integration task. Our first contribution is to extract and motivate quadratures from the recent literature on determinantal point processes (DPPs) and repelled point processes, as well as repulsive quadratures from the literature specific to the sliced Wasserstein distance.  We then numerically benchmark these quadratures. Moreover, we analyze the variance of the UnifOrtho estimator, an orthogonal Monte Carlo estimator. Our analysis sheds light on UnifOrtho's success for the estimation of the sliced Wasserstein in large dimensions, as well as counterexamples from the literature. Our final recommendation for the computation of the sliced Wasserstein distance is to use randomized quasi-Monte Carlo in low dimensions and UnifOrtho in large dimensions. DPP-based quadratures only shine when quasi-Monte Carlo also does, while repelled quadratures show moderate variance reduction in general, but more theoretical effort is needed to make them robust.

Published in Transactions on Machine Learning Research (02/2026)
Repulsive Monte Carlo On The Sphere For The Sliced Wasserstein Distance
Vladimir Petrovic
blcdamii-perrooee@unii--llllf-f
Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRIStAL, 59651 Villeneuve d'Ascq, France
Rémi Bardenet
remi.bardennttuunii-illlleef
Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRIStAL, 59651 Villeneuve d'Ascq, France
Agnès Desolneux
agnes. Sgnnesdeessheeuuquqersppsrssssss
Univ. Paris-Saclay, ENS Paris-Saclay, CNRS, Centre Borelli UMR 9010, 91190 Gif-sur-Yvette, France
Reviewed on OpenReview: https: // openreview. met/ forum? id= JScTmB6Ehu
Abstract
In this paper, we consider the problem of computing the integral of a function on the unit sphere, in any dimension, using Monte Carlo methods. Although the methods we present are general, our guiding thread is the sliced Wasserstein distance between two measures on Rd, which is precisely an integral of the d-dimensional sphere. The sliced Wasserstein distance (SW) has gained momentum in machine learning either as proxy to the less computationally tractable Wasserstein distance, or distance in its own right, due in particular to its built-in alleviation of the curse of dimensionality. There has been recent numerical benchmarks of quadratures for the sliced Wasserstein (Sisouk et al., 2025), and our viewpoint differs in that we concentrate on quadratures where the nodes are repulsive, i.e. negatively dependent. Indeed, negative dependence can bring variance reduction when the quadrature is adapted to the integration task. Our first contribution is to extract and mo- tivate quadratures from the recent literature on determinantal point processes (DPPs) and repelled point processes, well repulsive quadratures from the literature specific the sliced Wasserstein distance. We then numerically benchmark these quadratures. Moreover, we analyze the variance of the UnifOrtho estimator, an orthogonal Monte Carlo estimator introduced by Rowland et al. (2019). Our analysis sheds light on UnifOrtho's success for the estimation of the sliced Wasserstein large dimensions, as well as counterexamples from the literature. Our final recommendation for the computation of the sliced Wasserstein dis- tance to use randomized quasi-Monte Carlo in low dimensions and UnifOrtho in large dimensions. DPP-based quadratures only shine when quasi-Monte Carlo also does, while repelled quadratures show moderate variance reduction in general, but more theoretical effort is needed to make them robust.
1 Introduction
In Monte Carlo integration, introducing repulsion between the points at which the integrand is evaluated can bring significant variance reduction. In RR determinantal point processes (DPP) have for example been shown to yield central limit theorem with improved convergence rate over classical Monte Carlo, for compactly supported integrands (Bardenet & Hardy, 2020; Coeurjolly et al., 2021). In the same vein, even a modicum of negative dependence can reduce variance, e.g. applying a single step a gradient descent aimed at minimizing the Coulomb energy between the quadrature nodes (Hawat al., 2023). Beyond Euclidean spaces, Monte Carlo methods with DPPs have been considered over selected manifolds (Berman, 2024; Lemoine & Bardenet, 2024). One natural manifold to look at is the sphere Sd-1
Rd, however,

---

Published in Transactions on Machine Learning Research (02/2026)
beyond the case of S treated in Berman (2024), it is not yet clear whether DPPs and similar randomized quadratures with negative dependence can be a practical asset.
In machine learning, the problem of integrating over gd-1 naturally arises in recent applications of optimal transport. A central object in optimal transport is the so-called Wasserstein distance, an intuitive distance between probability measures with a host of theoretical properties (Peyré & Cuturi, 2018). On the negative side, numerically evaluating the Wasserstein distance between two measures typically starts with replacing these two measures by i..d. realizations, but the quality of the approximation rapidly degrades with the di- mension (Fournier& & Guillin, 2015). Moreover, even between two discrete distributions with M atoms each, the cost of a generic algorithm to compute the Wasserstein distance scales as M3 log(M), which becomes intractable for large M (Peyré & Cuturi, 2018). This has led to research on alternatives to the Wasserstein distance, one of which is the sliced Wasserstein distance (SW).
The sliced Wasserstein distance finds its roots in one-dimensional optimal transport (Bonnotte, 2013). The cost of computing the Wasserstein distance between two discrete distributions with their atoms on a line essentially boils down to sorting the abscissa of the atoms. In higher dimension d, the idea is hence to look at the projection of our discrete measures on given direction, compute the Wasserstein distance between these projected point clouds, and integrate the results along all possible directions. The corresponding quantity is an integral over the sphere, the integrand being a one-dimensional Wasserstein distribution, that defines a metric over the space of probability measures called the sliced Wasserstein distance. The SW distance preserves the main topological properties of the Wasserstein distance, while holding the promise to solve the aforementioned curse of dimensionality and tractability issues (Bayraktar & Guo, 2021; Nadjahi, 2021).
The SW distance has found many applications in machine learning, in gradient descent (Bonet et al., 2022), barycenter computation (Bonneel et al., 2015), generative models (Deshpande et al., 2018; Liutkus et al., 2019) or kernel methods (Kolouri et al., 2016). It also finds applications, for instance, in texture synthesis (Heitz et al., 2021), in regularizing autoencoders (Kolouri et al., 2019) or differential privacy (Rakotoma- monjy & Liva, 2021). More recently, it has also been applied to gradient descent on manifolds (Bonet et al., 2025). Finally, the SW distance has also been used as a proxy to the Wasserstein distance, when comparing the output of different sampling algorithms (Linhart et al., 2024). All the aforementioned applications use repeated evaluations of the sliced Wasserstein distance, e.g. between a reference measure and an iteratively refined approximation. Errors in the evaluations of the SW distance can thus deterioriate the performance of the underlying ML tasks, and it is natural to try to understand and reduce these evaluation errors.
When computing the sliced Wasserstein distance, two kinds of approximations are actually being made. The first one derives from the fact that one usually has access only to a finite sample from the underlying proba- bility distributions; this is known as the sample complexity problem (Manole et al., 2022). The second kind of error comes from the formulation itself of the sliced Wasserstein distance between two finitely supported measures as an integral over the sphere, which has to be approximated in general. The present paper focuses on that second kind of error, which is intuitively independent of the initial ML task. One typically relies on Monte Carlo algorithms on the sphere to get an estimate of the desired quantity. Although the cost of eval- uating the integrand is relatively cheap, stacking up a large number N of evaluations on as many directions on the sphere can still become computationally heavy, and the slow decay in 1-1/2 of the error of crude Monte Carlo integration will typically require such large N (Robert & Casella, 2004).
Several Monte Carlo methods have already been investigated to solve the integration task inherent to com- puting the SW distance. In particular, while we were working on this manuscript, a survey has appeared (Sisouk et al., 2025). Their conclusions are that for d E {2,3), quasi-Monte Carlo methods prevail, while in higher dimension (typically above d = 20), the so-called orthogonal Monte Carlo method (Rowland et al., 2019; Lin et al., 2020) is both more efficient than crude Monte Carlo and computationally cheap enough to be practical in ML applications. In the intermediate range, they do not provide clear guidelines but rather en- courage the reader to experiment. Besides also reviewing existing Monte Carlo methods for SW estimation, our contributions are twofold. First, we introduce and benchmark five randomized quadratures that have not yet been used to estimate the sliced Wasserstein distance. One of these is a natural importance sampling baseline. The four others are joint distributions with negative dependence that we draw and sometimes
2

---

Published in Transactions on Machine Learning Research (02/2026)
mildly adapt from the recent literature on repulsive Monte Carlo methods. Some of the resulting estimators already provably enjoy faster decaying variance than i.i.d. quadratures. To our knowledge, when considering the sliced Wasserstein distance, this has only been achieved by the estimator from Leluc et al. (2024). On top of the interest of computing the sliced Wasserstein distance, our numerical investigations are also meant to help us identify repulsive point processes that are useful for Monte Carlo integration on the sphere. Indeed, proving variance reduction result with negative dependence as in (Bardenet & Hardy, 2020; Hawat al., 2023) can be long and technical, so that it is important that the community focus their mathematical efforts on the most promising candidates. Precisely doing that, i.e. focussing our mathematical efforts on under- standing practically successful estimators, our second main contribution is to compute the variance of an estimator based on orthogonal Monte Carlo (Rowland et al., 2019; Lin et al., 2020). The latter has already been empirically shown to be successful for SW estimation in large dimensions, which our own experiments confirm. Our variance calculation sheds light on the situations where orthogonal Monte Carlo may (or may not) yield variance reduction.
The rest of the paper is organized follows. Section introduces background on repulsive point processes for Monte Carlo integration. Section 3 describes the main properties of the sliced Wasserstein distance, and reviews numerical quadratures that have already been implemented to estimate it. Section 4 presents new candidate estimators for the sliced Wasserstein distance, among which a natural importance sampling scheme and various repulsive point processes adapted to the spherical case. Section 5 presents our deriva- tion of the variance of an orthogonal Monte Carlo estimator known as UnifOrtho. All these methods are empirically evaluated and benchmarked in Section 6, Section 7 concludes the paper. The appendix presents supplementary background on spherical harmonics, additional details on the importance sampling baseline, as well as additional experiments.
2 Repulsive Monte Carlo
Monte Carlo methods are randomized algorithms for quadrature, i.e., numerical integration. The common idea is build linear combinations of a finite number of integrand evaluations well-chosen quadrature nodes (Robert & Casella, 2004). While classical Monte Carlo methods draw their nodes using independent random variables or a Markov chain, many recent works have tried to leverage negative dependence among nodes in Rd to obtain lower mean-square integration errors. One family of approaches focuses on settings where the target is hard to sample from or approximate, and introduces negative dependence into Monte Carlo Markov chains. For instance, the Wang-Landar algorithm an adaptive MCMC algorithm that forces the chain to move away from overexplored regions (Fort et al., 2015). A more recent example is (Ye al., 2020), where the authors introduce in Langevin diffusion term that pushes the state away from the past history of the chain, without altering the limiting distribution. However, so far, these algorithms only come with asymptotic guarantees, and the corresponding rates of convergence are still the same as without the negative dependence. A second family of approaches considers the more restricted setting where the target is easy to sample from or approximate and the dimension is low, and manages to improve on the convergence rate thanks to repulsiveness. In this setting, the bottleneck comes from the number of calls to the integrand, and is already the center of a vast and rich literature, see e.g. (Delyon & Portier, 2016; Leluc et al., 2025). We focus in this paper on this restricted setting where the target is easy to sample from, and we review two recent ways to introduce negative dependence in integration routines in R., determinantal and repelled point processes. We choose these two because they easily adapt to the sphere, we shall see in Section 4.
2.1 Determinantal point processes
Initially invented to model the arrival times of physical fermions in optics (Macchi, 1972), determinantal point processes (DPPs) have seen a recent surge of interest in probability (Soshnikov, 2000; Hough et al., 2006), statistics (Lavancier et al., 2015), and machine learning (Kulesza & Taskar, 2012). Formally, we shall only use projection DPPs, which can be defined as follows.
3

---

Published in Transactions on Machine Learning Research (02/2026)
Definition 1 (Projection DPP) Let X be separable complete metric space, and j be measure on its Borel sets. Let N 2 1, and oo..... ,N -1 be orthonormal functions in T.2()), Let
N ox(x) (x)
K y. H
(1)
Let ((]...,,x)))) be drawn from
1
det (()((K d.((x)...dd())) N!
(2)
Then, we say that the random set X = {XI..... NN} C X has for distribution the projection DPP of kernel K : X X IR and reference measure (4, and we write X ~ Ik)(
First, we note that (2) defines a bona fide probability distribution because K in (1) is a projection kernel, namely the kernel of the projection onto repulsive in the sense that the determinant in (2) favors configurations where the X,S spread evenly across X. Indeed, if K is smooth, two points close to each other correspond to two nearly identical columns in the Gram matrix ((K(s;, (1)),,)))))) and thus a small determinant. Third, DPPs with non-projection kernels can be defined (Hough et al., 2006), but we shall only be concerned by projection kernels in this paper. Finally, on top of having a relatively simple expression, a computational advantage of DPPs that makes them an ideal candidate for summarization tasks is that the chain rule for (2) can be simply expressed using Schur complements (Hough et al., 2006) (Propositin 19]. In more detail, to sample (2), it is enough to sample X from 1/N K (11121)dp(z1)) and for 2,.... N, iteratively sample X₃ from
see e.g. Hough et al. (2006). Second, DPPs are
F() K F 1) J KR (11.2- du(#z),
N k
where K( (1x-1) = K( (1x-1) is short for (K( K( rk-1)), and K
Individual sampling steps in (3) are typically implemented using rejection sampling. In Gautier et al. (2019b), the total number rejections for sampling an Orthogonal Polynomial Ensemble in Rd is estimated to be O(2 log(N)). A realization of this particular DPP in the square |-1, can be observed in Figure lc. In general, much is known on sampling DPPs, exactly or approximately (Gautier, 2020; Barthelmé et al., 2023).
2.2 Monte Carlo integration with DPPs
Monte Carlo methods relying on DPPs with specific kernels have been investigated when X = Rd and the target measure has a density w.t.. the Lebesgue measure, e.g. (Bardenet Hardy, 2020; Mazoyer et al., 2020; Belhadji et al., 2019; 2020). A general conclusion is that for the right choice of kernel, a DPP with cardinality N can integrate smooth functions with mean squared error in o(1/N), thus decaying faster than for classical Monte Carlo methods. For instance, the so-called multivariate orthogonal polynomial ensembles studied in (Bardener & Hardy, 2020) yield mean squared error in 1 1/N1-1/d for integrands that are continuously differentiable. In this paper, we rather consider integration on the sphere Sd-1 C RRd. Using change of variables such as spherical coordinates, is straightforward to adapt e.g. the quadratures proposed by (Bardenet & Hardy, 2020; Mazoyer et al., 2020) for [-1, to gd-1, at the price of an artificial accumulation of points. Closer to our interest for the sphere, Lemoine Bardenet (2024) show that for a compact complex manifold of complex dimension d/2 (and thus dimension when seen as real manifold), the right choice of kernel in (2) yields the faster rate 1/N1+2/d This applies to S?, where the corresponding DPP is called the spherical ensemble; see Section 4.2.2 for more details. However, this result does not easily generalize to gd-1 with d > 3. Still in the particular case d = 3, even finer results are available in (Berman, 2024), actually the first paper to explicitly investigate a DPP for integration on the sphere. Berman (2024) provides a theoretical analysis of the worst-case integration error the spherical ensemble for functions on the sphere in specific Sobolev classes. Finally, for integrands that are smooth enough to belong to a reproducing kernel Hilbert space (RKHS), DPPs (Belhadji et al., 2019) and mixtures of DPPs (Belhadji al., 2020; Belhadji, 2021) have been proven to yield fast-decaying mean squared errors.
4

---

Published in Transactions on Machine Learning Research (02/2026)
2.3 Quadratic-time alternatives to DPPs
Sampling DPP, while polynomial, still be long enough to discard DPPs when the of evaluating the integrand is low. In particular, needs to come up with rejection sampling routines to sample the conditionals (3) in reasonable time; see (Gautier et al., 2019a) for discussion. Aside from DPP sampling, there are ( algorithms that can still achieve a mean squared error decay in o(1/N). For instance, Delyon & Portier (2016) propose variant of importance sampling where the proposal PDF replaced by a kernel density estimator, with a fast error decay. A particularly natural repulsive strategy that does not require strong smoothness assumptions on the integrand is known as repelled point processes (Hawat et al., 2023). The idea is draw a computationally cheap randomized quadrature, and apply one step of gradient descent aimed minimizing the Coulomb energy of the configuration of quadrature nodes, as if they were identically charged particles. The result of such procedure can be observed in Figure 2. It is easy to come up with a similar algorithm for points on the sphere, as we shall do in Section 4. However, the main variance reduction result of Hawat et al. (2023) does not hold for the sphere, and we see our paper as an exploration which algorithms have promising empirical performances to motivate their theoretical study.
(a) 1000 points sampled i..d. uni formly in the square.
(b) The repelled configuration cor- responding to the 1000 points in la.
1000 points sampled from DPP in the square.
Figure 1: Realizations of three point processes
3 The sliced Wasserstein distance
Our motivating application for integration on the sphere is the computation of the sliced Wasserstein (SW) distance between two probability distributions.
Definition Let d N, p > O, and p,v be two probability measures on Rd. The sliced Wasserstein distance between μ and is
(w, -Tr- (W,(0 (nu) d0
1/p
a0)
(4)
where Qg denotes the push-forward measure of by the function as : x E Rd + 0Fr E R, W, is the one-dimensional p-Wasserstein distance, and the integral with respect to the uniform measure on the sphere.
Before discussing its evaluation cost, we give some motivating facts about the SW distance, and refer Nadjahi (2021) for an exhaustive reference.
3.1 Motivating properties
First, for all p 1, SW, metrizes the weak convergence on the space of finite p-moments probability mea- sures (Nadjahi, 2021) [Theorem 3.1]. On compact domains, the topology induced by the sliced Wasserstein
5

---

Published in Transactions on Machine Learning Research (02/2026)
metric is actually equivalent to the one induced by the Wasserstein metric since, if μ4, are supported on B(0, R),
SW,(u, & Wp(u, and (z(,, v) S where Cap a constant only depending on p and d (Bonnotte, PnoosPrroooittionn 5.1.3 and Theorem 5.1.15]. Second, and importantly to train e.g. generative models, it is also possible to perform gradient descent on the SW metric (Nguyen et al., 2024). Formally, if for X E _dx jux denotes the empirical mea- sure supported on the columns of X, then for a discrete measure v, the map X E Rdx M + SW2 E R is Ca (Bonneel et al., 2015) [Theorem 1].
SW, (w
(5)
Focussing now on its practical evaluation, one key advantage of the SW distance over the Wasserstein distance is its dimension-free sample complexity. Formally, for a measure H4, write its p-moment my(p) = f [idu(() and Pas for the empirical measure obtained from M ii.d. draws from p. For p E [1, 00), assume that y and both have finite moment of order g > p. Then
M-1/21 if q > 2p
BRISW,(RM, PM) SW (W,,,))) S / (1G x
M 1-1 log MM 1/p if 2p
M-(0-p)/p4 /44 if q E 2p)
where my (1,, v) my (p)) + mig m//((v) (Nadjahi, 2021) [TTheorem 7.14]. In particular, it is enough focus on evaluating the SW distance (4) between two empirical measures of support of cardinality M. This is an integral over the sphere, which cannot be solved analytically, thus requiring numerical quadrature. Fortunately, the integrand can be efficiently computed: for 2 1, it can be done exactly in time O(M log M+ Md). The Md part comes from the computation of the projection, as each point on which the measure is supported has to be projected onto a line. As for the M log(M) part, it comes from the theory of one- dimensional optimal transport, where one can show that the only computational bottleneck is essentially sorting the atoms of the involved (discrete) measures (Peyré & Cuturi, 2018, Remark 2.30).
The choice of numerical quadrature can be informed by the smoothness of the integrand. We know that the map
f 0 E gd-1 g + WP(0#M, 0gv),
(6)
is Lipschitz (Bayraktar & Guo, 2021) (Propooiioii 2.2], with Lipschitz constant
+ (v))
However, when both measures are discrete, while the map (6) is C outside of a set of measure 0 on the sphere, fails to be globally C² in general. This motivates the use of numerical quadratures that do not make strong smoothness assumptions on the integrand, intuitively dismissing methods that rely on the target integrand being in an RKHS such as (Belhadji et al., 2019; 2020). Moreover, the dependence in Md of the cost of evaluating the integrand -due to computing the projections in the push-forward measures- justifies searching for quadratures with a fast-decaying error if we are to estimate the SW between large datasets in high-dimensional spaces. Repulsive Monte Carlo methods such as DpP-based quadratures and repelled point processes thus seem natural to investigate. Before doing so, we quickly review the existing literature on advanced Monte Carlo techniques for the SW.
3.2 Existing Monte Carlo methods for the sliced Wasserstein distance
Besides the natural i.i.d. sampling on the sphere, several advanced Monte Carlo methods have been proposed that reduce the mean squared error in estimating (4), using either control variates or randomized grids.
3.2.1 Control variates
Control variates is standard variance reduction technique in Monte Carlo integration (Owen, 2013, Chapter X). In nutshell, consider yp; : 1d - R, i 1,...,&, such that J (f) = 0 for all Letting f gd-1 - R be a square-integrable function, and 01.... ., AN be drawn i..d. uniformly on the sphere, consider the ordinary least-squares (OLS) problem

---

Published in Transactions on Machine Learning Research (02/2026)
(f),
argmin RER.
a a- _ 3),9)(01))
To gain intuition, we note that for a fixed B E IR, minimizing the RHS of (7) in a yields the empirical mean of f(0,) - Ej=1 (O:), where i= I,...,N. N. This should in turn be close to f f(0)d since the yjj have integral zero. Optimizing over 3 further reduces the variance of Ig (f) at the cost of introducing small bias. The key decision to be made by the practitioner is the choice of and the control variates PI,...., (p. In particular, as both and N to infinity, if the space spanned by the control variates is large enough to allow reconstructing the integrand f, Portier & Segers (2019) obtain a central limit theorem with a squared error decaying faster than the Monte Carlo rate 1/N. We now present two choices of control variates that have been proposed in the specific case of the sliced Wasserstein integrand (6): the up/low method of Nguyen & Ho (2024) and the spherical harmonics in (Leluc et al., 2024).
Control variates "up" and "low". For two probabilities μ, v on Id with finite first and second moments my: My, En. En, we know (Peyré Cuturi, 2018, Remark 2.9) that the 2-Wasserstein distance satisfies
(2(,,,) v) = ml m,,11 2 + (2,
(8)
where ji, i are the centered versions of μe and u, ie i = (1) computing SW2, Nguyen & Ho (2024) thus suggest taking
# μ, where ((H) E Rd + x - my- When
control variate in (7), with Φi equal to
flow : 0 + (6T (m= - mm))) ā m - m
(9)
Note that flow is centered, and that it will likely have little impact when either * 2 or the target distribu- tions are already centered. In the same spirit, Nguyen & Ho (2024) also propose 1, with &i this time equal to
1 Pup 0 Flow (0) + 0" E D' d (T+(Ew) Ir(E)))
(10)
where the quadratic term on top of Flow upper-bounds the W2 distance between two centered Gaussians -hence the label up- and the remaining term guarantees a null integral, as required for a control variate. This time, the control variate is expected to pick up second-order information. Finally, note that while Plow and Pup use rather crude approximations to the integrand and are limited to the case p = 2, they are both cheap to compute and provide useful baselines.
Spherical Harmonics, Still based on (7), Lelue et al. (2024) rather propose to take 1, P1. $2.... to be spherical harmonics (Yg, k. 2 0, S k S he), ordered in the lexicographic order of (e, k). To wit, Yg = 1 is constant, and, for l 2 I, {Yg, 2 1} form an orthonormal basis of He, the hy-dimensional set of harmonic homogeneous polynomials of degree É restricted to gd-1. We refer to (Leluc et al., 2024, Section 4.1) or our Appendix A.1 for quick self-contained definition of spherical harmonics, but for now it suffices to say that 11, ....... is an orthonormal basis of T.2 (gd-). For fixed N, let 8 = be the number of spherical harmonics of degree at most 2Ln. Note that sn = * The estimator SHCV v) of the SW distance between two probability measures on Rd is then defined to be I in (7).
Leluc et al. (2024) prove that for d 2 2, p 2 1, n,v v having finite p-th moments, and when sn o(N2), so that LN N1/2(d-1)/ /EN for some sequence IN going arbitrarily slowly too when N grows,
JSHCV (14, v) SW v) = Or(EN (1/(6/NN-0/22/2(--1)))
(11)
demonstrating a reduction in the error rate compared to standard Monte Carlo. The whole procedure runs in O(Nw, + N )) where w the time complexity of evaluating f, so that the procedure is quadratic if snv = 2 as prescribed. It is expected that SHCV will be efficient when the integrand (6) appearing in the definition of the SW distance will be well-approximated by polynomials of degree lower than sN,
7

---

Published in Transactions on Machine Learning Research (02/2026)
and that it will outperform the control variates flow and (Pu soon as the degree is large enough, at higher computational price, however. Another caveat that we shall discuss again later is that the complexity estimate ignores the computational time spent evaluating spherical harmonics, which can be prohibitive in large-dimensiona settings; see also Appendix A.1.
3.2.2 Randomized grids
Letting N be the number of evaluations of the integrand (6) that one is willing to spend, and assuming for simplicity that & = N/d an integer, Rowland et al. (2019) propose to take k i.id draws from the Haar measure on the orthogonal group (d)) The columns of these matrices are then marginally uniformly distributed on the sphere, and the average of the integrand (6) over the reunion of these N kd columns is thus an unbiased estimator, called the UnifOrtho estimator in (Rowland al., 2019). Intuitively, since the columns of single Haar draw are orthonormal, they fill the sphere quite evenly, thus justifying our classification as a randomized grid. One could expect some variance reduction coming from this very uniform spread, but there appears to be no such theoretical guarantee so far. Rowland et al. (2019) even exhibit a counterexample of two empirical measures such that UnifOrtho yields a worse (i.e. higher-variance) SW estimator than crude i..dd Monte Carlo on the sphere. We clarify the situation with an explicit derivation of the variance of the UnifOrtho estimator in Section 5.
Quasi-Monte Carlo (QMC; Dick & Pilichshammer, 2010) methods are deterministic quadratures that be thought of as the computationally tractable higher-diieg gigher-dimensional version of a grid. Worst-case guarantees usually involve proving that the quadrature nodes have low discrepancy, and an additional randomization can help obtaín guarantees with more tractable constants. QMC methods for computing the SW distance have been numerically investigated in the dime--dimesiional setting in Nguyen et al. (2024). However, there no known low-discrepancy sequence on gd-1, as soon as d 2 3. An empirically promising alternative (Nguyen et al., 2024; Sisouk al., 2025) is to use the so-called Fekete points as quadrature, a notion from potential theory defined the set of points that minimize particular interaction potential over the sphere. We note however that the construction of Fekete points on the sphere in polynomial time is known to be hard problem, and in dimension 3 is even listed as Smale's 7th problem (Smale, 1998). A more straightforward alternative to low-discrepancy quadratures is to map a low-discrepancy sequence in [0, to gd-1, via some transformation such as using the inverse cumulative function of the normal distribution. Empirical results have been however less encouraging (Nguyen et al., 2024; Sisouk et al., 2025).
In this paper, we will consider randomized QMC benchmark in two and three dimensions, i.e. on S and $2. On S2, we use the generalized spiral points from Rakhmanov et al. (1994), which are easy to draw and have been proven to have low discrepancy, at least asymptotically (Brauchart et al., 2014). To wit, consider 2; (2i 1)/N for 1 iis N and
111 cos ((z), t2.2 1.8 //N mod(2
(12)
The generalized spiral points are the points on the sphere with spherical coordinates (æ,,,, 22.2). Note that the constant 1.8 is chosen arbitrarily, and is used to match the experimental setting of Nguyen et al. (2024). To randomize the quadrature and obtain an unbiased estimator, we simply apply a single uniformly drawn rotation to all points. In the two-dimensional setting, we will also include the regular grid on [- with a random rotation of uniformly drawn angle f ~ ~u[[-
4 New candidate estimators
We propose new estimators for the integral inherent to the sliced Wasserstein distance. Our first proposition is a natural importance sampling baseline, and the rest are repulsive methods: three DPPs, a repelled point process. For the DPPs, we select existing DPPs the probability literature and motivate our selection by applying existing theoretical results to the particular case of the sliced Wasserstein integrand. The novelty there is thus in the application of these DPPs, rather than in the creation of novel kernel or, say, a new central limit theorem. All estimators will be numerically compared in the experimental section.
8

---

Published in Transactions on Machine Learning Research (02/2026)
4.1 An importance sampling baseline
Crude Monte Carlo approximates the integral in (2) using i.i.d. samples from the uniform measure d0 the sphere. Importance sampling consists in rather drawing .1....,0 ON from a measure with density g with respect to d0, and then to define the estimator
K
1 ( - N g(0,)
f(0.)
(13)
It is unbiased by construction, and the choice of the proposal distribution g which minimizes Var(L (f)) is Sopt lffl (Robert & Casella, 2004) [Theorem 3.12]. Since this proposal is not available in practice, several schemes have been proposed to approximate it using part of one's computational budget in evaluations of the integrand. For instance, limiting ourselves to proposal distributions in the parametric family
1 Gomf {1 vmf(--l,,) + vmf[- le le, vins(-le, x)) ((w))xxp(==7()
<> 0, E ES-1
of symmetrized von Mises-Fisher distributions, we spend a fixed fraction (0, of our N evaluations of the integrand to find the PDF g" in Gvuf that minimizes an estimate of the KL. divergence between g and gopt; this is the so-called cross-entropy method (Kroese et al., 2013); see Appendix A.2 for numerical details on how we perform the fit.
4.2 Three determinantal point processes
4.2.1
Orthogonal polynomial ensembles on spherical coordinates.
Representing points on the sphere by their spherical coordinates, we can obtain DPP on the sphere by x [0, rr;; changes of coordinates are C' (C-ddiffeomorphisms and thus preserve DPPs (Lavancier et al., 2015, Proposition A.1.). As DPP baseline, thus blindly follow Bardenet & Hardy (2020), who a projection DPP (Definition 1) with eigenfunctions (dx) in (1) being the products of Legendre polynomials, orthogonal with respect to the uniform distribution. Efficient rejection sampling routines that implement the chain rule (3) are available in the Python package DPPY Gautier et al. (2019b). A central limit theorem for a simple estimator built on such a DPP is available in Bardenet & Hardy (2020), thus potentially helping us obtain asymptotic confidence intervals. However, our integrand (6) not regular enough, nor is compactly supported within the interior of as required in the results of Bardenet & Hardy (2020). Intuitively, we should rather use DPPs that handle both the manifold structure of the sphere and allow for less smooth integrands.
mapping a DPP on X = [0, 2m 2m)d-2 d-2
4.2.2
The spherical ensemble.
In the specific setting 3, another projection DPP over & available in the probability literature, the so-called spherical ensemble. The spherical ensemble comes from random matrix theory, with dedicated sampling algorithm by construction.
Definition (Spherical ensemble, Theorem 3 in Krishnapur (2009)) Let A and be standard i..d. N x complex Gaussian matrices. Consider : S2 [North} + C the stereographic projection (North being the North pole), and Ai...., Ay the eigenvalues of the random matrix A B. Then Sn = ( (A1)...., -1 is DPP with respect the uniform measure do the sphere, of almost sure cardinality N.
[XN)]
This point process is naturally repulsive as can be observed on Figure 2.
There are several results that support using the spherical ensemble for Monte Carlo integration on the sphere. Berman (2024) showed that, akin to quasi- Monse Carlo designs, it has low discrepancy with high probability, thus yielding fast-decaying worst-case integration error for smooth functions. In more Monte Carlo vein, there exist fast central limit theorems for the spherical ensemble under weak smoothness assumptions.

---

Published in Transactions on Machine Learning Research (02/2026)
Indeed, for our integrand (6), Theorem 1 in Rider & Virag (2007) and Theorem 2.5 in Marzo Sánchez et al. (2024) imply
Var
Σ
(121001
+ + -
(14)
VOEN
and
N
1 N
E
(2 - 122500 Law N (0,
( (VTKE12))
(15)
The only smoothness assumption needed if for the variance in (15) to be finite. In our case, this follows from our integrand being Lipschitz continuous, so that it has an almost-evens everywheee bounded gradient by Rademacher's theorem (see e.g. Cheeger (1999); although when both measures are discrete, supported on M points, Rademacher's theorem can be replaced by noting that the integrand is C except on a finite union of great circles). The estimator in (15) has the fastest converging mean mean-square error in results for the estimators of the sliced Wasserstein discussed in this paper, beating the rate in 1/N1-1/2 1/N3/2 associated to the control variates in (11). We thus expect the spherical ensemble dominate Monte Carlo estimators in d 3. A major downside of the spherical ensemble is that it is hard to generalize in higher dimensions; see Beltrán & Etayo (2019) and Lemoine & Bardenet (2024). There is however a close cousin to the spherical ensemble that generalizes to any dimension.
among known
4.2.3 The harmonic ensemble.
Following the formulation given in Marzo Sánchez et al. (2024) and Beltrán et al. (2016), let H4 be the space of homogeneous harmonic polynomials in I& of degree f, restricted to the sphere gd-1, and he dim (He). The harmonic ensemble is the DPP with respect to the uniform measure on the sphere and with kernel
K(x,V) v) = (I+tad11/3 1)/2)
HL p((d-1)/2.(4 1)/2,1 11/2 (/2,))//ddd)))22-)) (x y)
(16)
where a2 ho +... + h2 and pf(d- p( 11/12-2 is Jacobi polynomial (Gautschi, 2004). One can show that it is a projection DPP in the sense of Definition 1, where the eigenfunctions (dxx) are given by spherical harmonics {Y 2 0, < k < ha]; see Appendix A.1. The harmonic ensemble can be sampled using the chain rule (3), although for d = 2, there is a simpler random matrix model, as the harmonic ensemble is known in this particular case as the Circular Unitary Ensemble (CUE), which the law of the eigenvalues of Haar-distribured unitary matrix; see e.g. Remark 4.1.7 in Anderson et al. (2010). A realization of this specific point process on 32 can be observed in Figure 2.
Like the spherical ensemble, a strong motivation for using the harmonic ensemble is the availability of a fast central limit theorem that translates into small asymptotic confidence intervals for Monte Carlo integration. Indeed, letting SN = (01...., Ory } be the harmonic ensemble, Theorem 2.2 in Marzo Sánchez et al. (2024) implies that for our integrand (6),
lim N
Var
Σ (LL)(O) = f
N
L
Moreover,
N
(1)
/ (122(0)d6)
lay
N(O.
|| | |??
N
BESx
where II-IIl g is a specific semi-norm on the Sobolev space H (gd-1) that is equivalent the semi-norm
If1 : J1
xx / [s(x) - f(v)z drdy, f L2,
(17)
7( yd
where 7) is the geodesic distance on gd-1; see Marzo Sánchez et al. (2024). To wit, H (Sd-1) is the function space for which this quantity is finite. The Lipschitz continuity of (/ ensures that s(6) E (dd(111))
10

---

Published in Transactions on Machine Learning Research (02/2026)
that Marzo Sánchez et al. (2024) [Theorem 2.2] applies and gives the aforementioned central limit theorem.
Note that this definition of the harmonic ensemble constrains us to sample a specific number of points, T,, It can be interesting to look at what happens in intermediary levels i.e. to consider incomplete harmonic ensembles. This has been implemented but the runtime becomes quite large when the number of points grows.
4.3 Repelled point processes on the sphere
To further reduce the computational cost of repulsive Monte Carlo compared to DPPs, quadratic-time alter- natives have been considered for integration on R", such as the repelled Poisson process of Hawat et al. (2023) that we recall in Section 2.3. We propose a straightforward adaption the sphere. More precisely, let X be a finite point configuration on the sphere gd-1, and r E X. Define
Fs xx) = Σ
(18)
x y|l*
yE y
which we think of a repulsive force exerted on x by the other points of the configuration. Like Hawat et al. (2023), unless otherwise specified, we take = d. We consider the repelled configuration
cF,, cF,, x()) | c x [I)]
i X
x
X
(19)
where, unlike Hawat et al. (2023), we need to project back onto the sphere. Letting the original configuration be a Poisson point process of intensity p 0, tentatively extending the results of Hawat et al. (2023), we expect the estimator
PPPP n. x (KLL
1 P
Σ (z) xEñ
(20)
to be an unbiased estimator of the sliced Wasserstein distance between μ and v, with reduced variance compared to a sum over X, at least for e > O small enough. Similarly, we expect the same properties to hold if the initial point process set N i.i.d. draws from the uniform measure on the sphere.
Note that in Hawat et al. (2023), choice of independent of f, and proportional to pー1 is suggested. Our empirical findings (see A.5) suggest that this should be the correct magnitude for our in the case d. Note also that the whole procedure only requires the computation of all the pairwise distances and hence runs in O(N2), as it is the case in the Euclidean setting, where N is the number of projection directions to be sampled. Overall, we mainly focus our study to a binomial point process X with N points. It is also possible to apply this repelling step to all the other methods presented here. This leads in various cases to a significant variance decrease at relatively cheap computational cost as we will experimentally show.
5 On the variance of the UnifOrtho estimator
UnifOrtho, as introduced by Rowland et al. (2019) and recalled in Section 3.2.2, is recommended by the recent (Sisouk et al., 2025) for SW estimation in large dimensions. Anticipating on our own experimental results in Section 6, we will recommend as well. However, a theoretical understanding of the variance of the UnifOrtho estimator is lacking, and its proponents even identified cases where the variance might exceed that of crude Monte Carlo estimator based on i.i.d. samples (Rowland et al., 2019). We contribute here new derivation for the variance of the UnifOrtho estimator, which sheds light integrands for which it brings variance reduction. This behavior comes from fundamental properties of the spherical harmonics.
For the sake of clarity, recall from Section 3.2 that the spherical harmonics are defined as follows. For E 1, let He be the set of harmonic homogeneous polynomials of degree & restricted to s and he its dimension. Take Yg to be constant. Then, for 1, take YE, 1 < k S he, to be any orthonormal basis of He. The family of all Yf, 2 0, < k < he, is an orthonormal basis of (Sd-1). For more details about spherical harmonics, Appendix A.1 or (Dai & Xu, 2013).
11

---

Published in Transactions on Machine Learning Research (02/2026)
(a) N 1000 ii.d. points
(b) N = 1000 ii.. points repelled
(c) = 999 points from UnifOrtho
(d) N 1000 points spherical orthog- onal polynomial ensemble
(e) N 1024 points harmonic ensem- (f) N 1000 points spherical ensem- ble
ble
Figure 2: Various point processes over the sphere
Proposition 4 f be a continuous function on gd-1, and (X1l... (Xa) be matrix drawn from the Haar measure on the orthogonal group (d). Let f(e, k) = Jg-- f(x) (x)dr denote the spherical coefficients of f. Then
Var
ΣΣ ((xi) d
1
Var(f(xx))
1 to0
1)- 1-1 {2e
(21)
1=1
E=1
d-1 1 Var Var(f(Xi))
)1-2((16-2(1)) az(-1Hu((f)),
(22)
lu where Hz(f) Σ f(e,k)) ar 1=1
26 +
I( ( 26+1
28 + 1- and A2s
2
Note that in spite of the freedom to choose any basis of Hj when defining the spherical harmonics, the coefficients ue(f) are well-defined, as they are independent the choice these bases.
Proof: Expanding the variance and using the invariance by rotation of the Haar measure yields
1 d
1 s(x))
-EE((X)++++ + d d =(xx((x))((xa) E[f(x))
Var
d i=1
12

---

Published in Transactions on Machine Learning Research (02/2026)
By construction of the Haar measure, conditionally on X1, X₂ follows the uniform measure on S n X (i.e., the d - 2-dimensional Hausdorff measure Hd-2, normalized to have mass 1) (Meckes, 2019, chapter 1.2). In particular,
= -IX(X)))))))] = -If(x;)
(23)
where Ff(u) fea- (w) is the Funk transform of f. Now, combining Theorem 3.4 and Example 3.12 in Rubin (2024) shows that the spherical harmonics are eigenvectors of the Funk transform. More precisely, for all E N and 1 < k < hi, FY 0 and
FY k = --1)
We note in passing that this is analogous to the classical Funk-Hecke formula Dai & Xu (2013) [Theorem 2.9] and comes from the reproducing property of the spherical harmonics kernel
he
Zz(x,y)
hy k(),k)y and reporting into (23) yields 0
for é 22 1. Finally, decomposing f as f f
Σ
Σ
B((()))(Xa)
Σ 1-1) A2e
8=0
Now Ho(f) [[f(X])]?, À = 1, and standard properties of the Gamma function show that Age = 22/22-2. Combining these facts gives the result.
Proposition 4 calls for comments. The first term in both (21) and (22) is the variance of the crude Monte Carlo estimator, and (21) and (22) are two different expressions for the difference in variance between UnifOrtho and that crude Monte Carlo estimator. First, is clear from e.g. (21) that one can get either a decrease an increase in variance from UnifOrtho, depending integrand ʃ. This explains the observed increase in variance in an example of Rowland et al. (2019). To make another more extreme example, note that X2- = =Ld 1, so that
the "energy profile" (Hze(f))EN of the
Var Σ
(2(x.)) =
for all k. In contrast, integrating Y! leads to an increase in variance compared to crude Monte Carlo. Second, we note that the sum in (21) is alternating: each nonpositive term in the sum is followed by a nonnegative term. In P2z some high even frequency l can be responsible for variance increase in (21). When the dimension grows, the generalized Stirling formula yields A2r = ( ), so that only the first terms in either (21) or (22) carry significant weight. The interest of (22) is to show the effect of dimension growth on sequence of integrands with the same spectral profile (fz{(f) throughout dimensions: as a for fixed l decreases as 1/d, nonnegative terms get attenuated more and more, and the variance overall decreases. Third, note that the Funk transform sends all odd-degree spherical harmonics to zero, and in particular UnifOrtho has the same variance as crude Monte Carlo for an odd integrand. The integrand (6) in the sliced Wasserstein distance is even, hence only decomposes onto even harmonics, in coherence with UnifOrtho's success.
2, Aze 1 for all f, so that each term carries the same weight, and single large isolated
6 Experiments
In this section, we numerically illustrate the repulsive Monte Carlo estimators of Section 4. The methods we compare are often referred to using acronyms.
i.i.d. is classical Monte Carlo with i.i.d. uniform points the sphere; it is the default baseline.
13

---

Published in Transactions on Machine Learning Research (02/2026)
ISVMF is short for importance sampling with von-Mises Fischer proposal; see Section 4.1.
UnifOrtho refers to the union of independent Haar-distributed bases introduced in Rowland et al. (2019); see Sections 3.2.2 and
CV up (resp. CV low) is short for Control Variates "up" (resp. "low") as in Nguyen Ho (2024); see Section 3.2.1, Equation 10 (resp. 9).
SHCV is short for Spherical Harmonics control variates Sliced Wasserstein, as introduced by Leluc et al. (2024) and described in Section 3.2.1, Section 3.2.1.
Repelled is described in Section 4.3, specifically in (19), while Repelled SHCV corresponds to using
spherical harmonics control variates (see Section 3.2.1) built on the repelled points. To be more precise, we first sample points according to the repelled point process from (19), to which we then apply the SHCV process from Section 3.2.1.
The three DPPs from Section 4.2 are denoted as OPE for the stereographic projection of the mul- tivariate Jacobi orthogonal polynomial ensemble from Section 4.2.1, Harmonic for the harmonic ensemble from Section 4.2.3, and Spherical for the spherical ensemble from Section 4.2.2. Note that CUE (short for Circular Unitary Ensemble), is the 2-dimensional version of the harmonic Ensemble.
Spherical SHCV, only present when d 3, consists in applying spherical harmonics control variates to the spherical ensemble.
Finally, QMC or Randomised regular grid corresponds to the randomized quasi-Monte Carlo grids in d E (2 described in Section 3.2.2 in Equation 12.
Following Sections 3.2, and 5, we expect witness the following decay rates. If N is the number of projections considered, the classical variance decay rate of i..dd points is 1/N. Some methods enjoy faster than i..dd decay rate, such as the harmonic ensemble from Section 4.2.3 or the SHCV method from Section 3.2.1 at 1/N1+1/(d-1), In the specific case of = 3, the fastest known method is the spherical ensemble, from Section 4.2.2 at a rate 1/N1+2/2 = 1/N². We cannot say much about the decay rate of the other methods, except from UnifOrtho, which still at 1/N but with a different constant as stated in Section 5.
We consider three different experimental settings. The first one is a toy example where we compute the SW distance between two independent Gaussian samples. In order to see how our algorithms behave when comparing more realistic point clouds, we then compute the SW₂ distance between pairs of datasets from a database of sitree-iimensioall point clouds (Chang et al., 2015) used in previous papers on the SW distance (Leluc et al., 2024; Nguyen & Ho, 2024). Finally, to generate a different kind of realistic point clouds, we place ourselves in the position of researcher who wants to compare the outputs of various MCMC algorithms, a task for which the SW has recently been used (Cardoso et al., 2023; Linhart et al., 2024). This time, we focus on SW, rather than SW2, since the former corresponds to worst-case integration error, a natural figure of merit compare MCMC algorithms.
The implementations used are available here httt::://ggttuub.c/////aaapppae
6.1 Gaussian toy example
For any given dimension d, we sample two independent vectors mx, my from (00, Ia), and, independently, two matrices U, V from (0, Taxa). Consider then Σx = U and Σy = VTV. Finally, sample 21,..., (resp. yl, y1,...,yrr) i.i.d. from N(mx, Σx) (resp N (mmy )), and define
M
M
M
M
For each dimension and number of projections, we consider 100 independent realizations of each estimator. In d 2 and d = 10, for SHCV, maximal degree of for the spherical harmonics is fixed, as in the original
14

---

Published in Transactions on Machine Learning Research (02/2026)
0.6
Unbketaa
SBCV
10-3
Pepdaa
Pry-dal SWOV
CUE
-012
10
-0.6
10-4
1000
2000
400
550
700
50C
900
1000
(a) MSE vs computing time for Gaussian toy example, (b) Errors vs number of projections for Gaussian toy example,
2
=
0.5
104
0.2
102
ISVMP
Bipeldd
0.0
Bayil SEDVV
100
-012
10
1400
(c) MSE vs computing time for Gaussian toy example,
(d) Errors vs number of projections for Gaussian toy example,
d 10
d
10
0.6
10
SBEV
10
10:2
-0.2
0..2
10
1.00
22.01
350
2000
7101
1500
1000
11500
1.201
1450
(000
(e) MSE computing time for Gaussian toy example, d 20
(f) Errors number of projections for Gaussian toy example, d = 20
Figure 3: Results for the Gaussian toy example, across 2, 10, 20. The actual value of the 2-sliced Wasser- stein distance is estimated using Monte Carlo integration with 10 projections.
paper (Leluc al., 2024). For = 20, this maximal degree reduced to 2. Empirically, up to 1600 projection points in d 20, the number of control variates corresponding maximal degree of 3 is indeed too large for the estimator to get near the (known) value of sliced Wasserstein. Note that similar phenomenon is
15

---

Published in Transactions on Machine Learning Research (02/2026)
observed in d = 10 on 100 projections 250 projections (see Figure 3d). This is related to the requirement fixed in Equation 11 for the estimator to be consistent.
The results are given in Figure 3, with the left panel showing estimated mean-squared errors vs. computing time, and the right panel showing boxplots of the integration errors. The reference values are computed with a comparatively long Monte Carlo run.
For d 2, Figure 3a highlights that the randomized regular grid far outperforms any other method in terms of MSE. The determinantal point process CUE stands as second, and all the other methods stand in the same range in terms of MSE. Things are different in the 10- and 20- dimensional settings, where the randomized grid and the DPPs do not feature anymore among the leading methods. In d = 10, as per Figures 3c and 3d, the differences between the methods are less sharp, but UnifOrtho dominates, closely followed by CV low, as well as SHCV and Repelled SHCV, once there are enough projections for the linear systems for consistency show. In d 20 dimensions, the only relevant methods seem be UnifOrtho and CV low, which far outperform any other method. These conclusions coherent with the ones presented in Sisouk et al. (2025).
Overall, repulsive methods are among the leading methods in each dimension, but no single repulsive method uniformly dominates: as expected, a randomized grid or a well-chosen DPP are adequate in low dimension, while higher dimensions seem to favor UnifOrtho. Maybe surprisingly we note that repelling the points seems to have only a moderate effect on the MSE. This effect is not even guaranteed to be a decrease in the MSE, and we will investigate this more quantitatively in the Appendix A.5.
6.2 chres-aimennional point clouds
We now consider three-dimensional point clouds Chang et al. (2015) in the Shapenet database. They are configurations of points that cover shapes that range from simple cylinders to planes or benches. We arbi- trarily consider four point clouds from the database, and compute the difference between point clouds #2 and #34, and the distance between point clouds #3 and #35. The point clouds are shown in Figure 4.
(b) Point cloud #34: cylin-
(a) Point cloud #2: table
der
(c) Point cloud #3: sofa
(c) Point cloud #35: chair
Figure 4: Two-dimensional projections the various point clouds used in Section 6.2.
The typical integrand in the SW2 between two such point configurations looks quite different from the toy Gaussian case of Section 6.1. In particular, can be multimodal; see Figure 9. The results of our experiments are shown Figure 5.
The results are comparable with those of the two-dimensional Gaussian toy example of Section 6.1, com- forting the conclusion that in dimension d < 3, for the integrands and the regimes we consider, randomized grids should be the default quadrature: they are both cheap to sample and provide significantly more ac- curate integral estimators than sophisticated Monte Carlo methods such as SHCV or DPPs like the spherical ensemble.
Among the other methods, three seem to be of similar performance: SHCV, repelled SHCV and the spherical ensemble. As the number of projection directions grows, the spherical ensemble gains an edge over the other two methods, in accordance to its faster variance decay (15). A further improvement is obtained by
16

---

Published in Transactions on Machine Learning Research (02/2026)
0.004
Utalertao ug
Rigalled SERDY
0.000
1017
Value
0.000
0.002
-0.044
10
256
576
1024
Number of projections
(a) Errors for the SWa between point clouds #2 and #34.
0.004
0.002
Value
0.000
0.002
-..04
256
11 576
784
1024
Number projections
Errors for the SWa between point clouds #3 and #35.
Figure 5: Boxplots of the errors for three-dimensional point clouds. The boxplots are centered around reference value of the sliced Wasserstein estimated using QMC with 10 points.
combining the spherical ensemble with SHCV, i.e. evaluating the SHCV estimator on spherical ensemble realization rather than i.i.d. points. Finally, we note that ISVMF does not necessarily reduce the MSE of the i..d. estimator: this is likely due to the multimodality of the integrand, which is not reflected in the von Mises-Fischer proposal.
To understand the behavior of the UnifOrtho estimator in the specific case, we used a QMC sequence to estimate the spectral profile (Hze(f)), where f is the integrand the sliced Wasserstein between two point cloud, we refer to Proposition for further details. Figure 6 shows in both cases fast decay of these coefficients, with sharper slope for comparing point clouds #2 and #34. This explains the higher gain of UnifOrtho in this case, seen in Figure 5.
17

---

Published in Transactions on Machine Learning Research (02/2026)
10-4
al C31
10
10
10-1
0
10
Figure 6: Evolution of fe for the two integrands appearing in SW2 in Section 6.2.
6.3 Comparing MCMC kernels
To provide a realistic use case of the SW distance in high and arbitrary dimension, we consider the numerical validation of an MCMC kernel. Numerically assessing that an MCMC kernel targets the expected distribution, as well as comparing MCMC kernels in terms of integration errors with respect to a target distribution, are natural tasks in computational statistics and machine learning. In that context, the sliced Wasserstein between realization of MCMC history and a known target distribution can be used as a figure of merit, done e-g. in Cardoso et al. (2023). To see why, note first that the 1-Wasserstein distance worst-case integration error. Indeed, the classical dual formulation of the W1 distance reads
/ fdy - /
W Wi(, = sup
(24)
f-Lip(h)
where Lip(f) is the Lipschitz constant of the (Lipschitz) function f; see e.g. (Peyré & Cuturi, 2018). Second, SWi can be used as proxy for W1, in the sense of the equivalence in (5). Hence, the law of the SW, distance between a (random) MCMC history and the target distribution provides information on the integration error incurred by the MCMC kernel.
More formally, assume that the MCMC algorithm targets a distribution μı, and outputs random configuration of points (X1....,Xp) Xr). Call
T
MCMC He
8x.
T
the corresponding (random) empirical measure. Note that evaluating directly SW, (JMCMC 14) is not pos- sible, but if it is possible to sample Y,, i.idd from μ (as is often the case when testing sampling algorithms on simple targets), we can use
M
5%
M
as proxy for p. The triangular identity indeed guarantees
SW1 MCMC .4) < SW1 MCMC ( + SWi( n).
(25)
The second term of the right-hand side can be controlled via results involving the sample complexity (Manole et al., 2022), and scales as 1/VM. We thus focus on estimating SW, (WNNCMC (id), using Monte Carlo integration over the sphere. Note that, since our goal is to illustrate various quadratures on the sphere, we will consider a single realization of μd MCMC per dimension, but an MCMC practitioner wanting to estimate the quality of an MCMC kernel should repeatedly sample independent MCMC histories and consider the distribution of the obtained SW distances.
18

---

Published in Transactions on Machine Learning Research (02/2026)
We consider d E (10,30). Our target distribution is the banana-shaped target that is classically used to demonstrate the ability of gradient-based based MCMC samplers, such as Hamiltonian Monte Carlo (Duane et al., 1987), to make long-range jumps and thus reduce the asymptotic variance of the corresponding MCMC estimators. Formally, the banana-shaped target is the distribution of the image of a Gaussian vector X ~ W(O, Ta) by the map f: IR + T& defined by [2)+1(X) = -22+1 and (X)+ = -222+2 + (x25+1 - 5)2 for j 2 O. We further fix the number of projections in the SW estimators to N = 103, the number of points on which the reference measure is supported to M 104.
To obtain realizations of ff wXCMC, we consider four MCMC kernels from the PyMC v5.23.0. library (Abril-Pla et al., 2023), namely four variants of Hamiltonian Monte Carlo (HMC; (Duane et al., 1987)). HMC has several hyperparameters, such as a stepsize and a mass matrix parameters, and PyMC offers different options to tune them. Our first kernel (henceforth referred to as regular HMC) is the default automatic tuning in PyMC. Our second kernel (broken HMC) corresponds to us manually blocking the adaptation of the mass matrix, and setting it the identity matrix. Our third kernel (regular NUTS) is the No-U-Turn adaptive HMC sampler of Hoffman et al. (2014), as implemented again in PyMC. Our fourth kernel (broken NUTS) is NUTS, but with us manually blocking the online adaptation of the stepsize parameter -which renders the analysis of the Markov chain difficult- and setting it to fixed value. The objective it to observe the practical relevance of the hyperparameter tuning mechanisms in HMC.
In d = 10, we consider five estimators of the SWi distance, namely iidd, UnifOrtho, Repelled iidd, SHCV, and Repelled SHCV. For SHCV, we were able to set the maximum degree of spherical harmonics to 4. d = 30, we keep ii.d, UnifOrtho, and Repelled iidd. Note that we do not consider CV up and CV low, since they were specifically designed for SW2.
Our results for d = 10, 30 are respectively shown in Figures 7 and 8. We show the average of 1, 000 inde- pendent realizations of each estimator, with the 95% Gaussian confidence interval for the mean, Bonferroni corrected across the 5 x 4 estimators (respectively 3 x 4) corresponding to each value T. Strictly speaking, one can thus statistically compare all confidence intervals for any given value of T, but we should refrain from comparing across values of T. Since our objective is to compare the accuracies of various quadratures, this seemed a natural correction.
T
ii.d.
Repelled
UnifOrtho
SHCV
Repelled SHCV
4.955 + 0.093
4.957 + 0.062
4.958 + 0.028
4,960 = 0.049
4.977 = 0.037
4.959 + 0.101
4.960 + 0.068
4.958 + 0.026
4.937 = 0.052
4.957 0.040
10
4.957 + 0.089
4.957 + 0.073
4.958 + 0.026
4.956 0.060
4.959 0.039
4.805 + 0.102
4.811 + 0.062
4.811 = 0.023
4.812 0.050
4.812 0.041
0.741 + 0.007
0.741 + 0.005
0.741 + 0.002
0.741 + 0.004
0.741 0.003
0.847 + 0.012
0.848 + 0.008
0.848 + 0.003
0.847 + 0.007
0.848 0.006
100
4.051 1: 0.082
4.050 1: 0.054
4.051 = 0.020
4.049 = 0.044
4.051 = 0.029
0.579 1: 0.008
0.579 + 0.006
0.580 = 0.002
0.580 - 0.005
0.379 = 0.004
0,270 + 0.003
0.270 + 0.002
0.270 = 0.001
0.270 = 0.003
0.270 = 0.002
0.487 + 0.009
0.487 + 0.006
0.487 + 0.002
0.487 = 0.005
0.487 0.004
000
0.374 + 0.004
0.374 + 0.003
0.374 + 0.002
0.374 + 0.005
0.374 0.004
0.242 +: 0.004
0.242 + 0.003
0.242 1 0.001
0.242 = 0.003
0.242 0.002
0.158 + 0.002
0.1580 0.001
0.1580 + 8 10
0.158 0.002
0.1580 0.001
0.097 1 0.001
0.0968 7
0.0968 + -0
0.097 = 0.001
0.0968
10000
0.790 1: 0.016
0.790 1: 0.010
0.790 = 0.004
0.790 - 0.008
0.790 = 0.006
0.008 1 0.001
0.0984 10
0.0983 + 1000-
0,0953 10
0.0953 10-4
Figure 7: Averaged SW, and asymptotic confidence intervals in d = 10. The color code is Blue for broken HMC, Red for regular HMC, Green for broken NUTS, and Purple for regular NUTS.
The first conclusion is that UnifOrtho consistently yields smaller confidence intervals than the other methods, in both dimensions, which is why we display the corresponding column in bold. Similarly, the second conclusion is that repelled versions of each algorithm reduce the size of the confidence intervals in 10, but the improvement is less perceptible in d = 30. This is to be taken with a pinch of salt, however, as we do not provide confidence interval on the variance of the estimator.
19

---

Published in Transactions on Machine Learning Research (02/2026)
T
Li.d.
Repelled
UnitOrtho
2.135 + 0.017
2.135 + 0.017
2.135 = 0.002
1.919 + 0.010
1.918 = 0.010
1.918 = 0.001
10
4.057 +: 0.067
4.063 = 0.078
4.063 = 0.018
4.717 + 0.089
4.718 = 0.086
4.719 = 0.022
0.726 + 0.009
0.725 = 0.009
0.725 = 0.004
0.849 + 0.013
0.848 + 0.011
0.848 = 0.004
100
4.242 +: 0.077
4.251 + 0.079
4.247 = 0.020
0.529 + 0.007
0.529 + 0.008
0.529 = 0.003
0.363 + 0.005
0.363 + 0.004
0.363 = 0.001
0,288 + 0.003
0.288 + 0.004
0.288 = 0.001
1000
0.236 + 0.003
0.236 = 0.003
0.236 = 0.001
0.215 + 0.003
0.215 + 0.003
0.215 = 0.001
0.109 + 0.002
0.169 + 0.002
0.1693 6-
0.134 + 0.002
0.134 + 0.002
0.1341 + 10
10000
0.0645 8- 10
0.0645 10
0.0646 + 10
0.0638 7
0.0638
0.0638 3-
Figure 8: Averaged SW1 and asymptotic confidence intervals in d = 30. The color code is Blue for broken HMC, Red for regular HMC, Green for broken NUTS, and Purple for regular NUTS.
As to our mock goal to compare algorithms, Figure/ 7 shows first, for instance, no statistical gain in using NUTS rather than regular HMC when T = 10, 000. In d = 30, a similar phenomenon can be observed in Figure 8 when comparing broken NUTS and regular NUTS when T = 10, 000. In that case, only UnifOrtho has a small enough variance to yield a statistically significant difference in performance, in favor of regular NUTS, as expected. This time, the pinch of salt comes from our use of a single MCMC run for each pair of values of T and d. Still, as quadrature algorithms are concerned, UnifOrtho is to be preferred.
7 Discussion
Our empirical findings suggest that, when working in small dimensions (d E (2,3)), the lowest variance is obtained by randomizing simple deterministic quadratures. Indeed, the randomized spiral points in d = 3, and the classical grid in d outperform most sophisticated random methods, at cheap computational cost. When the dimension grows, these methods become unavailable, and more inherently random quadra- tures become attractive. Crude Monte Carlo, using i..d. uniform samples, quickly gets outperformed by most of the presented methods. Among them, DPPs are competitive in smaller dimensions, but their sam- pling cost becomes prohibitive as dimension increases. This is especially true for the harmonic ensemble, whose cardinality is bound to be exponential in the dimension, while sampling intermediate levels requires extensive calls the spherical harmonics and rejection sampling phase with loose rejection bound. On the other hand, the Repelled processes are cheap alternatives to DPPs that lead to (small) variance reduc- tion. Yet, their behavior is less well understood. While Appendix A.5 suggests that some of the intuition on tuning the repulsion carries over from the Euclidean case, combining the repulsion operator with e.g. control variates leads to unstable behavior. Turning to control variates methods, they consistently lead to variance reduction in our benchmark, yet they come with restrictions: CV up and CV low are limited to SW2, while SHCV requires the computation of spherical harmonics. Finally, a clear cost-efficient algorithm outperforms every other method in higher dimension: UnifOrtho. This is interesting, mator, yet with limited negative dependence due to the fact that it is the union of many independent small repulsive point processes. We contributed to the understanding of the success of UnifOrtho by providing an expression for the variance the corresponding estimator in terms of the spherical harmonics coefficients of the integrand, which also explains why variance can actually increase if applied to integrands with specific spectral profiles. As stated in our introduction, our observations are independent of the ML task that will use the SW distance, and the important factors affecting the performances of the methods are the dimension of the problem and the number of points on which the measures are supported. Another avenue for future work is to combine UnifOrtho and control variates provide a uniform decrease in variance. Similarly,
it is repulsive Monte Carlo esti-
20

---

Published in Transactions on Machine Learning Research (02/2026)
understanding the spectral profile of the integrand in the SW distance, in terms of easy-to-estimate features of the two involved distributions, would help choosing the right estimator.
Acknowledgments
We acknowledge support from ERC grant BLACKJACK ERC-2019-STT--5D7 8551666 and ANR Al chair BACCARAT ANR-20-CHIA-0002.
References
O. Abril-Pla, V. Andreani,, C. Carroll, L. Dong, C. J. Fonnesbeck, M. Kochurov, R. Kumar, J. Lao, C. C. Luhmann, and O. A. Martin. PyMC: a modern, and comprehensive probabilistic programming framework in Python. PeerJ Computer Science, 2023.
G. W. Anderson, A. Guionnet, and O. Zeitouni. An introduction to random matrices, volume 118. Cambridge university press, 2010.
R. Bardenet and A. Hardy, Monte Carlo with determinantal point processes. Annals of Applied Probability, 2020.
S. Barthelmé, N. Tremblay, and P.-O. Amblard. A faster sampler for discrete determinantal point processes. International Conference on Artificial Intelligence and Statistics, pp. 5582-5592. PMLR, 2023.
E. Bayraktar and G. Guo. Strong equivalence between metrics of Wasserstein type. Electronic Communications in Probability, pp. 1-13, January 2021.
A. Belhadji. An analysis of Ermakov-Zolotukhin quadrature using kernels. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
A. Belhadji, R. Bardenet, and P. Chainais. Kernel quadrature with determinantal point processes. In Advances in Neural Information Processing Systems (NeurIPS), 2019.
A. Belhadji, R. Bardenet, and P. Chainais. Kernel interpolation with continuous volume sampling. In Inter- national Conference on Machine Learning (ICML), 2020.
C. Beltrán and U. Etayo. A generalization of the spherical ensemble to even-dimensional spheres. Journal of Mathematical Analysis and Applications, 475:1073- July 2019.
C. Beltrán, J. Marzo, and J. Ortega-Cerdà Energy and discrepancy of rotationally invariant determinantal point processes in high dimensional spheres. Journal of Complexity, pp. 76-109, December 2016.
R. J. Berman. The spherical ensemble and quasi-Monte-Carlo designs. Constructive Approximation, 59: 457-483, April 2024.
C. Bonet, N. Courty, F. Septier, and L. Drumetz. Efficient Gradient Flows in Sliced-Wasserstein Space. Trans- actions on Machine Learning Research, 2022.
Clément Bonet, Lucas Drumetz, and Nicolas Courty. wnnsersssse distances and flows on cartan- hadamard manifolds. Journal of Machine Learning Research, 26(32):1-76, 2025.
N. Bonneel, J. Rabin, G. Peyré, and H. Pfister. Sliced and Radon Wasserstein Barycenters of Measures. Journal of Mathematical Imaging and Vision, pp. 22-45, January 2015.
N. Bonnotte. Unidimensional and Evolution Methods for Optimal Transportation. Phd thesis, Université Paris Sud Paris XI; Scuola normale superiore (Pise, Italie), December 2013.
J. Brauchart, E. Saff, I. Sloan, and R. Womersley. QMC designs: optimal order quasi Monte Carlo integration schemes the sphere. Mathematics of computation, 838290112811-2811, 2014.
21

---

Published in Transactions on Machine Learning Research (02/2026)
G. V. Cardoso, Y. J. El Idrissi, S. Le Corff, and E. Moulines. Monte Carlo guided diffusion for bayesian linear inverse problems, 2023.
A.X. Chang, T.A. Funkhouser, L.J. Guibas, P. Hanrahan, Q. Huang, Z. Li, S. Savarese, M. Savva, S. Song, H. Su, J. Xiao, L. Yi, and F. Yu. ShapeNet: An Information-Rie -iich 3D Model Repository. CoRR, abs/1512.03012, 2015.
J. Cheeger. Differentiability of Lipschitz Functions Metric Measure Spaces. Geometric & Functional Analysis GAFA, 9:428-517, June 1999.
J.-F. Coeurjolly, A. Mazoyer, and P.-O. Amblard. Monte Carlo integration of non-differentiable functions on [0, 1J', 1, ...... using single determinantal point pattern defined on [0, 1]a. Electronic Journal of Statistics, 15()):6228-624) 2021.
F. Dai and Y. Xu. Spherical harmonics. In Approximation theory and harmonic analysis on spheres and balls, pp. 1-27. Springer, 2013.
B. Delyon and Portier. Integral approximation by kernel smoothing. Bernoulli, 2016.
1. Deshpande, Z. Zhang, and A. Schwing. Generative Modeling Using the Sliced Wasserstein Distance. 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3483-3491, June 2018.
J. Dick and Pilichshammer. Digital Nets and Sequences. Discrepancy Theory and Quasi-Monte Carlo Integra- tion. Cambridge University Press, 2010.
S. Duane, A. D Kennedy, B. J Pendleton, and D. Roweth. Hybrid Monte Carlo. Physics letters B, 195(2): 216-222, 1987.
V. Dutordoir, N. Durrande, and J. Hensman. Sparse Gaussian processes with spherical harmonic features. In International Conference on Machine Learning, pp. 2793-2802. PMLR, 2020.
G. Fort, B. Jourdain, E. Kuhn, T. Lelièvre, and G. Stoltz. Convergence of the wang-landau algorithm. Mathe- matics of Computation, 84(295):2297-33377 2015.
N. Fournier and A. Guillin. On the rate of convergence in wasserstein distance of the empirical measure. Probability theory and related fields, 162(3):707-738, 2015.
G. Gautier. On sampling determinantal point processes. PhD thesis, Centrale Lille Institut, 2020.
G. Gautier, R. Bardenet, and M. Valko. On two ways to use determinantal point processes for Monte Carlo integration. In Advances in Neural Information Processing Systems (NeurIPS), 2019a.
G. Gautier, G. Polito, R. Bardenet, and M. Valko. DPPy: DPP Sampling with Python. Journal of Machine Learning Research, pp. 1-7, 2019b.
W. Gautschi. Orthogonal polynomials: computation and approximation. Oxford University Press, USA, 2004.
D. Hawat, R. Bardenet, and R. Lachièze Rey. Repelled point processes with application to numerical integra- tion. In revision for Scandinavian Journal of Statistics, 2023.
Eric Heitz, Kenneth Vanhoey, Thomas Chambon, and Laurent Belcour. A sliced wasserstein loss for neural texture synthesis. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9412-9420, 2021.
M. D Hoffman, A. Gelman, et al. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res., 15(1):1593-1623, 2014.
J. Hough, M. Krishnapur, Y. Peres, and B. Virág. Determinantal processes and independence. Probability surveys, 2006.
S. Kolouri, Y. Zou, and G. K. Rohde. Sliced Wasserstein kernels for probability distributions. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 5258-5267 2016.
22

---

Published in Transactions on Machine Learning Research (02/2026)
Soheil Kolouri, Phillip E. Pope, Charles E. Martin, and Gustavo K. Rohde. Slisee-wwssersteinn auto-encoders. In International Conference on Learning Representations (ICLR), 2019.
M. Krishnapur. From random matrices random analytic functions. The Annals of Probability, January 2009.
D. P. Kroese, R. Y. Rubinstein, and P. W. Glynn. Chapter 2 The Cross-Entropy Method for Estimation. In C. R. Rao and Venu Govindaraju (eds.), Handbook of Statistics, Handbook of Statistics, pp. 19-34. Elsevier, January 2013.
A. Kulesza and B. Taskar. Determinantal point processes for machine learning, Foundations and Trends in Machine Learning, 2012.
F. Lavancier, J. Møller, and E. Rubak. Determinantal point process models and statistical inference. Journal of the Royal Statistical Society, 2015.
R. Leluc, A. Dieuleveut, F. Portier, J. Segers, and A. Zhuman. Sliced-Wasserstein Estimation with Spherical Harmonics as Control Variates. International Conference on Machine Learning, February 2024.
R. Leluc, F. Portier, J. Segers, and A. Zhuman. Speeding up Monte Carlo integration: Control neighbors for optimal convergence, 2025.
T. Lemoine and R. Bardenet. Monte Carlo methods on compact complex manifolds using Bergman kernels. arXiv preprint arXiv:2405.09203, 2024.
H. Lin, H. Chen, K. M Choromanski, T. Zhang, and C. Laroche. Demystifying Orthogonal Monte Carlo and Beyond. In Advances in Neural Information Processing Systems (NeurIPS), pp. 8030-8041. Curran Associates, Inc., 2020.
J. Linhart, G. V. Cardoso, A. Gramfort, S. Le Le Corff, and P. LCC Rodrigues. Diffusion posterior sampling for simulation-based inference in tall data settings, 2024.
A. Liutkus, U. Simsekli, S. Majewski, A. Durmus, and F.-R. Stöter. -asserseass flows: Nonparametric generative modeling via optimal transport and diffusions. In International Conference on machine learning, pp. 4104-4113. PMLR, 2019.
O. Macchi. Processus ponctuels et coincidences Contributions à l'étude théorique des processus ponctuels, avec applications à l'optique statistique et aux communications optiques. PhD thesis, Université Paris-Sud, 1972.
T. Manole, S. Balakrishnan, and L. Wasserman. Minimax confidence intervals for the sliced Wasserstein distance. Electronic of Statistics, 16(1):2252-2345, 2022.
J. Marzo Sánchez, M. Levi, and Ortega Cerdà. Linear statistics of determinantal point processes and norm representations. International Mathematics Research Notices, 2024, vol. 2024, num. 19, p. 12869-12903, 2024.
A. Mazoyer, J.-F. Coeurjolly, and P.-O. Amblard. Projections of determinantal point processes. Spatial Statis- tics, 38:100437, 2020.
E. S. Meckes. The Random Matrix Theory of the Classical Compact Groups. Cambridge Tracts in Mathematics. Cambridge University Press, 2019.
K. Nadjahi. Sicedd-Wasersteein distance for large-scale machine learning theory, methodology and extensions. Phd thesis, Institut polytechnique de Paris, November 2021.
K. Nguyen and N. Ho. Sliced Wasserstein Estimation with Control Variates. International Conference on Learning Representations, February 2024.
K. Nguyen, N. Bariletto, and N. Ho. Quasi-Monte Carlo for 3D Sliced Wasserstein. International Conference Learning Representations, February 2024.
23

---

Published in Transactions on Machine Learning Research (02/2026)
Art B. Owen. Monte Carlo theory, methods and examples. 2013.
G. Peyré and M. Cuturi. Computational Optimal Transport. Poundations and Trends in Machine Learning, 11 (5-6):355-2000 2018.
F. Portier and J. Segers. Monte Carlo integration with growing number of control variates. Journal of Applied Probability, 56(4):1168-1188) 2019.
E. A. Rakhmanov, E. B. Saff, and Y. M. Zhou. Minimal Discrete Energy on the Sphere. Mathematical Research Letters, pp. 647-662, 1994.
Alain Rakotomamoniy and Ralaivola Liva. Differentially private sliced wasserstein distance. In International Conference on Machine Learning, pp. 8810-8820. PMLR, 2021.
B. Rider and B. Virag. Complex Determinantal Processes and H Noise. Ciccronni of Probability, pp. 1238-1257, January 2007.
C. P. Robert and G. Casella. Monte Carlo statistical methods. Springer, 2004.
M. Rowland, J. Hron, Y. Tang, K. Choromanski, T. Sarlos, and A. Weller. Orthogonal Estimation of Wasser- stein Distances. In Proceedings of the Twenty-Second International Conference Artificial Intelligence and Statistics, pp. 186-195. PMLR, April 2019.
B. Rubin. On the injectivity of the shifted Funk-Radon transform and related harmonic analysis. Journal d'Analyse Mathématique, 153:777-800, September 2024.
K. Sisouk, J. Delon, and J. Tierny. A user's guide to sampling strategies for sliced optimal transport, 2025.
S. Smale. Mathematical problems for the century. The Mathematical Intelligencer, 20:7-15, March 1998.
A. Soshnikov. Determinantal random point fields. Russian Mathematical Surveys, 55:923-975, 2000.
S. Sra. A short note on parameter approximation for von Mises Fisher distributions: and fast implementa- tion of Is(x). Computational Statistics, pp. 177-190, March 2012.
Mao Ye, Tongzheng Ren, and Qiang Liu. Stein self-repulsive dynamics: Benefits from past samples. Advances in Neural Information Processing Systems, 33:241-252, 2020.
A Appendix
A.1 Spherical harmonics
The spherical harmonics are class of functions that play an important role when approximating functions on the sphere. We refer to Dai & Xu (2013) for a comprehensive introduction, from which we isolate a few points here for completeness.
Let d 2, the simplest definition that the spherical harmonies on $d-1 are the homogeneous harmonic polynomials of Rd restricted to the sphere gd-1, Alternately, if A is the Laplace- Belttrai operator gd-1 and X; f(E+de d d-2), 2), the spherical harmonics of order & N can be defined as the elements of the eigenspace He of A corresponding to eigenvalue - One can then show that He is the set of harmonic homogeneous polynomials of degree restricted to gd-1 as expected. Furthermore,
m
θ He
is the space of harmonic polynomials in IRd restricted to Sd-1 of degree up to L. We also note, following Marzo Sánchez et al. (2024), that
T == dim 2L + (d - 1) (d (d- ()
2 - + o(T.d-)).
(26)
d-1
L
[(d)
24

---

Published in Transactions on Machine Learning Research (02/2026)
For a given f, let h2 dim(H)) and { I h I he) be any orthonormal basis of He. Then Dai & Xu (2013)
[Theorem 1.2] state that the elements of (Ygle E N, < k S he) are centered functions for the uniform
measure on the sphere, as soon as / 2 1, which form a Hilbert basis of L3)(s-1).
From a computational standpoint, it is often useful to note the following addition formula Dai Xu (2013)
[Theorem 2.6],
hy ((
X (a,,)))
vz, y E 1 Zz(x.)) =
(27)
where O2 is the Gegenbauer polynomial of degree l and x = d-2 This leads to the following definition.
Definition A set of points f- Thee gd-1 is said to be fundamental if the matrix C₁
(G)( is invertible.
Fundamental sets are particularly interesting since, (z1,-.,,,,, is fundamental set, then (C)((, S < he] is basis of He (Dai Xu, 2013) [Theorem 3.3]. This theorem is at the heart of the library' developed by Dutordoir al. (2020) to compute spherical harmonics. Their method consists in greedily building fundamental set that is likely to lead to stable Cholesky decomposition C₁. This is done by iteratively adding point that maximizes the determinant of (C2(, iskjsh Then, through a Cholesky decomposition of C₁ they obtain the (C)((, si<hi ht). In other words, they obtain an orthonormal basis of Hi. orthonormalization of
The greedy construction of a fundamental set is computationally heavy, although one has to only run once only. As computationally cheaper alternative and at the price of stability of the Choleky decomposition, the point sets which are not fundamental lie in {det(C,) so is of measure zero. Hence almost every set of points is a fundamental set (Dai & Xu, 2013). However, cost that cannot be avoided is that all the spherical harmonics of a given level have to be computed, which comes down to finding the Cholesky decomposition of an fig x h& matrix, and he grows as gd-2
0}, which is an algebraic hypersurface of g(d-1)xh
A.2 More on the importance sampling scheme
For completeness, and because fitting a von-Mises-Fishee distribution is not straightforward, we provide here pseudocode in Algorithm for our fitted importance sampling estimator.
A.3 Pseudocode for competing methods
For completeness, add pseudocode to describe the methods compare to. Algorithm 2 describes the classical HKPV algorithm to sample from DPP (see Section 4.2), Algorithm 3 describes the SHCV method (see Section 3.2.1), Algorithm 4 describes the crude Monte Carlo with a single control variate algorithm (see Section 3.2.1), Algorithm 5 describes the UnifOrtho sampling strategy (see Section 3.te that, the SHCV implementation described here highlights the fact that the weights of the quadrature are independent of the function to integrate.
A.4 Discussion the shape of the integrand
We include in Figure various plots of the integrand (6) of the sliced Wasserstein distance in three di- mensions. In Figure 9a, we show the integrand corresponding to two Gaussians with random means and covariances, specified in Section 6.1. In Figure 9b, we examine the integrand (6), but this time between two empirical measures based on respective i.idd draws from the same two Gaussians. The integrands in Figures 9a and 9b are visually similar, as expected. Moreover, they seem to be unimodal up to symmetry. In particular, we expect importance sampling with a fitted symmetrized vMF proposal to yield low variance.
Figures 9c and 9d show the integrand (6) for the point clouds used in Section 6.2. Here the landscape seems more erratic, and the regularity well the number of modes is less straightforward determine. Yet,
25

---

Published in Transactions on Machine Learning Research (02/2026)
Algorithm 1 A cross-entropy-fitee importance sampling estimator.
1: Input: Measures and v, Number N of points to sampled, Budget fraction r E (0, 1) to allocate to estimating the proposal.
2: Sample , .._mJ ii.d. from the uniform measure on the sphere. Evaluate the integrand fup in 6 on them. Define
i mmx = (fg & [ra]].
Define ) = [111( mmax x) > 0], and evaluate the quantities
IrN]
Σ (X.)x
Σ P. (x
[|NN]
BINNI = L
(x.),.l
Σ (x.)
Σ
[irN] (d R. [
4: Let [[rN]
as in Sra (2012).
R² [nN]
5: Sample X [rN]+I..... X from ( ((NN) + vmf(-E ([(rN]))
Return
|N Σ
1-r
N
(X:) ( (((r)) + vmf(X,] - ELNNJ K(NN))
(x.) +2
N
NM
IINN +1
Algorithm 2 The HKPV algorithm for sampling the harmonic ensemble of size N Hough et al. (2006).
Input: N, the kernel K.
Sample 0g d0, on the sphere for 2 < k N do
Use rejection sampling to draw
K(0,, 0k) 10 01:--1 1) 1 (01: dAk
0x
-
where K ((0,, 11k-11 K(O,, 01 (1)..,0)....,.... 1T) and K;-1 = (K(,,, (iok,,;)))ssh 1j)
end for 8: Return (fI....,0r 1N}
Algorithm 3 The SHCV method Lelue et al. (2024)
Input: N, spherical harmonics (fi))=
Sample f1,... ON i.i.d. uniformly on the sphere.
Compute the matrices & (p1( H TT'
Compute the vector fy (Sa. (0.))))s)
5: Let ll = H HH H" and 1 = (I, 1I,
Return SHCV (p, (T)-
(1N-
intuitively the resulting integrands should have a relatively sparse decomposition in the bases of spherical harmonics, confirmed by Figure 6.
26

---

Published in Transactions on Machine Learning Research (02/2026)
Algorithm 4 Crude Monte Carlo integration with single control variate Nguyen & Ho (2024)
1: Input: N, single control variate vp Sd-1 + IR
2: Sample .1,...,0N NN i.i.d. uniformly on the sphere.
N
3: Compute JN = (Fun(O:)) UN
Σ (0:)
N
4: Compute 3
(0)))())-
5: Return N (((()) Bx((.))
Algorithm 5 Sampling from the Haar measure (Meckes, 2019, Paragraph 1.2)
Sample XI,... Xa i.i.d. (O, Compute the QR decomposition of (X1l... Xa) via the Gram-Schmidt algorithm.
3: Return the matrix Q.
A.5 A few words on repelled point processes
The repelled point processes behave inconsistently in the experiments of Section 6. Intuitively, a small repulsive perturbation, meaning a small > in (19), should lead to some variance reduction, yet the magnitude of seems depend on the number N of points to repel as well as on their distribution. this section, we experimentally investigate this optimal choice for 6, guide future theoretical investigations. Following the seminal case of a homogeneous Poisson process in Rd expect variance reduction to happen when is of the order of /N, where N is the cardinality of the configuration to repel.
To assess the influence of e, we sampled 100 independent realizations of the repelled estimator for each of discrete set of values of e and a choice integrands and initial point processes, and reported the x? confidence interval for the variance of each estimator. We correct these confidence intervals with Bonferroni correction across the finite set of values for c, reach simultaneous confidence level of 0.969.
Figure 10 shows the results for an initial point process made of i.idd uniform points on the sphere, for choice of indicators and values of N. The red line is reference corresponding to 0, placed at an arbitrary low value of for comparison. The gray dashed line corresponds to c 1/N. We observe that 1/N is indeed a sensible choice, leading to variance reduction by a factor up to 2, even for a non-smooth integrand like the indicator of a half-sphere.
Using the i.i.d. repelled points build the SHCV estimator, as denoted by Repelled SHCV in Section 6, the situation becomes much more unstable, shown in Figure 11. Looking closely, it seems that choice of slightly under 1 1/N seems consistently lead to some variance reduction, but small variation in can have drastic consequences, observed on Figure 11b. Note also that, unlike the vanilla repelled i.i.d. estimator, the optimal choice for seems to be dimensioon-eependen,, as shown by the difference in the dips between Figure 10b and Figure 11e. Hence, although repelling the points in the use of SHCV has the potential to diminish the variance, further theoretical investigation are required to correctly tune e.
Finally, in line with the experimental observations of Hawat et al. (2023), repelling structured points can also lead to straight-up increase in the variance. This is the case when the starting configuration is randomized grid, as in the QMC or the UnifOrtho estimators in Section 6, see Figure 12. Yet, Figure 12b suggests that Repelled UnifOrtho can actually lead to dip in variance when integrating an indicator. This was an unexpected behavior and further theoretical work is needed to understand this phenomenon.
27

---

Published in Transactions on Machine Learning Research (02/2026)
-8
E.M
in (a#t)
(atto
2
2
(a) Heatmap of the true integrand between two Gaus- sian measures
(b) Heatmap of the integrand between two point clouds sampled from these Gaussian measures
0.05
0.12
0.04
-..10
00.08
0.03
-0.06
0.02
-0.04
-0.01
0.02
(c) Heatmap of the integrand between point clouds #2 and #34
(d) Heatmap of the integrand between point clouds #3 and #35
Figure 9: Heatmaps of the integrand in 3D for various distributions.
28

---

Published in Transactions on Machine Learning Research (02/2026)
0055
DUISS
10000
0015
2014
2002
0.000
0.030
10-6
10mm
10
10-9
10
(a) Integration of the indicator of half-sphere with N = 100 li.d. points, d =
Integration of the indicator half-sphere with N = 1000 i..d. points, d = 3
0.45
0.35
0.22
0.10
0.25
0.08
0.00
1015
-0"1
107
(c) Integration the sliced Wassertein between sampled Gaussian supported 100 points with N 100 i.idd points, 3
(d) Integration the sliced Wassertein between sampled Gaussian supported on 100 points with N 1000 i.i.. points, d 3
0.22
0.07-
0.70
0.06
0.05
0.01
0.14
10
10" 20
10 19
207*
13"*
10
10
(e) Integration sliced Wassertein between sampled Gaussian supported on 100 points with N 100 i.idd points, 10
(f) Integration of the sliced Wassertein between sampled Gaussian supported on 100 points with N 1000 i.id. points, d 10
Figure 10: Confidence intervals for the variance of the estimator with a Bonferroni corrected confidence level of 0.969, where the repelled points are sampled i.i.d.
29

---

Published in Transactions on Machine Learning Research (02/2026)
0030
0078
0070
0,06-
0.02
0.020
0010
10-6
10
10-9
(a) Integration of the indicator of halfrsphere with N = 100 and SHCV, d = 3
Integration of the indicator half-sphere with N = 1000 and SHCV, d =
20007-
1022
00244
20006
1022
30025
0016-
0016
anic
0012
0.005
0031
10
-0"
13
107
(c) Integration the sliced Wassertein between sampled Gaussians supported 100 points with = 100 and SHCV, = 3
(d) Integration the sliced Wassertein between sampled Gaussians supported 100 points with N = 1000 and SHCV, = 3
0.045
0.080
0.001
0.005
10
10
20
(e) Integration of the sliced Wassertein between sampled Gaussians supported on 100 points
with N 1000 and SHCV, d =
Figure 11: Confidence intervals for the variance of the estimator with a Bonferroni corrected confidence level of 0.969 where the SHCV estimator is built on repelled i.i.d. uniform points.
30

---

Published in Transactions on Machine Learning Research (02/2026)
0.022
0.001
0.000
0.015
0.005-
00006
000024
(a) Integration of the indicator of half-sphere with N = 1000 RQMC points, d = 3
0.12
0008
0.10
1006
0.05
1
0.00
I0022
0.02
0008
70.2
10"
10"
1/P
(c) Integration of the sliced Wassertein between sampled Gaussian supported on 100 points with N 1000 UnifOrtho points, 3
(b) Integration of the indicator of half-sphere with 1000 UnifOrtho points, 3
Confidence intervals for the variance of the estimator with a Bonferroni corrected confidence level of 0.969, when repelling points sampled using either UnifOrtho or RQMC d = 3.
Figure 12:
31


Reviews:
[{"review_id": "TOeKgHBFlK", "reviewer": "Reviewer_XjpC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors evaluate the Monte-Carlo method for computing the sliced Wasserstein distance (SWD), which can be used as a proxy to Wasserstein distance that is widely used in OT but always intractable especially in the high-dimensional setting. The contributions include benchmark randomized quadratures for estimating SWD; after identifying the empirically most successful MC based estimator, the authors further provide theoretical understanding of of the concentration in terms of the dimension $d$ and the basis $Y^\\ell_k$.", "claims_and_evidence_yn": "Yes", "claims_explanation": "I have the following concerns and questions;\n\n1. Can you explain more on the notations in Proposition and make it self-contained? For instance, the definition of $l, h_l,k$ in describing the basis.\n2. My main concern is on Proposition 4, to me it is essentially a bias-variance tradeoff, where the proof starts with the variance-covariance-bias expression of $Var(\\frac{1}{d}\\sum^d_{i=1} f(X_i))$, in the proof, there are some steps not clear to me:\n\ni) How the author obtain the decomposition of $f(x) = \\sum_{\\ell = 0}^{\\infty} \\sum_{k = 1}^{h_\\ell} \\widehat{f}(\\ell,k)\\, Y_k^{\\ell}(x)$, any assumption is made during the decomposition?\n\nii) by integrating the decomposition of $f$ to (21), can the authors provide more details on the simplification? Because to me, it is not that obvious to get $\\mathbb{E}\\big[f(X_{1}) f(X_{2})\\big] = \\sum_{\\ell = 0}^{\\infty} (-1)^{\\ell} \\lambda_{2\\ell} \\, \\mu_{2\\ell}(f)$ directly.\n\n3. The comment on the Proposition 4 is intuitive and can also be supported by some simulation results. When the authors justify that 'When the dimension $d-2$ grows, the generalized Stirling formula yields  $\\lambda_{2\\ell} = \\mathcal{O}(\\ell^{-\\frac{d-2}{2}})$, so that only the first terms in either (19) or (20) carry significant weight.', if I understand correctly, when dimension increase, the variance of UnifOrtho estimator will be close that that of crude MC, can this also be supported in the simulation?", "audience_yn": "Yes", "audience_explanation": "The introduction part of the paper is very informative to the readers in the related area. The simulation part is also comprehensive in comparing the existing MC based approaches for computing SWD, which can be used as a guideline to pick the good strategy.", "requested_changes": "Each section and paragraph especially when previous approach in the literature is introduced, the notation as well as its explanation is supposed to be self-contained, the current version is a bit hard to follow since the notations are repeatedly used across sections, and many important results are being referred to the original paper. It would make the paper more readable if the authors can summarize and compare the existing MC based approach in terms of the variance reduction rate to give a more straightforward comparison and high-level understanding of these approaches.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cvVW6dIODp", "reviewer": "Reviewer_TXUM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes new empirical methods of estimating the Sliced Wasserstein distance and does a detailed comparative study of it.", "claims_and_evidence_yn": "No", "claims_explanation": "**strengths**\n\nThe experimental study is fairly detailed.\n\n**weakness**\n\nThe paper is quite badly written - particularly so, considering that its a submission to a ML journal and hence the baseline expectation would be to expect a clear connection to ML being made. \n\nApart from the general issue above, the specification of the different algorithmic techniques being presented is quite vague. Its hard to distinguish between which parts are general commentaries and which parts are specification of the algorithm to be tested.", "audience_yn": "Yes", "audience_explanation": "Although this paper is not my area specialization, from the introduction given I am fairly confident that this topic can be of interest to many people - if the paper is extensively revised to be a much better writing.", "requested_changes": "All my requests below are essential to be addressed for securing acceptance.\n\n- The second paragraph of page-3 gives some details on how computing the sliced-Wasserstein (SW) distance impacts other conventional aspects of ML. But this is too thin a description and motivation for readers like me, who are ML theory/optimization theory specialists, to understand the need and scope of this question. My first and foremost suggestion is that this needs to be expanded into a detailed section explaining the interface between usual ML questions and SW distance. The beginning of the paper needs to be a lot more ML-friendly. \n\n  What things in ML can improve if we have better estimators of SW distance?\n\n- The section 4 on new candidate estimators is quite unclear. For a start we need an explicit definition of the von-Mises-Fisher distribution. Can this be organized into separate pseudocodes for computing SW via each of the proposed new methods of sampling points on the sphere? \n\n- The situation with the UnifOrtho seems a bit unclear. The discussion in section 3.2.2 seems to suggest that UnifOrtho is not a good estimator. Then given all these discussions why do we have a separate Section 5 dedicated to discussing it? Isnt that discussion in Section 5 better suited as an appendix?\n\n- The beginning of Section 6 needs a lot of reorganizing. I think this bullet point list is very hard to parse, particularly because multiple methods are being referred to the same previous subsection. What we need are separate numbered equations and pseudocode for each of the methods being tested and those be referred to from this list to indicate the methods being tested. \n\n-  Lastly, for each of the experiments, we need to see a final comparison table of performance between the different techniques which makes it clear as to who won that race and on what metrics. And there needs to be a final reflection on what methods are being recommended by the authors for the readers. It is to be expected that the final recommendation would depend on the ML question of the end-user, but that dependency is hardly becoming clear.", "broader_impact_concerns": "N/A", "additional_comments": "N/A", "extra_scores": {}, "extra_text": {}}, {"review_id": "HwqGAsKriQ", "reviewer": "Reviewer_AZwZ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper investigates the use of repulsive Monte Carlo methods for the numerical integration required to compute the sliced Wasserstein distance. The core contribution is a unified comparison of various quadrature techniques, focusing on those that enforce \"repulsion\" between sampling points on a high-dimensional sphere. The authors introduce and benchmark methods based on DPPs and repelled point processes. Another contribution is a novel theoretical analysis of the variance of the UnifOrtho estimator, which leverages spherical harmonics to explain its strong empirical performance in high dimensions.\n\n---\n\n**Strengths:**\n\n- The paper is well-written and provides a clear guide for approximating the integrand.\n- The empirical study is thorough, comparing a wide and relevant range of methods across multiple dimensions and tasks.\n- The new variance analysis for the UnifOrtho estimator is a strength. It provides a deeper theoretical understanding of the variation reduction feature of the UnifOrtho algorithm.\n\n---\n\n**Weaknesses:**\n\n- Although the paper is framed around repulsive Monte Carlo, its survey and experimental comparisons are not exhaustive. The authors focus primarily on DPPs and a specific potential-based repulsion scheme. Other relevant algorithms from this domain are not discussed or benchmarked, such as methods based on Wang-Landau, Stein's method, or adaptive biasing potential method.\n- While the investigation is valuable, the newly benchmarked \"Repelled\" processes are shown to offer nearly negligible performance difference among them, e,g., in Figure 7 and Figure 8.\n- The DPP-based methods are demonstrated to be computationally prohibitive in higher dimensions, limiting their practical applicability for this task. Thus, the scope of this paper sounds a bit narrowed, ignoring the comparisons to other repulsive Monte Carlo methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "Proposition 4 provides a precise, closed-form expression for the variance of the UnifOrtho estimator. This formula directly links the variance to the spherical harmonic coefficients of the function being integrated.\n\nThe central, practical claim, i.e., practitioners should use randomized QMC in low dimensions and UnifOrtho in high dimensions, is backed by a comprehensive set of experiments in Section 6.", "audience_yn": "Yes", "audience_explanation": "The findings of this paper would be of significant interest to a substantial portion of the TMLR audience, particularly those who work with generative models, optimal transport, and high-dimensional probability distributions.\n\nThe paper directly addresses a practical, foundational problem: how to efficiently and accurately estimate the sliced Wasserstein, which is a popular loss function for training GANs, VAEs, and one-dimensional OT.", "requested_changes": "1. The paper frames its contributions around \"repulsive Monte Carlo\" but primarily focuses on DPPs and a specific potential-based repulsion scheme. The work would be strengthened by briefly acknowledging the broader literature on repulsive and adaptive sampling in the introduction or related work section (Section 2). Mentioning other prominent approaches, such as [1,2,3,4], and briefly justifying the paper's focus on DPPs would make the survey feel more complete and avoid questions about why certain methods were excluded from the benchmark. This does not require new experiments, but rather a few sentences to better posit the work.\n\n> [1]. Mao Ye, Tongzheng Ren, and Qiang Liu. \"Stein self-repulsive dynamics: Benefits from past samples.\" Advances in Neural Information Processing Systems 33 (2020): 241-252.\n>\n> [2]. Gersende Fort, et al. \"Convergence of the Wang-Landau algorithm.\" Mathematics of Computation 84.295 (2015): 2297-2327.\n>\n> [3]. Michel Benaïm, and Charles-Édouard Bréhier. \"Convergence of adaptive biasing potential methods for diffusions.\" Comptes Rendus. Mathématique 354.8 (2016): 842-846.\n>\n> [4]. Gersende Fort, et al. \"Convergence and efficiency of adaptive importance sampling techniques with partial biasing.\" Journal of Statistical Physics 171.2 (2018): 220-268.\n\n2. The finding that the simple repulsive method provides only moderate variance reduction is an interesting result. The paper would benefit from a short discussion in the conclusion speculating on why this might be the case. Is it a fundamental limitation of this specific energy function, or could more sophisticated repulsion schemes potentially yield better results? This would add valuable insight for future work in this direction.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)