A User's Guide to Sampling Strategies for Sliced Optimal Transport
Authors: Keanu Sisouk
Julie Delon
Julien Tierny
Venue: Accepted by TMLR ()

Abstract: This paper serves as a user's guide to sampling strategies for sliced optimal transport. 
We provide reminders and additional regularity results on the Sliced Wasserstein distance.
We detail the construction methods, generation time complexity, theoretical guarantees, and conditions for each strategy. Additionally, we provide insights into their suitability for sliced optimal transport in theory. Extensive experiments on both simulated and real-world data offer a representative comparison of the strategies, culminating in practical recommendations for their best usage.

Published in Transactions on Machine Learning Research (06/2025)
A User's Guide to Sampling Strategies for Sliced Optimal Transport
Keanu SISOUK
Keanu.Sisouk@liip..f
CNRS, LIP6
Sorbonne University
Julie Delon
MAP5
Paris Cite University
Julien Tierny
julien. berruyggBorboonnn--iievrsss./.
CNRS, LIP6
Sorbonne University
Reviewed on OpenReview: https: // openreview. net/ forum? id- CC
Abstract
This paper serves as a user's guide to sampling strategies for sliced optimal (ransport [Rabin et al., 2012; Bonneel et al., 2015]. We provide reminders and additional regularity results on the Sliced Wasscrstein distance. We detail the construction methods, generation time com- plexity, theoretical guarantees, and conditions for each strategy. Additionally, we provide insights into their suitability for sliced optimal transport in theory. Extensive experiments on both simulated and real-world data offer representative comparison of the strategies, culminating in practical recommendations for their best usage.
1 Introduction
The Wasserstein distance is acclaimed for its geometric relevance in comparing probability distributions, Having gathered lot of theoretical work [Santambrogio, 2015; Villani, 2008], it has also proved to be relevant in numerous applied domains the last fiftcon years, such as image comparison [Rabin et al., 2009], image registration [Feydy et al., 2017], domain adaplation [Courly al., 2016], generative modeling [Arjovsky al., 2017; Guirajani et al., 2017; Salimans et al., 2018], inverse problems in imaging [Hertrich et al., 2022] topological data analysis [Edelsbrunner & Haror, 2009; Pont al., 2021], to name just few. The computational demands of the Wasserstein distance are, however, quite high, since evaluating the distance between two discrete distributions of N samples with traditional linear programming methods incurs a runtime complexity of O(N3 log N) [Peyré et al., 2019]. This computational burden has motivated the development of alternative metrics sharing some of the desirable properties of the Wasserstein distance but with reduced complexity.
The Sliced Wasserstein (SW) distance [Rabin et al., 2012; Bonneel et al., 2015], defined by slicing the Wasser- stein distance along all possible directions the hypersphere, is one of these efficient alternatives. Indeed, the SW distance maintains the core properties of the Wasserstein distance but with reduced computational overhead. For compactly supported measures, Bonnotte [Bonnotte, 2013] showed for instance that the two distances are equivalent. Again, it has been successfully applied in various domains, such as domain adap- tation [Lee et al., 2019], texture synthesis and style transfer [Heitz et al., 2021; Elnekave & Weiss, 2022]. generative modcling [Desppande et al., 2018; Wu al., 2019], regularizing autocncoders [Kolouri et al., 2018], shape matching [Le et al., 2024], and has even been adapted on Riemaniann manifolds [Bonet et al., 2024].

---

Published in Transactions on Machine Learning Research (06/2025)
The SW distance between two measures i and V can be written the expectation of the one dimensional Wasserstein distance between the projections of and on line whose direction is drawn uniformly on the hypersphere. It benefits from the simplicity of the Wasserstein distance computation in one dimension. In practice, computing the expectation on the hypersphere is unfeasible, so it is cstimated thanks to numerical integration. The most common method for approximating the SW distance is to rely on Monte Carlo approximation, by sampling M random directions uniformly on the hypersphere and approximating the integral by an average on these directions. Since the Wasserstein distance in ID between two measures of N samples can be computed in ( log N), computing this empirical version of Sliced Wasserstein has a runtime complexity of O(MN log N). This complexity makes it a compelling alternative to the Wasserstein distance, especially when the number N of samples is high.
As a Monte Carlo approximation, the law of large numbers ensures that this empirical Sliced Wasserstein distance converges to the true expectation, with a convergence rate of O Hr) This convergence speed is slow but independent of the space dimension. However, it is important to keep in mind that to preserve some of the properties of the SW distance, the number M of directions should increase with the dimension. For instance, it has been shown that for the empirical distance to almost surely separate discrete distributions (in the sense that if the distance between two distributions is zero then the two distributions are almost surely equal); the number of directions M must. be chosen strictly larger than the space dimension (Tanguy et al., 2023].
Classical Monte Carlo with independent samples is not always optimal, since independent random samples do not cover the space efficiently, creating clusters of points and leaving holes between these clusters. Tn very low dimension, quadrature rules provide efficient alternative mothods to classical Monte Carlo. On the circle for instance, the simplest solution is to replace the M random samples by the roots of unity 0 S k S M - 1}: since the function that we wish to integrate is Lipschitz, this ensures that the integral approximation converges at speed O(H) However, such quadrature rules are unsuitable for high-dimensional problems, they require exponential number of samples to achieve a given level of accuracy.
{e ak*
Another alternative sampling stratogy is to rely on quasi-Monte Carlo (Q.M..)) methods, which use deter- ministic, low-discrepancy sequences instead of independent random samples. Traditional Q.M.C. methods are designed for integration over the unit hypercube [0, 1|6. The quality of a Q.M.C. sequence is often mea- sured by its discrepancy, which measures how uniformly the points cover the space. A lower discrepancy correlates with a better approximation, according to the Koksma- Hlawka incquality [Brandolini et al., 2013]. Examples of dis-diccyypancy sequences for the unit cube include for instance the Halton sequence [Halton, 1964], and the Sobol sequence [Sobol, 1967], and different approaches have been investigated to project such sequences on the hypersphere. While quadrature rules are recommended for very small dimensions (d = 1 or 2 for instance), Q.M.C. integration is particularly effective in low to intermediate dimensions. A variant of low-discrepancy sequence is one where randomness injected in the sequence while preserving its "low discrepancy' property. Such a sequence called a randomized low-discrepancy sequence, and this is the foundation to randomized quasi-Monte Carlo (R.Q.M.C.) methods [Owen, 2019]. Q.M.C. methods do not only rely on cow-discepaancy sequences, but can also use point sets of given size directly optimized to have low-discropancy, such as s- s-Riesz point configurations on the sphere [Götz, 2003]. However Q.M.C. and R.Q.M.C. methods on the sphere have a strong practical downside: they suffer from the curse of dimen- sionality. Indeed the higher the dimension the harder it is to generate samples with Q.M.C. and R.Q.M.C. approaches. Moreover, the higher the dimension, the slower the convergence rate, and the more regular the integrand needs to be to ensure fast convergence. The recent paper [Nguyen et al., 2024] already proposes an interesting comparison of such Q.M.C. methods to approximate the Wessse--issr distance in dimension 3, showing that such methods could provide better approximations that conventional M.C. in this specific dimensional setting.
All the sampling strategies mentioned above are designed to provide a good coverage of the space. However, they do not take into account the specific structure of the integrand, which is the Wasserstein distance between the one dimensional projections of the two measures μ and v. More involved methods to improve Monte Carlo efficiency include importance sampling, control variates or stratification [Asmussen & Glynn, 2007]. Such variance reduction techniques strategics can also be used in conjunction with quasi-Monte Carlo
2

---

Published in Transactions on Machine Learning Research (06/2025)
integration. Control variates have been explored for Sliced Wasserstein approximation in [Nguyen & Ho, 2023] and [Teluc et al., 2024], showing interesting improvements in intermediate dimensions over classical Monte Carlo.
The goal of this survey is to provide a detailed comparison of these different sampling strategies for the computation of Sliced-Wassersterr in various dimensional settings. Tt is intended as a user-guide to help practitioners choose the appropriate sampling strategy for their specific problem, depending on the size and dimension of their data, and the type of experiments to be carried out (whether or not they need to compute numerous SW distances for instance). We will also look at the particularities of the different approaches, some being more appropriate than others depending on whether given level of accuracy is desired (in which case an approach allowing sequential sampling is preferable to one requiring optimization of a point set) or, on the contrary, a given computation time is imposed. We will mainly focus on sampling strategies which are independent of the knowledge of the measures μ and v, such as uniform random sampling [Asmusson & Glynn, 2007], orthonormal sampling [Rowland et al., 2019], low-discrepancy sequences mapped on the sphere [Halton, 1964; Sobol, 1967], randomized low-discrepancy sequences mapped on the spheres [Owen, 2019], Fibonacci point sets [Hardin et al., 2016] and Riesz configuration point sets [Götz, 2003]. For the sake of completeness, we also include in our comparison the recent approach [Leluc et al., 2024], which appears to be the most efficient among recent control variates approaches proposed to approximate Sliced Wasserstein.
The paper is organized as follows. Sec. 2 introduces some reminders on the Sliced Wasserstein distance such as its definition and some regularity properties. Sec. 3 explores all the sampling methods considered in this paper, hightlighting their theoretical guarantees, the conditions under which they can be used, and identifying which methods suffer from the curse of dimensionality. Then Sec. 4 provides comparison of each sampling method's experimental results different datasets. Finally, in Sec. 5 we offer detailed recommendations for choosing and using these sampling methods effectively in practice.
2 Reminders on the Sliced Wasserstein Distance
2.1 Definition
In the following, write { .> the Fuclidean inner product in Rd, || the induced norm, S (x E Rd |||||| 1} the unit sphere of Rd. For 0 E gd-1, we write To Rd + R the map x + (0lx), sa-1 the uniform measure over Sd-1, We also denote # the fush-forward operation !.
For two probability measures z and v supported in Rd and with finite moments of order 2, the Sliced Wasserstein Distance between μ and v is defined as
/ W3 Hout)/oaaa
SW = Eauu(sa- (ootuus
Tot)l SSa-1
(1)
This distance, introduced in [Rabin et al., 2012], has been thoroughly studied and used as dissimilarity measure between probability distributions in machine learning [Bonneel el al., 2015; Nadjahi, 2021; Kolouri et al., 2018], and more generally as an alternative to the Wasserstein distance. Its simplicity stems from the fact that the Wasserstein distance between two probability mcasures in one dimension has an explicit formula. Indeed, for two probability measures P1 and P2 on the line, the Wasserstein distance W2(P1, P2) can be written
W2(P),P2) = -1
(2)
where Fi and F2 are the cumulative distribution functions of P1 and p2, and FT' and Fz are their respective generalized inverses (see (Santambrogio, 2015) Proposition 2.17). For two one dimensional discrete measures
The furs-forward of a measure / Rd by an application T: Rd R² is defined measure T#p on Rt such that for all Borel sets E B(RR), #H)(B) ((T ((B)..
3

---

Published in Transactions on Machine Learning Research (06/2025)
N
N ogen this distance becomes
Σ P1 N Oap and p2 N k=1
Σ
N Σ Ea(N) - Vr(2a)
(P1.
(3)
where a and are permutations of [] which respectively order the sets {I1...., AND and {si..... yN on the line.
(a) On the left, we see the two discrete distributions (blue points) and v (yellow points). On the right, have their projections Ta#t (blue circles) and тe#t (yellow circles) along the direction One then takes the increasing ordering of no # and To# to obtain the corresponding matchings (green arrows) and computes the cost following Eq. 3.
100
3.5
0.2
4004
(b) On the left, have plot of 0 tTo #o4 in polar coordinates, with the distributions z and from Fig. 2.1a (top). The grey lines represent the angles where + W (mo #w) is not differentiable, the magenta line the line of angle right, we have 1D plot of + ( Ma #V), here the hashed area represents SW (1, v) and again the vertical grey lines represent the values where (TT #o#))
_ and the blue dot a specific value W2 (me #e with the same angle. On the
differentiable.
Figure 2.1: On the first row, Fig. 2.1a illustrates the computation of W2 (motH mo# for fixed 0. On the second row, Fig. 2.1b gives geometrical illustration of SW2 (p, v2) with u.v taken as in Fig. 2.la.

---

Published in Transactions on Machine Learning Research (06/2025)
N >
As a consequence, the Sliced Wasserstein distance between two discrete probability measures
day
N
and v N Σ Eye on Idd (i.e. with (#.)) N. (yx (- Rd) can be rewritten:
((0k Vrona, (()), (H, N Σ Se-1 ((Foo()) - O)) N ga-1
1
1
E/
(4)
where dg and To denotes respectively permutations which order the one dimensional point sets ((zx, 8))x=I,...,N and (Wy,, 0))==,,,)). Fig. 2.1 illustrates the computation of (2 To for two discrete measures in two dimensions (Fig. 2.1a), and shows how this quantity varies when 0 spans [0, 2 (Fig. 2.1b).
Since the permutations dg and To depends on the direction 0, the integrals in Eq. 1 and Eq. do not have closed forms. For this reason, practitioners rely on Monte Carlo approximations of the form:
M
(Ton (1) Mra, (k), 0,))
(5)
NM
k=1j=1
where .1 are i.idd and follow uniform distribution on the sphere. Classically, the convergence rate of such Monte Carlo estimations to SW ((7))) [Hammersley & Handscomb, 1964]. In this context, it is natural to question the optimality of sampling methods to approximate SW efficiently different scenarios.
2.2 Regularity results on 0 + ( #H, To#t)
The efficiency of sampling strategies used in numerical integration is highly dependent on the regularity of the functions to be integrated. For this reason, in the following we give some properties of the function (Fig. 2.1b):
f - W2 (To##, #o#
(6)
the hypersphere Sd-1, We first look at classical regularity properties of f. Proposition Lipschitz on 5t-1
Proof- Let a and be two probability measures with finite moments of order 2, and 01, sd-1, The triangular inequality W2 yields
[W2 #1, ne, #v) We(To, ne, #v)l S W2 #V. ne, no, #p) W2( #V, no, #v).
We also have
Wi Hu, #1,
inf E [[K0,, Y) < inf E [K0] 02, 02, xy < 110 xx- X-p
X-uuY~w mm
We show similarly that WT Hu, mo #V) < ||0 (Vxx=
(|xx))
Thus
(T #p, No,
Wy (To, # ne, #v)l < ||0 -2ll
Since ff is Lipschitz continuous, it is differentiable almost everywhere. However the previous result does not give us the set where is non differentiable. In the following we give more complete proof when μ and are discrete following the notations introduced in Sec. 2.1.
Proposition 2 When a and finite discrete measures, f piecewise C000 (Cpc) and Tipschitz Sd-1
5

---

Published in Transactions on Machine Learning Research (06/2025)
N N > 8z and V = N
N
Proof. For discrete measures
by on RR, f can be rewritlen as
f(0) min fo f0), where (0) aeEN
Σ (ok Mo(k) 10)2
(7)
where En is the set of permutations of []] We assume that the { (resp. {yy}) are all distinct. In the following, we study the regularity of f as function of Rd and deduce the regularity properties of its restriction Jise-1. Observe that each fo defines a quadratic function on IRd and f, as a minimum of a finite number of such functions, is continuous and also piecewise C on Rd. Since is continuous on RR, its restriction to gd-1 is also continuous. To show that this restriction to Sd-1 is also in Cpu it is enough to observe that the set of points of Rd where f is not differentiable is included in the finite union of hyperplanes (U,,, Span(a,, U (pannS) - )))+), since these hyperplanes are the locations where the minimum in Eq. jumps from a permulation a to another one (see Fig. 2.2 an illustration of those hyperplanes). Each of these hyperplanes intersect Sd-1 on great circle, and we call U/ the sphere minus this finite union of great circles. The open set U (which is dense in Sd-1) can be written as the union Uk V& of finite number of connected open sets V,, such that on each Vi, the permutation a which attains the minimum in Eq. 7 is constant and unambiguous. We write this permutation di. On each Vi, fisa-1 = fa. thus is C on Vi and its derivative can be obtained as the projection of Vfa; on the hypersphere. For U, writing dg the permutation which attains the minimum in Eq. for the direction 8, this derivative can be written
(d-1)](0)=
Σ (TTR (oo( Noo(k))) - (ax - -oocu)10)20)
Since these derivatives are upper bounded on the compact set Sd-1, it follows that f is also Lipschitz on 5d-1.
In the case where several x4 (or yj) are equal, several of the functions fo coincide. For instance, x1 22 the values of (1) and (2) can be exchanged without modifying fr. By eliminating all the redundant functions, make the same reasoning as before to show the same regularity results on f. In this case, all the pairs (z) with zj should be removed when constructing the set of great circles dividing the hypersphere.
The following proposition will also be useful in the next sections. Proposition 3 3: f E H (Sd-1), where, for a E N, the Sobolev space (Sd-11) is defined as [Hebey, 1996]
H = {h E 13(5d-1) æ L² 15l sa).
with j multi-index and اiا the partial mixed derivative of order [ on 1d-1
Proof. We have seen previously that f is continuous and piecewise C piecewise quadratic to be more precise. Thus its weak derivative is piecewise linear with discontinuities on a finite union of hyperplanes, which is L2
3 Reminders on sampling strategies on the sphere and their theoretical guarantees
In this section, we present the different sampling methods for numerical integration Sd-1 considered in this paper, before comparing them experimentally in Sec. 4. This paper addresses three main types of sampling: random sampling, discrepancy-based sampling, and a control variate approach. The first type includes the classical Monte Carlo (M.C.) method ([Hameersley & Handscomb, 1964], [Lemieux, 2009]) on the sphere and its variant called orthonormal sampling [Rowland et al., 2019]. The second one relies on a concept called the discrepancy ([Lemieux, 2009), & Pillichshammer, 2010)) of a point set, which represents the number of points in a unit of volume, and can be divided into two categories: low-discrepancy sequences (or digital nets) and point sets (or lattices). Among the former category, also investigate

---

Published in Transactions on Machine Learning Research (06/2025)
33
-05
-10
-1.0
Figure 2.2:
Illustration of the open subsets Uh=1 V/x and their intersection with the hyperplanes (Ji,Span(z;; - ( +)) in the specific case of two measures made of two diracs
" 8z: with 31, 42 (1,0,0)* 1,, and Σ Sy with v1 (0,0,1)", (0,0,
i=1 The hyperplanes divide the sphere into the colored sections where and To are constant.
method based on a spherical sliced- Wasserstein type discrepancy [Bonet et al., 2023]. The last type details a control variates method [Lemieux, 2009] using spherical harmonics [Müller, 1998] for this purpose [Leluc et al., 2024]. We will also determine which method, and under which conditions, is theorctically suitable based the regularity properties established in Sec. 2.2. Tab. presents a taxonomy of all the sampling methods explored in this paper. It details which method's convergence rate result is independent from the dimension (i.e. the dimension does not appear in the asymptotic rate), which one can be computed independently (i.c. cach sample can be generated independently from the others), and which one can be computed and stored in advance.
Sumpling types
Dimension independence
Independent computation
Possible pre-computattion
Random Sampling
Uniform Sampling
Oothoonlaal Sumpling
Based discrepancy
Riesz Point Set Riesz Point Set. Randomized
Fibonacci Point Set Fibonacci Point Set Randomized
Sobol / Sobol Randomized mapped on gd-1
Halton/ / Halton Randomized on Sd
Spherical Sliced Wassenteia Discrepancy
Control variates
Spherical Harmonics Coutrol Variates
Table 1: Taxonomy of the three types sampling methods investigated this paper.
Tab. gives summary of the convergence rate and computational complexity of each sampling method
explored in this paper. In this table nı of -MM
3.1 Random samplings
We first explore classical strategies for randomly generating points on the sphere: uniform sampling [Ham-
mersley & Handscomb, 1964] and orthonormal sampling [Rowland et al., 2019). These strategies are the
7

---

Published in Transactions on Machine Learning Research (06/2025)
Sampling types
Theoretical couvergenoe
Time complexity
Space complexity
Random Sampling
Uniform Sampling
00//M)
((M) ((M)
O(M)
Orthonormal Sampling
None
O(M)
Based on discrepancy
Riesz Point Set / Riesz Point Set Randomized
1/M on S' Not applicable otherwise
O(M")
((M)
Fihonacci Point Set / Fihonacci Point Set Randomized
Not. applicable
((ar)
O(M)
Sobol / Sobol Randomized mapped on 5 Halton Halton Randomiaed on S4-1
None
0 (Adlog (MF)
O(M)
Nome
0 /Mllgg (M)
O(M)
Spherical Sliced Wasserstein Discrepancy
None
O O(Aflog(
O(M)
Control variates
Spherical Harmonies Control Variates
(L/(mNNV))
O(M))
O(M)
Table 2: Convergence rate, time complexity and spacial complexity (w.r.t the sampling number) summary of the sampling methods studied in this paper.
most commonly used for estimating SW2, and their convergence rates do not depend the dimension of the input measures.
3.1.1 Classical Monte Carlo
The classical Monte Carlo method uses uniform random sampling to generate the projection angles. For (0x) MBM))eN- i.i.d samples of sd-1 2, we write the Monte Carlo Estimator
M Σ ((0)) with M EEN.
XM = M
(9)
The law of large numbers ensures that XM converges to (w, Eg [f()]] as M goes to infinity. Moreover, the rate of convergence for this unbiased estimator is given by
V[X1 M
VVxM
(10)
VM
where o² v((() ga-1 (00)saa-100) (W, D) < +oc This convergence rate in Eq. 10 does not
depend on the dimension of the input measures. In order derive confidence intervals for SW (1, v), we can rely on the Central Limit Theorem [Fischer, 2010] which states that
VM XM SW2(k, C + NO,,
M-+00
This allows us to compute confidence intervals for SW3 (1, v) by using the quantiles of the standard normal distribution. This means that for M large enough, P [Xm SW2 (p,
0q1-0/2
091-0/2
+1
VM
/M
M-+00
a, with a in [0, 1] and qı-a/2 the quantile of level 1 - a/2 of (O, 1). One strategy for choosing M is taking M such that 291-a/2 VM s with 2 0 chosen precision. The value of being unknown, possibility is to plug a consistent estimator of o2, such as
Σ
f(0.)² -x²
M
Xianliang & Zhougyi [2022] provide an alternative criteria for choosing M, however it is quite impractical as it requires compute the Wasserstein distance between y and v.
IIn practice to simulate random variable ~ one takes normal random variable 2 Ni0,, La) o 0 and chooses 12T 84-1 [Asmussen Glynn, 2007].
8

---

Published in Transactions on Machine Learning Research (06/2025)
3.1.2 Orthonormal sampling
A variant of the uniform sampling covered in Sec. 3.1.1 was introduced by [Rowland et al., 2019], which presents a simple variant for the previous Monte Carlo estimator XM by sampling random orthonormal bases. This method is inspired by variance reduction techniques known as stratification [Lemieux, 2009]. Let O(d) be the orthogonal group in R For POP) ~ (O())) denoting OI...., 0s all the columns of M
1
the matrices .,...., , we define Yu = M E f(0.). It is easy to show that each 0 follows the uniform
distribution on $d-1 [Rowland et al., 2019]. As a consequence, the estimator Yι is still unbiased. Although it is not possible to show that YM has a smaller variance than XM in general, this estimator is most of the time more efficient than XM our experiments and show an equivalent or better rale of convergence practice. This might be due the fact that the diversity of the samples is increased by the orthonormality constraint.
Remark 1 Other fully random point processcs on [0, 1]2 or S suitable for Monte Carlo integration are studied in the literature. Among them, we can mention Determinantal Point Process (DPP). Recent works, such [Feng et al., 2023], have proposed DPP methods directly on the sphere S Unfortunately, due to the lack of publicly available implementations, we could not experiment efficiently with these methods.
3.2 Sampling strategies based on discrepancy
We examine this section two different types of deterministic sampling based on discrepancy: low- discrepancy sequences (digital nets) and low-discrepancy point sets (lattices). They were developed to replace random sampling, expecting to have better convergence rate than the classical Monte Carlo method.
3.2.1 Low-discrepancy sequences
Quasi-random sequences, better known as low-discrepancy sequences (I..D.S),, are sequences mimicking the behavior of random sequences while being entirely deterministic. To date, these sequences are only defined on the unit hypercube [0, 11d. We introduce below a first definition of discrepancy ([eemieux, 2009],[Dick Pillichshammer, 2010)).
Definition 1 The discrepancy of a set of points P = {uI,.... uM } in 1j is defined as
PNl
(M) sup TET
M
d
where |A| denotes the cardinal of set A, 2d is the d-dimensional Lebesgue measure and I = ([] a4 byl 0
i=1
a < 6, 1}. The star-discrepancy Dy (P) is defined the same way with I"
111 [0, bil
0 6, 1}.
We can now provide a definition of a Low discrepancy sequence (L.D.S.).
Definition 2 : Let (um)men- be a sequence in [0, 1jd. Denoting Pa = . uM for any M E NN, u T.D.S.
Dis --0 0.
M++-00
The notion of discrepancy is important because it is related to the crror made when approximating an integral on the hypercube by its Monte Carlo approximation. This relation is made explicit by the Koksma- Hlawka inequality ([eemieux, 2009]; [Dick & Pillichshammer, 2010]; [Brandolini et al., 2013]). This inequality requires introduce the notion of Hardy-Krann variation V,, of function h on [0, 1] [Aistleitner et al., 2016], which is out of the scope of this paper, but can be broadly understood as measure of the oscillation
of h on the unit eube [0, 1]
9

---

Published in Transactions on Machine Learning Research (06/2025)
Proposition 4 (Koksma-Hlawe) inequality) Let h : [0, 1]d - R have bounded variation V₁, on [0, 1J" in the sense of Hardy-Krause [Aistleitner et al., 20. Then for (ui,..., war) a point set in [0, 1]4, we have
M h(u))) - - -J h(z)d10(()) S (u
(11)
M
The proof of this inequality and basic results on discrepancy theory can be found in(Kuipers & Niederreiter, 2012] and [Dick Pillichshammer, 2010]. Eq. 11 shows that the absolute error made by the Monte Carlo approximation is upper bounded by a term depending only on h and the star discrepancy. Compared to the Central Limit. Theoreom, this inequality is not probabilistic and not asymptolic, the bound being valid for every M IN An important limitation is the term Vh. which is impractical to compute directly. When d = 1, this term is exactly the total variation of h, but in general, it is only upper bounded by the total variation. In the case of our function f involved in the estimation of SW, Vy < +oo holds since f is Lipschitz continuous. Another limitation of the previous bound is that the rate of convergence of the star discrepancy Di of a sequence is most of the time not explicit and difficult to compute [Owen, 2005].
Nevertheless, this proposition ensures that the rate of convergence of the star discrepancy of a sequence is better than O (())) ): for M large enough the approximation of the quasi Monte Carlo approximation using this sequence will oulperform the one of classical Monte Carlo.
In the following, we present two L.D.S. defined on the unit square [0, 1] and see how their star discrepancy decreases with M. We then focus practical methods to map these sequences from the hypercube to the hypersphere Sd-1.
3.2.1.1 Halton sequence
The Halton sequence (ut)ee E (Rd)N [Halton, 1964] a generalization of the von der Corput sequence [van der Corput, 1935]. In the following, we write, for any integer ci(i) the coefficients from the expansion of i in base b, and we define the radical-inverse function in base b
400
Σ 9016-''), WEN. 1-
The Halton soquence in dimension then defined as
w (0b, (h. (,)
where by is chosen the i-th prime number.
3.2.1.2 Sobol sequence
This sequence uses the base b= 2. To generate the j-th coordinate of the i-th point uy a Sobol sequence [Sobol, 1967], one needs primitive polynomial of degree n, in Z/22[X],
n 01, x + 02,
+ any -1j + 1.
This polynomial is used to define a sequence of positive integers (mz.j) by recurrence, with +z/27 the inner law of 2/22:
mkj = 201, M +2/22 2 02,,; +2/22 ... +2/22 2 my-nnid +2/22 m3-nn,
The values maji for 1 k < nřj can be chosen chosen arbitrarily provided that each odd and less than 2². Then one generates what is called direction numbers:
ma.j
Uk.j =
2
10

---

Published in Transactions on Machine Learning Research (06/2025)
The j-th coordinate of u, is then obtained as
0:0
Wi = Σ cs (i)oR..
3.2.1.3 Convergence rate of Halton and Sobol sequences Both sequences (Halton and Sobol) have a star discrepancy which converges to (which means that they are indeed L.D.S.). The convergence rate is given by the lollowing property [Niederreiter, 1988] [Owen, 2019].
Proposition 5 Let (Um)mEN- be either the Halton sequence or Sobol sequence in [0, 1] Then for M 2 1. wc have
(M)
D (ui,..... uas) ca
M
where ca is constant that depends only on the dimension.
Thanks to Eq. 11, for any function h such that V, +00 (which the case for our function f), this implies a convergence rate of the Monte Carlo estimator using these sequences in O (("or(M)") which means ((M-1++) for every > 0. This convergence rate is better than the one of classical Monte Carlo with i..d. sequences, even if the rate of convergence slows down when the dimension increases, because of the term log(M)".
Remark 2 Note that L.D.S. are designed mimic the behavior of random uniform sampling in [0, 1)4 while being completely deterministic. This deterministic behavior leads to patterns in the sampling; because of those patterns, the higher the dimension, the harder it is for those to fill the "gaps" in [o, 1] Morcover, the term log( MM implies that one needs M to be very large (exponential) to get the same level of space coverage in high dimension than in low dimension.
Remark 3 : Observe that both for Sobol and Halton sequences, generating M values has a complexity in OMllgg? (M)), where b is the base (or smallest basis for Halton) chosen.
3.2.1.4 L.D.S. on the sphere
To our knowledge, there is no true L.D.S. on the unit sphere Sd-1 for d 3, this question remaining an active research area. Practitioners typically map L.D.S. from the hypercube the hypersphere, using one of the methods described bclow:
Equal area mapping |Aistleitner et al., 2012]: this method is only defined for mapping points in the unit square to S Denoting (2) E [0, 1f2, one gets a point = (1 1 - 222) on S2 with:
v(1.3) (Vi (o)(s)), VI (Psin()))))) .)), n, 3 10.1|
(12)
Spherical coordinates [Arfken et al., 2011]: This method maps the points from an L.D.S. in 1]d-1 to Sd-1 by using the spherical coordinates. Unfortunately, we found that the resulting sampling is usually not competitive compared other sampling methods.
Normalization onto the sphere [Basu, 2016]: An L.D.S. is generated in the d-hypereube [o, 1]d and mapped to Rd using the inverse cunulative distribution function of the standard normal distri- bution (separately each dimension). Then each point in the resulting sequence is normalized by its norm to map it onto Sd-1.
Specific case of $2.
In the specific case of $2 it has been shown by Aistleitner et al. [2012] that if is an L.D.S in [0. 1]2 and & the equal area mapping defined in Eq. 12, the spherical cap discrepancy DL..M (P(P) (see definition 3 in the next section) of the mapped sequence in O (Zz17) However, their experiments showed that the correct order seems rather to be O (Mgg(()) for 1/2 < e
M/74
11

---

Published in Transactions on Machine Learning Research (06/2025)
3.2.2 Deterministic point sets S#-1
This section details different methods to design well distributed point sets on 5 Contrary to the L.D.S. defined above, these point sets are defined directly on the sphere, in order to be approximately uniformly distributed on gd-1, To measure this uniformity, we can rely on the notion of spherical cap on the sphere: a spherical cap of center c e Sd-1 and t E [-1, 1] is defined as
C(c, t) (x S (z, t}.
(13)
In other words, a spherical cap is the intersection of portion of the sphere and a half-space (sce Fig. 3.1 for illustration).
Figure 3.1: Illustration of spherical cap on S². The circle represents the intersection of the plane (2,,) with the sphere, and the purple colored area is the cap C(c, t) as noted in Eq. 3.2.2.
To the best of our knowledge, there is no equivalent to the Koksma-Fflawka inequality for the sphere in full generality [Brauchart, 2011]. A sequence of points {u,} on 5d-1 is said asymptotically uniformly distributed on 5d-1 for every spherical cap C, the proportion of points inside the cap, converges to the measure of the cap (C-1 It can be shown that this assumption is equivalent to assume that for every continuous function h, the Monte Carlo approximation 1 EM h(u))) converges to E [k()]-
M
In order to gel a non asymplotic notion of the uniformity of a point set on Sd-1, we can rely on different notions of spherical cap discrepancy on the sphere, defined as follows.
Definition 3 The spherical cap max-discrepancy of a point set PM size M is defined as [Marzo & Mas, 2021]:
{
PM (C(c,t) -c-,c(())))
Dman PPaa(PM))
sup -E54-1,E([11
M
The spherical cap L2 di-diccryyancy of a point set PMr of size M is defined as [Brauchart, 2011]:
(PMs) - [[
PM PM
8a (C(e, daa- dea- (c) }
gd-1
M
where t(c,t) is spherical cap of center and height f.
Again, the idea is compare the proportion of points Par that fall inside spherical cap with the measure of the cap. This comparison is done for all possible caps on the sphere, and Dmax represents the worst error over all possible caps, while D2 represents the average squared error over all possible caps.
12

---

Published in Transactions on Machine Learning Research (06/2025)
When using Q.M.C. the hypersphere to approximate the integral of functions h, another notion often used in the literature is the worst-case (integration) error (W.C.E.) on Banach space of functions, which is the largest possible error made by the method on the space. For instance, on (Sd-))),
Definition 4 For PM- {u1...., uM): for a E N
1 M
WCF PM, H (5
Σ m=
M(mm)
[ Sd-1(SA-1) SS h(w)dss-1(w))
sup (EH-(S(-11
M
Under some regularity condition, a sufficient and necessary one being a 2 + for (a( Brauchart & Dick [2013] show that optimizing the spherical cap La-discrepancy is equivalent to optimizing the W.C.E. thanks to the Stolarsky's invariant principle [Stolarsky, 1973]. the case of our function f, we have seen that f regular enough in the specific case of St, since f (St with dimension larger than this result does not hold anymore since f does not belong to any Sobolev space (S4) with > 1.
+ 3.2 However in
3.2.2.1 Fibonacci point set on S
Denoting o the polar angle and X the azimuthal angle forming the geographical coordinates (p, x); we retrieve the Cartesian coordinates (æ, y, 2) using the spherical coordinates (see Fig. 3.2 for an illustration). Noting o Ity5 the golden ratio, the m-th point = (Pmm Xm) of the Fibonacci point set is given by
2m
Pm arccos
2M +
2mn -2
Xm
It is simple and efficient way, convergence rate wise, to generate points S for the quasi-Monte Carlo method but it is only defined on S². The complexity of the generation is linear in M, and according to [Maqques, 2013], the corresponding convergence rate for the W.C.E. and the Læ-spherical cap discrepancy is in O( () For an extensive list of other popular point configurations on $2, see [Hardin et al., 2016].
Z
sin cos( y sin (p) sin(2); cos(x)
y
P
Figure 3.2: Illustration of the spherical coordinates in R for points on the sphere S?
3.2.2.2 Equi-distributed points generated by the discrete s-Riesz energy
Another classical way to define equi-distributed point sets on the hypersphere is to rely on optimization. In such methods, the point set PM is defined as the minimizer of a cortain energy functional Eg,
P
arg Es uM).
13

---

Published in Transactions on Machine Learning Research (06/2025)
The most common energy functional is the s-Riesz energy, which is defined as follows.
Definition 5 For & 2 and Par = {W1,... uMr) a set of points on Sd-1, the s- Riess energy of P is defined as
if 8 > 0,
uzll
(.(PM)
log
if 8 = 0.
||u Ugl
Engy
The resulting point set is called a minimal s-energy configuration. The s-Riesz energy can also be defined for < O, in this case the point set Pas is obtained as the maximizer of E, = Σ |lu uyll [Brauchart,
2011]. Minimising E, is non trivial, the functional being not convex, and the problem becomes more complex when the dimension increases. Minimal energy configuration points for Eg called Fekete points and it is known that for 0 S 8 < d, these sets are asymptotically uniformly distributed with respect to the normalized surface measure sd-1, which means that Monte Carlo estimates using the Fekete points converge to the integral against 8d-1 [Mazo & Mas, 2021].
The spherical cap L=-discrepancy
a point configuration is minimal and only if the sum of distances in the configuration is maximal. This would correspond maximizing a s-Riesz energy for s -1 [Brauchart, 2011]. However, the link between the configurations of minimal s-Riesz energy and the ancies of these configurations is in general not straightforward, see [Brauchart, 2011], [Marzo & Mas, 2021], [Götz, 2003]. For 0 S < d, and PM a minimizer of size M of the Riesz s-energy on Sd-1, the authors of [Marzo & Mas, 2021] show that
2, discrep-
D Dmox (PM) & O max [M -Maaaai 2 M 2(d-014)
0, but the speed of convergence degrades with the dimension which means that the uniformity of these configurations likely to suffer from the curse of dimensionality. Fig. 3.4 shows an example of s-Riesz points and Fibonacci points on S2 with 500 points.
This implies that Dmaz (PM)
M+++0
Remark 4 Since computing Riesz point configurations involves optimization (with non linear complex- ity), the time neoded to generate those points can be impractical. Note that generally the generation of the 8- &-Riesz configuration points has a runtime complexity of O(TM²), where T the number of iterations of the optimization loop.
In the specific case of 51, the Fekete points are unique up to rotation, and are the M-th unit roots (see [Götz, 2003] and see Fig. 3.3 for an illustration):
{e
x k O M M-1]
This explains why for 2D discrete measures, uniform grid on St gives better results than any other sampling method for computing SW2 as we will see in Sec. 4.
3.2.3 Random Quasi Monte-Caria
The principle of Randomized Quasi-Monte Carlo (R.Q.M.C.) methods is to reintroduce stochasticty in Q.M.C. sequences. Indeed, Q.M.C. methods such as the ones described in Sections Sec. 3.2.1 and Sec. 3.2.2 are deterministic. For a given M, the estimator given by one of these methods is always the same. As such, we cannot easily estimate the error or the variance of the Monte Carlo approximation. Besides, while results such as the Koksma- Hlawka inequality ensures that they converge at a certain rate, the different quantities involved in the incquality are much more complex to estimate than the one involved in the Central Limit theorem. Random Duasi-Monts Carlo methods were especially designed to recover this ability estimate the error easily. These sequences usually defined on [0, 1]d.
14

---

Published in Transactions on Machine Learning Research (06/2025)
102
175
502
.075
-199
Figure 3.3: Plot of the 10-th unity roots, i.e solutions to the equation 210 1.
Illustration of s-Riesz points (on the left) and Fibonacci points (on the right) on S², with 500 Figure 3.4:
points for both configurations.
Definition ([Owen, 2019]) Let (ūz))21 be sequence of points in [0, 1jd. said to be suitable for R.Q.M.C. if Vi, 02 U(0, 18 and if there exist finite e > 0 and K > such that for all M 2 K,
Dir(Ps)) log"(M) M
1,
where Pas = fa...... km}.
M
Denoting XM = M Σ M((.)) the empirical estimator of (E [h(0)), the assumption â; M((0, 1)4) implies that XM is unbiased. Besides, the previous inequality implies that if (0 is suitable for R.Q.M.C., then the variance of XM bounded by v2 (M)) For functions such that V, < co, this yiclds a convergence rate in O
M2
similar to the one of low discrepancy sequences.
Once a randomization method is chosen (such that it provides suitable R.Q.M.C. sequences), the process can be repeated several times to obtain K independent random estimators X1. ......... X of Eg-sa [b(0))
K
The agregated estimate XM.K
One of the key advantages of this approach is that this variance (or confidence intervals) can be estimated by the empirical variance of the K independent estimators.
x has variance decreasing in
K
R=1
There are several ways to generate sequences from low discrepancy sequences on [0, 1] in order to make them suitable for R.Q.M.C.. One of the most simple methods consists applying the same random shift all points in the sequence, and taking the result modulo 1 componentwise [Lemieux, 2009]. More involved methods, such Digital shift Scrambling, described in [Lemieux, 2009] and [Owen, 2019].
15

---

Published in Transactions on Machine Learning Research (06/2025)
However, to the best of our knowledge, there is proper R.L.D.S. on the sphere, stated by Nguyen et al. [2024]. In practice, R.I.D.S.. on the unit cube are mapped onto the sphere by the methods described in paragraph 3.2.1.4. Another possibility, as done in Nguyen et al. [2024], is to generate a random rotation matrix and apply directly on point configurations on $d-1, such as the oncs described in Sec. 3.2.2.
3.3 Spherical Sliced Wasserstein
A sampling method based on a. Serred--assessteein type discrepancy on the sphere S was developped by Bonct ot al. [2023] for d 2 3. We denote Ca22 the set of great circles of sd-1, a great circle being the intersection between a plane of dimension 2 and 5 [Jung et al., 2012]. The authors of Bonet et al. [2023]
define a pseudo distance, called Spherical Sliced Wasserstein distance, between two probability measures θ, = defined on Sd-1:
( W2(
(14)
where for all x E Sd-1, =c(x)) arg min dga-1 (z. with dga-1(x, (z, v) arcos( (2r, v)) [Fletcher et al., 2004] and { is the uniform distribution over Ca.2-
As shown in Bonet et al. [2023], this distance can be used to sample points on Sd-1 by minimizing SSW2 M
Σ 50, and the uniform measure E = sa-1 on sd-1, To this aim, for CI...., Ct I. independent great circles, they approximate SS,2(0, =) by its Monte Carlo approximation 1 Σ (TC To #3) Then, they note that TC, #s 31 [Jung, 2021] for each and (0, I 1=1
between a discrete mcasure θ = M
derive closed form for (TC,#O, s1) based Delon et al. [2010]. The final distance SS,?(0, then be optimized with respect to the point positions 0, with a projected gradient descent.
Remark 5 There are cases in which SSW is metric:
Based on Quellmalz et al. [2023], SSW is a metric between any two probability measures on S?.
A result from Liu et al. [2024] also shows that SSW is a metric between any two absolutely continuous probability measures with continuous donsity functions on 5d-1 for d 2 3.
Remark 6 Noting T' the number of iterations for the gradient descent algorithm, and L as above, then the time complexity this method is in L(L.
Remark Notice that SSW's form is similar to the La-spherical cap discrepancy, where instcad of averaging the 'error' made by the sampling spherical cap, it averages the "error' made by the sampling a great circle.
3.4 Variance reduction
All methods described so far are based on the idea generating points on the sphere such a way that these points are sufficiently well distributed be used for Monte Carlo integration, and ideally yield faster convergence than M.C. with i.i.d. sequences. These point sequences point sets are defined independently of the function to be integrated.
More involved approaches, such as importance sampling or control variates, use the knowledge of the function to be integrated to improve Monte Carlo estimators by decreasing their variance. Recently, two control variates based methods have boen developped to estimate the Slicod Wasserstein distance. control variate is a centered random vector Y E IR, easy to sample, with finite second moments. Assume we want to estimate (f Writing 01...., ,0 i..d. samples of 0 sa-1 and Yi...., YM Al independent copies of the random centered vector Y, consider the following estimator
M
Σ f(())) - -8"
M
16

---

Published in Transactions on Machine Learning Research (06/2025)
where 3 E RP is constant vector to be determined. The variance of this estimator is proportional to Var(f(0) BTY). It follows that if we write f the parameter minimizing this variance, then the pair (L(f(0)),3*) is solution of the least square problem
min El((()) - (G
An empirical version of this quadratic problem sample (01,.... 0xs) writes
(E(f(0)) Mr BM) = arg min | - --1 - 1 C
(15)
where F (f(01)) 1M (I....,,)) 1)T RM, and Y Y= (Y,)) i= E RMP..
Recently, Nguyen & Ho [2024] introduced a Sliced Wasserstein distance estimation using Gaussian control variates and Leluc et al. [2024] developped a method using spherical harmonics control variates. We focus only on Lelue et al. [2024] here, since their method yiclds much better experimental results. In their work, Lelue et al. [2024] chose Spherical Harmonics [Müller, 1998] as control variates. Spherical harmonics are functions which form orthonormal basis (d:) of the Hilbert space T.2 8a-1). In this setting, the random variable Y is thus chosen as Y (d)(0)): with 0 sa-1. In practice the number p
is chosen as p = Ln,d = Σ (d, 2(), the number of spherical harmonics of even degree up to 2n, with N(d, n) (2n d 2) (n+d-3)! (d- 2)hm! the number of spherical harmonics of degree n in dimension d.
Leluc et al. [2024] then compute the solution (SHCVM Sa of (15) on a sample (01.... 0xr) and use the control variates estimator SHCV? M.n as estimator of the (squared) Sliced Wasscrstein distance.
They prove the following convergence property.
Proposition 6 Let v.v be two discrete measures in Rd with finite moments of order 2 and let d 2 2. For any sequence of degrees n = (Mas) such that nAs = o M (2(d-1)))) as M - t00, we have
SHCV SW2 - Op nM1122
(16)
where the notation X,, (p()n) means that the sequence is stochastically bounded 3,
Notice that since nM = o M (2(d-))))) in high dimensions d the global convergence rate is similar to that of the classical Monte Carlo method described in Sec. 3.1.1.
4 Experimental results
This section presents experimental results from all the different sampling strategies presented in Sec. 3, on a variety datasets. To provide representative results, we select datasets spanning a range of dimensions going from 2 to 28 28. Those include toy dataset and three 'real-life" ones. We first present results on Gaussian mixtures in the following dimensions (2, 5, 10, 20, 50), the six ground truths (true distances) are estimated using 108 projections. Secondly,
show some dimensionality reduction results 12 different datasets of persistence diagrams (for the case of 2 dimensional discrete measures). Then show some convergence results in the specific case of measures in 3 dimensions. Specifically, we compare different estimations of the Sliced Wasserstein distance between 3D point clouds taken from the ShapeNetCore dataset, see [Chang et al., 2015]. Finally we compare different dimensionality reduction results on the MNIST datasct [LcCun, 1998]. For the experiments on the Gaussian mixtures we compare the listed stralegies with the following sampling
The notation X. Op(an) that > there exists finite 0 and N > such that P[[X.] > Ka.] <
all > N.
17

---

Published in Transactions on Machine Learning Research (06/2025)
numbers (100,300, 500, 700, 900, 1100, 2100, 3100, 5100, 6100, 7100, 8100, 9100, 10100). Otherwise, we use the following sampling numbers (100, 1100, 2100,3100 4100, 5100, 6100, 7100, 8100, 9100, 10100). Tab. 3 displays the acronyms of all the sampling methods compared in the following experiments. For each sampling method from Tab. 3, thore are two variants finishing with the term "Area Mapped' and two variants finishing with the term "Normalized Mapped'. The first one means that we applied the equal area mapping detailed in paragraph 3.2.1.4. The second one means we normalize each point generated by those methods, this normalization method is also detailed in paragraph 3.2.1.4.
Name
Legends
Dimensions
Riesz Randomized
R.R.
20,
Uniform Sampling
U.S.
2,3,5,10,20,50
Othoruormal Sampling
O.S.
1,,,,,0
Halton Area Mapped
H.A.M.
2,3
Halton Raudomized Area Mapped
H.R.A.M.
Halton Normalized Mapped
H.N.M.
5,10,20,50
Halton Randomized Normalized Mapped
H.R.N.M
5,10,20,50
Fihonacei Point Set
F.P.S.
Fibonacci Randomized Point Set
F.P.P.S..
Sobol Area Mapped
S.A.M.
Sobol Randomized Area Mapped
S.R.A.M.
Sobol Normalized Mapped
S.N.M.
5,10,20
Sobol Raudomized Normalized Mapped
S.R.A.M.
5,10,20
Spherical Harmonics Control Variates
S.H.C.V.
3,5,10220
Spherical Sliced Wasserstein Randomized
S.S.W.R.
3,5,10,20,50
Table 3: For each method used in this experimental part, associated acronym, and list of dimensions where this method is used.
4.1 Implementation of the sampling methods
This section provides details the implementations used for the sampling methods, and specifies how the parameters are set. The implementations used are grouped and are available here https://githin /omm Keanu Seanu-Sisouu/SK--Sappp Guide.
For both methods we used python included functions to sample a Gaussian variable and to sample orthogonal matrices in dimension. For sampling orthogonal matrices we use the following python library scipy.stats.ortho group https://docs.sct org/ toc//sccip/rrteerenc//ggeeeee wttpaa/eeeeeeeeeeeeeeeeeeeeeeeee scipy. stats. stats.ortho ort.o cachoogroup.hii html.
Classical M.C. methods:
Halton & Sobol sequences: In dimension 3 and less, we use python implementatiions from the library scipy.stats.quua (attps://docs wwtpsassccccccccccoeeeieoooeeeee ottp/////.ccc///teeeee////////// stats. qmc. Hal.on html & attps://docs.sceeeee stats.qme Soo.c.oblll..e. version. For high dimensions, we Leluc et, al. [2024]|s implementation available here https:
http/////cccceeeeeeeeeee///eee/e As for the parameters we set 'scramble' to True to get the randomized
//github. com Htmm/RemmLLUU/SSCCI
Remark For the Sobol soquence, we noticed that the implementation provided by Lelue et al. [2024] cannot be used in dimension higher than 20.
Riesz point configuration: We use a code provided by François Clement (https://sites. math. washington, edu/-fclement/ implementing projected gradient descent method, where choose the number of iterations T = 10, the gradient step as and 0.1. The function can be found in the riesz noblur py script in the repository https://githis wwwe//m///////eee /KennnSisouu Sii-Samlling--oi Guide.
Spherical Sliced Wasserstein: We used the following implementation from Bonet et al. [2023] that can be found in POT library (Python Optimal Transport) httpps:/pyythnno... gith...
18

---

Published in Transactions on Machine Learning Research (06/2025)
auto_ auto_exaaaaaa www.ppeesbaacknddd////eeaeeeeeee http:/pottsss_uiiitotrrro teres_tocchhhttt For the prper-parameters we set the num- ber of iteration T- = 250, the learning rate e 150 and the number of great circles T. = 500. For the initialization, we generate OM ~ sd-1 following the method described in Sec. 3.1.1.
Spherical Harmonics Control Variates: We use the implementation provided by Lelue et al. [2024], available in https://gi com/Rem1LLLC//SHOV They provide two possible functions SHCV and SW CV, both functions return value of SW estimate. These functions differ in the way they implement the optimization of Eq. 15. Depending on the number of control variates, one of the functions is much more efficient than the other. For this reason, in our experiments, we use both functions and always keep only the minimal error among the two.
4.2 Gaussian data
This part details the experiments on toy dataset chosen because it is simple to replicate and simple to understand. We compare different estimates of SW2 (pa,
(Wa( Pa) for d E (2, 5, 10, 20, 50}. We pick up Lelue
N
Σ Sy with 21, EN N2x X), y1..... UN ~ (y, Y), where N = 1000. The means are drawn as x, ~ W(1,, La) and the covariances are X, Y.Y = Eg DT Eyzy where all entries of the matrices are drawn using the standard normal distribution. In Fig. 4.1, we show convergence curves generated by all the different sampling strategies all the dimensions listed above. Fig. 4.2 reports the distance estimation error as a function of computation time (in seconds). In both figures, both axes are log scaled. We can see in Fig. 4.1 that up to dimension 5, Q.M.C. methods are preferable convergence wise, then the orthonormal sampling is preferable in dimension 20 and 50. In contrast, we can see in Fig. 4.2 that for dimensions less than 10, the S.H.C.V. method has a better error, with similar running time. For higher dimensions, however, the orthonormal sampling is much faster, for a given error target.
et al. [2024]'s setting, using pa = ト dr. and Va
Remark 9 : One may notice in Fig. 4.1b that both the S.H.C.V. method and the Q.M.C. method with the s- RRiesz points (R.R.) reach a plateau at around 103 projections. Our hypothesis is that both methods have a better estimation of SW3 than the simple random sampling with 10st projections that we use as a ground truth. We test this hypothesis in a simple case where SW2 (v. can be computed explicitly. We define a p= 160,, + )22] and v = _[0,, + &yul, with 31, 22 = (1, (0, -,, and // 1/2 (0, 0, 00,0, Simple computations yield (W2 2x-v2) Knowing the true value of SW2 (1, v),, we find that with 10 points, the Q.M.C. method with the s-Riesz points configuration and the S.H.C.V. methods already have errors order smaller that ones made by uniform sampling with 10 points.
3t
Remark 10: : Note that for the running time curves, we do not include the s- s-Riesz points configuration
starting from the dimension 3 because it takes around 102 seconds to generate 10 points and 9 x 10 seconds to generate 10a points. However, observe that those points, once generated, can be stored once for all to compute other SW2 distances any other Monte Carlo estimation problems for functions defined the unit sphere. This means that these configurations should not be discarded by default. For practical applications where the number of SW3 distances to compute is large, the computing time for these configurations can be factorized by the nunber of distances to compute and hence could become a negligible factor when the sampling number is moderate.
Remark 11 Recalling the rumning time complexity O(T in paragraph 3.2.2.2 and the rumning time results above, this shows that one needs to spend 9 x 107 seconds to generate 10 points. This demonstrates the limitation of this sampling method in terms of scalability, in other words when one needs a very large sampling number.
4.3 Persistence diagrams reduction dimension score
The goal of this section is to evaluate the relevance of the sampling methods studied in Sec. 3, in the context of concrete use case, involving two-dimensional real-life datascts. For that, we focus in this section on persistence diagrams, a popular object used in Topological Data Analysis [Edelsbrunner & Harer, 2009]. Persistence diagrams data abstractions encapsulating the features of interest of complex input datasets
19

---

Published in Transactions on Machine Learning Research (06/2025)
I.R.
S.H.C.V.
10
0000000
R.R.
O.S.
10
Error
10
H.R.A. M.
10
P.P.S.
P.R.P.S.
10
S.A.M.
S S.R.A.M.
102
10
10
10a
10
10
Number of samples
Number of samples
(a) 2D
(b) 3D
10
S.N.C.V.
S.I.C.V.
R.R.
R.R.
U.S.
U.S.
10-1
O.S.
O.S.
10
H.N.M.
H.N.M.
Error
S.R.N.M
10
S.R.N.M.
10
S.N.W.R.
S.S.W.R.
10
3
102
10
10
102
10a
104
Number of samples
Number samples
(c) 5D
(d) 10D
S.H.C.V.
10
U.S.
O.8.
10
O.S.
H.N.M.
H.N.M.
H.R.N.M.
Bevor
6-90
H.R.N.M.
10
S.S.W.R.
S.N.M.
10
S.R.N.M.
10-2
10
102
10
10
10
10a
104
Number of samples
Number of samples
(e) 20D
(f) 50D
Figure 4.1: Comparison of convergence rate results for the studied sampling methods (Gaussian data, Sec. 4.2).
(e.g- scalar fields) into simple two-dimensional representations. Specifically, we consider an input dataset represented as a piecewise linear (PL) scalar field, namely a function f : M + R defined on a PL (dm)-
manifold M with dm s 3. Take a value a R, we denote f=l(a) = f-1( oc, the sub-level set of f at a. While increasing a, the topology of (a) changes at the critical points f in M. Those critical points are classified by their index I: for minima, for 1-saddles, dı for (dM - 1))saddlles and dı for maxima. Following the Elder rule [Edelsbrunner Harer, 2009], a topological feature of f= (a) (connected component, cycle, void) is associated with pair of critical points (c, c') such that f(c) < f(c))
20

---

Published in Transactions on Machine Learning Research (06/2025)
I.R.
S.H.C.V.
10
0-000000--000000 0.0000000
U.S.
O.S.
10-1
Error
10
P.P.S.
10-2
P.R.P.S.
S.A.M.
S.R.A.M.
10
S.5.W.R.
100
10-2
10
10
10
10
10
10
10
10
Timer
Timer
(a)
(b) 3D
10
S.H.C.V.
S.H.C.V.
U.S.
10-1
10
M.R.N.M.
H.R.N.C
Error
S.N.M.
S.N.M.
S.R.N.M
S.R.N.M.
S.S.W.R.
10
S.S.W.R.
10
10
10
10
10-2
10
10
10
10a
10a
10
10
10
10a
Timer
Timer
(c) 5D
(d) 10D
S.H.C.V.
10
H.N.M.
10
H.N.M.
H.R.N.M
H.B.N.M.
S.N.M.
Eriu
S.N.M.
10-1
9-00
S.R.N.M
10
S.S.W.R.
10-2
10
111.11
m0I
111I
10
10
10"
10
10 3
10
10
10
10a
10a
Timer
Timer
(e)
20D
(I)
50D
Figure 4.2: Distance estimation error function of computation time (seconds). Computation times include the point generation as well the SW3 distance approximation.
and I. I 1. This pair corresponds to the birth and death of the topological feature during the sweep of the range from -oo to +oo by a, and it is called a persistence pair. As an example, when two connected components of / (a) merge at a critical point é', the younger one (created last) dies to let. the older one (created first) live on. Then those persistence pairs are represented as 2D points where the horizontal the birth of a topological feature (noted b f(c)) and where the vertical one corresponds its death (noted f(c)). The lifespan of feature is called persistence and is simply encoded as b d. This representation is called the Persistence Diagram, and its popularity in topological
coordinate corresponds
21

---

Published in Transactions on Machine Learning Research (06/2025)
data analysis is explained by its stability to the addition of noise. See Fig. 4.3 for simple example of persistence diagram.
Remark 12 Two persistence diagrams can have different number of points, so to make balanced
Ni
No
transport problem one has to augment them. Formally, denoting d1 Ni doz d2 No 6yz the
Σ
N1 Σ (n)())): Ad = No k=1
N1
diagrams, and noting Ad, = Ni k=l
Σ (na(m)) their projections on the diagonal A, one
considers N [N1d + N24d, and N [N2d2 + N14d, as input measures with N = N₁ + Then
the Sliced Wasserstein distance can be used to compare persistence diagrams as detailed by Carrière el al. [2017].
Figure 4.3: A simple example of persistence diagram issued from gaussian mixture (left). On the right you see that the persistence diagram is stable to the addition of noise.
We present dimensionality reduction results on ensembles of persistence diagrams [Pont et al., 2021] described in [Pont et al,, 2022], which original scalar fields include simulated and acquired 2D and 3D ensembles from SciVis constests [Organizers, 2004]. The dimensionality reduction techniques used are MDS [Kruskal & Wish, 1978] and t-SNE [van der Maaten & Hinton, 2008] applied distance matrices obtained by the SW estimations between the persistence diagrams. For given technique, one quantifies its ability to preserve the cluster structure of an ensemble by running the k-means algorithu in the resulting 2D- layouts. Then one evaluates the quality of the clustering with the normalized mutual information (NMI) and adjusted rand index (ARI), which should both equal to for a clustering that is identical to the classification ground-truth. Tab. 4 shows the average clustering scores of both MDS [Kruskal Wish, 1978] and (-SNE [van der Maaten & Hinton, 2008). First we take the average from distance matrices made by each SW3 estimates all sampling number {100, 1100, 2100, 3100, 4100, 5100, 6100, 7100, 8100, 9100, 10100). Then we average again over all the 12 different ensembles of persistence diagrams. One can see that all the methods are quite similar. But overall the s- s-Riesz points configuration, which are just the M-th unity roots up to a rotation, is slightly better.
Table 4: Average NMI and ARI scores for over all 12 ensembles of persistence diagrams.
Method
MDS NMI
t-SNE NMI
Method
MDS ARI
t-SNE ARI
Riesz
0.74
0.65
Riesz
0.64
0.51
Uniform
0.74
0.59
Uniform
0.64
0.44
Orthonormal
0.75
0.63
Orthonormal
0.64
0.48
Halton
0.74
0.58
Halton
0.63
0.41
4.4 3D Shapenet 55Core Data
This part details convergence results on a 3D datasel commonly used as a benchmark when sludying shape comparison techniques. So as in [Nguyen et al., 2024] and [Leluc et al., 2024], we took three 3D point clouds issued from the SharenctCore dataset introduced by [Chang al., 2015]. Among the different shapes in the dataset, we took one lamp, one plane and one bed; with all three of them having N 2048 points. Fig. 4.4 displays the three datasets considered for this experiment.
22

---

Published in Transactions on Machine Learning Research (06/2025)
Figure 4.4: The three point clouds taken from the ShapenetCore dataset (a plane, lamp and bed).
10
10
S.H.C.V.
10
0
10
10
10
I.R.A.M
F.P.8
F.I.P.S
S.B.W.R.
10
10-6
10-7
10a
10a
10a
10a
10a
10
10a
10a
10
Number samples
Number samples
Number of samples
Figure 4.5: Comparison of convergence rate results from the different sampling methods. The first plot shows errors made with respect to the SW? distance between lamp and a plane. The second one between a plane and bed. The last one corresponds SW? between plane and bed.
Fig. 4.5 shows different convergence curves of Sliced Wasserstein estimates between the three point clouds. As in Sec. 4.2, the methods dominating are the Q.M.C., R.Q.M.C., S.S.W. and S.H.C.V. methods, especially the s-Riesz points configuration and the Spherical Sliced Wasserstein sampling.
4.5 MNIST reduction dimension score
The goal of this section is twofold. First, it evaluates the practical convergence of the studied sampling methods on real-life high-dimensional datasets. Second, it describes an application of the SW distance for high-ddimensional data, namely, dimensionality reduction. For this, we select the classical MNIST dataset [LeCn, 1998]. To construct our dataset, we represent each digit image as a point in 228x28 For each class {0, 2, 3, 4,5, 5, 6, 7,8, 9), select randomly 600 digit images and divide them into groups of 200. This results in 30 point clouds of 200 points each, in R28x28, with 10 ground-truth classes. Fig. 4.6 illustrates the 30 x 30 matrix of SW distances between all point clouds in the dataset. We use MDS and t-SNE to produce 2D layouts from the distance matrices generated by the various sampling methods with different sample sizes. We then apply a clustering algorithm to these 2D layouts and average the clustering scores (NMI and ARI, see Sec. 4.3) on all sampling numbers for all the studied sampling strategies. Results are provided in Tab. 5. In such high dimension (d 784), we that the performance of L.D.S. collapse, the three sampling methods standing out being the s- Riesz points configuration, the uniform sampling and the orthonormal sampling.
5 Recommendation & conclusion
In this paper, we have studied several sampling strategies on the sphere for computing an approximation of the Sliced Wasserstein distance.
23

---

Published in Transactions on Machine Learning Research (06/2025)
Figure 4.6: Sliced Wasserstein distance matrix of our dataset using 10 projections. All 10 classes, 7, of 3 members each are well represented in the matrix.
Table 5: Average NMI and ARI scores with standard deviation. Higher scores correspond to better clustering.
Method
MDS NMI
t-SNE NMI
Method
MDS ARI
t-SNE ARI
Riesz
1 0.
0.98 1 2c-2
Riesz
1 + 0.
0.95 1 7e-2
Uniform
1 1 0.
0.97 + 4e-2
Uniform
1 + 0.
0.91 le-1
Orthonormal
1 1 0.
0.98 i 3e-2
Orthonormal
1 = 0.
0.94 + 8e-2
Halton
0.91
le-1
0.91 9e-2
Halton
0.75 = 2
0.76 2e-1
S.S.W.
1 + 0.
0.98 + 4e-2
S.S.W.
1 = 0.
0.94 + le-1
Regarding theoretical guarantees, this study highlighted the following limitations. The classical i.ddd sam- pling benefits from theoretical guarantees with convergence rate in O(1/VM) and a time complexity linear in the number M projections. Orthonormal sampling and L.D.S. such as Halton or Sobol lack conver- gence rate guarantecs on the sphere (these guarantees being only obtained for sequences on hypercubes for L.D.S). As for deterministic point generation methods (like Riesz), the Sliced Wasserstein integrand also lacks sufficient regularity guarantee results in dimensions higher than
While lacking theoretical guarantees in terms of convergence, the experimental study suggests that Q.M.C methods D...D... or s-Riesx points) provide competitive results in small to intermediate dimension, while having a similar convergence rate to classical random sampling methods in intermediate to higher (for Riesz) dimensions. These results scom to indicate that, while f is not rogular enough for the convergence guaran- tees detailed in this paper, there may some non-proven convergence results requiring weaker regularity conditions that would be applicable SW.
Now, considering computation limes, as shown by Fig. 4.2 and Tab. 2, classical i.i.d. sampling remains the slowest method in all our experiments. While orthonormal sampling lacks theoretical guarantees, it seems to be one of the most efficient methods whatever the dimension, and is particularly competitive in high dimensions, with very reasonable increase of computation time. L.D.S. methods also remain competitive in pratice for small dimensions. s-Riesa points, while competitive in terms convergence rate, have a prohibitive time complexity O(M3), which makes them completely unsuitable for large number of projections.
The experiments also suggest that the S.H.C.V. method is very competitive in intermediate dimensions, while becoming less efficient when d increases.
Based on the different experimental results provided in this paper, we make the following recommendations:
For small dimensions (less than 3). Q.M.C. methods such as s-Riesz points or L.D.S. mapped onto the sphere be privileged with respect to uniform sampling,
24

---

Published in Transactions on Machine Learning Research (06/2025)
For high dimensions (greater than 20), the orthonormal sampling method emerges the most suitable choice. It is also one of the simplest methods to implement, which makes it particularly attractive in practice.
For intermediate dimensions (between 5 and 10), choosing an appropriate method should depend on the experimental requirements. Spherical harmonics are an excellent option if computational resources are limited and if the number of SW distances to be computed is low. However, it is worth noting that some Q.M.C. strategies, being independent of the input measures, have the advantage of allowing the generaled points to be reused and of allowing an independent computation in M (except the Riesz points). This should be particularly beneficial when high number of projections is roquired and large number of SW distances must be computed. In such cases, suggest to store the samples to factorize the computing time across experiments.
Acknowledgments
This work is partially supported by the European Commission grant ERC-2019-COG "TORI" (ref. 863464, https: //erc-tori, seno--ttr.ggtub..... github.io/). JD acknowledge the support of the "France 2030" funding ANR-23-PEIA-000: ("PDE- AI") and the support of the funding SOCOT- ANNR22--E40-0017.
References
C. Aistleitner, J. S. Brauchart, and J. Dick. Point sets on the sphere S² with small spherical cap discrepancy. Discrete and Computational Geometry, August 2012. ISSN 1432-0444. doi: 10.1077/40454001294511.. URL /ttp://fo....aoo/0000.0000000000
Christoph Aistleitner, Florian Pausinger, Aune Marie Svane, and Robert F, Tichy, On functions of bounded variation, 2016.
George B. Arfken, Hans J. Weber, and Frank E. Harris. Mathematical Methods for Physicists: A Comprchensive Guide. Academic Press, 2011. ISBN 9780123846356. URL https://books google.
Stt/bookaPiddJJjHH.J-----
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. Tn International conference on machine learning, pp. 214 223. PMLR, 2017.
S. Asmussen and P.W. Glynn. Stochastic simulation: Algorithms and analysis. Springer Verlag, 2007.
Kinjal Basu. Quasi-Monte Carlo Methods in Non-Cubical Spaces. Phd thesis, Stanford University, 2016.
Clément Bonet, Lucas Drumetz, and Nicolas Courty. Sliced-wwsserstein distances and flows on cartan- hadamard manifolds. arXiv proprint AXXiv:2403.06600 2024.
Clément Bonet, Paul Berg, Nicolas Courty, François Septier, Lucas Drumetz, and Minh- Tan Pham. Spherical sliced-wassersteinn, 2023. URL attps://arrae www.oogg//////eeeoooo/ooooeee 08780.
Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and Radon Wasserstein baryeen- (ers of measures. Journal of Mathematical Imaging and Vision, 51()):22 45, 2015.
Nicolas Bonnotte. Unidimensional and evolution methods for optimal transportation. PhD thesis, 2013.
Luca Brandolini, Lconardo Colzani, Giacomo Gigante, and Giancarlo Travaglini. On the inequality. Journal of Complexity, 29(2):158 172, 2013.
J. S. Brauchart. Optimal discrete riesz energy and discrepancy, 2011.
25

---

Published in Transactions on Machine Learning Research (06/2025)
Johann S. Brauchart and Josef Dick. A characterization of sobolev spaces on the sphere and an extension of stolarsky's 's invariance principle to arbitrary smoothness. Constructive Approximation, 38(3):397 444, October 2013. ISSN 1432-0940. doi:
URL http://dx.dee doi. www.tooo...o.oooooo 1007/
600365-013-9217-2
Mathieu Carrière, Marco Cuturi, and Steve Oudot. Sliced wasserstein kernel for persistence diagrams, 2017.
URL https://arxiv. http:///xxxvvoogggaee/ http/abs/17000 033588.
Angel X. Chang, Thomas Funkhouser, Leonidas Guibas, Pat Hanrahan, Qixing Huang, Zimo Li, Silvio Savarese, Manolis Savva, Shuran Song, Hao Su, Jianxiong Xiao, Li Yi, and Fisher Yu. Shapenet: An information-rich 3d model repository, 2015. URL https://arxiv: wttp:///////....aaaa// 10bpp//111........
Nicolas Courty, Rémi Flamary, Devis Tuia, and Alain Rakotomamonjy. Optimal transport for domain adaptation. TEEE transactions on pattern analysis and machine intelligence, 59(9):1853-180)) 2016.
Julie Delon, Julien Rabin, and Yann Gousseau. Transportation distances on the circle and applications, 2010. URL https://arxiv www.l/voorg/////eeeeeeeee 5499.
Ishan Deshpande, Ziyu Zhang, and Alexander G Schwing, Generative modeling using the sliced wasserstein distance. In Procoedings the IEEE conference on computer vision and pattern recognition, pp. 3483- 3491, 2018.
Josef Dick and Friedrich Pillichhaammer, Digital nets and soquences: discrepancy theory and quasi-Mons Carlo integration. Cambridge University Press, 2010.
H. Edelsbrunner and J. Harer. Computational Topology: An Introduction. American Mathematical Society, 2009.
Ariel Elnekave and Yair Weiss. Generating natural images with direct patch distributions matching. 03 2022.
Renjie Feng, Friedrich Götze, and Dong Yao. Determinantal point processes on spheres: multivariate linear statistics, 2023. URL
Jean Feydy, Benjamin Charlier, Saançois-Xavier Vialard, and Gabriel Peyré. Optimal transport for diffeo- morphic rogistration. Medical Image Computing and Computer Assisted Intervention- MICCAI 2017: 20th International Conference, Quebec City, QC, Canada, September 11-13, 2017, Proceedings, Part 20, pp. 291 299. Springer, 2017.
H. Fischer. A History of the Central Limit. Theorem: From Classical to Modern Probability Theory. Sources and Studies in the History of Mathematics and Physical Sciences. Springer New York, 2010.
ISBN 9780387878577. URL https://books. betooks.googge www.goo...boooooeeoooooooooo wtt.://ooooss.dd//////aaaaaia..
P.T. Fletcher, Lu Conglin, S.M. Pizer, and Joshi Sarang. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Transactions on Medical Imaging, 23(8):995 1005, 2004. doi: 10.1109/TMI.200 2004. 831793.
Ishaan Guirajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville. Improved training of wasserstein gans. Advances in neural information processing systems, 30, 2017.
M. Götz. On the riesz energy of measures. Journal of Approximation Theory, 122(1):62--8, 2003. ISSN 0021-9045. doi: http://ceceeaarrcccc///////20000
URL https: httpes//www.siiencetiree.. sci.nceedeee.......
John H Halton. Algorithm 247: Radical-inverse quasi-random point sequence. Communications of the ACM, 7(12):701-702, 1964.
J.M. Hammersley and D.C. Handscomb. Monte Carlo Methods. Methuen's monographs on applied proba- bility and statistics. Methuen, 1964. ISBN 9780416523409. URL https: ttttps//bbokke /btoks/oooks///e//o/ google fr/books7id KK4OAAAAQAAJ.
26

---

Published in Transactions on Machine Learning Research (06/2025)
D. P. Hardin, T. J. Michacls, and E. B. Saff. A comparison of popular point configurations on S?, 2016.
Emnanuel Hebey. Sobolev Spaces on Riemaniann Manifolds. Springer, 1996.
Eric Heitz, Kenneth Vanhoey, Thomas Chambon, and Laurent Belcour. A sliced wasserstein loss for neu- ral texture synthesis. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp. 9412 9420, 2021.
Johannes Hertrich, Antoine Houdard, and Claudia Redenbach. Wasserstein patch prior for image superres- olution. IEEE Transactions on Computational Imaging, 8:693 704, 2022.
Sungkyu Jung. Geodesic projection of the mises misos-fisher distribution for projection pursuit of directional data. Electronic Journal of Statistics, 15, 01 2021. doi: 10.1214/21-EJS1811
Sungkyu Jung, Ian L. Dryden, and J. S. Marron. Analysis of principal nested spheres. Biometrika, 99 (3):551-566 568, 07 2012. ISSN 0006-3444. doi: 10.1093/bioni /as//m///////222...... URL https://dol http:////oo.org//...1o 1099//.1993/ biomet/ass022.
Soheil Kolouri, Phillip E Pope, Charles E Martin, and Gustavo K. Rohde. Sliced wassenstein auto-encoders. In International Conference on Learning Representations, 2018.
B Kruskal and M Wish. Multidimensional Scaling. In SUPS, 1978.
Lauwerens Kuipers and Harald Niederreiter. Uniform distribution of sequences. Courier Corporation, 2012.
Tung Le, Khai Nguyen, Shanlin Sun, Nhat Ho, and Xiaohui Xic. Integrating efficient optimal transport and functional maps for unsupervised ised shape correspondence learning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 23188 23198, 2024.
Yann LeCun. The mnist database of handwritten digits. http://yann. lecun. com/exdb/mnist// 1998.
Chen-Ye Yonn-Yulo Lee, Tanmay Batra, Mohammad Haris Baig, and Daniel Ulbricht. Sliced wasserstein discrepancy for unsupervisod domain adaptation. In Procoedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 10285-10295, 2019.
Rémi Lelue, Aymeric Dieuleveut, François Portier, Johan Segers, and Aigerim Zhuman. Sliced-wasssrsteiin estimation with spherical harmonics as control variates, 2024. URL actps://arxivv srtrvvv.orr/abe//420.0044.
Christiane Lemieux. Monte carlo and quasi-monte carlo sampling. 2009.
Xinran Liu, Yikun Bai, Rocío Diaz Martin, Kaiwen Shi, Ashkan Shahbazi, Bennett A. Landman, Catie Chang, and Soheil Kolouri. Linear spherical sliced optimal transport: A fast metric for comparing spherical data, 2024. URL https: http://axxvvvvv//e wwtp:os//////////ge/ ab4/ 06055.
Ricardo Marques. Bayesian and Quasi-Monte Carlo spherical integration for global illumination. Phd thesis, Université de Rennes, 2013.
Jordi Marzo and Albert Mas. Discrepancy of minimal riess energy points. Constructive Approximation, 54 (3):473 506, 2021.
Claus Müller. Introduction. Springer New York, New York, NY, 1998. ISBN 978-1-4612-0588-44 doi: URL http:::/ddoi.org//1//
11.-007977---46120555
Kimia Nadjahi. woreedd-sssrsttein distance for largo-scale machine learning: theory, methodology and extensions. PhD thesis, Institut polytechnique de Paris, 2021.
Khai Nguyen and Nhat Ho. Control variate sliced wasserstein istioos 05 2023.
Khai Nguyen and Nhat Ho. Sliced wasserstein estimation with control variates, 2024. URL https:
arg/abs/2300 00402.
27

---

Published in Transactions on Machine Learning Research (06/2025)
Khai Nguyen, Nicola Bariletto, and Nhat Ho. Quasi-monte carlo for 3d sliced wasserstein, 2024.
Harald Niederreiter. Low-discrepancy and low-dispersion sequences. Journal of number theory, 30(1):51 70,
1988.
Organizers. The IEEE SciVis Contest. http://ssclviscontes.. wttp://iwwioottteeteeeeees....s
2004.
Art B Owen. Multidimensional variation for quasi-monte carlo. In Contemporary Multivariate Analysis And Design Of Experiments: Celebration of Professor Kai-Tai Fang's 65th Birthday, pp. 49 74. World Scientific, 2005.
Art B Owen. Monte carlo book: the quasi-monte carlo parts, 2019.
Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends in Machine Learning. 11(5-6):355 607, 2019.
Mathieu Pont, Jules Vidal, Julie Delon, and Julien Tierny. Wasserstein Distances, Geodesics and Barycenters of Merge Trees Ensemble Benchmark. 2021.
Mathieu Pont, Jules Vidal, Julie Delon, and Julien Tierny. Wasserstein Distances, Geodesics and Barycenters of Merge Trees. JEEE TVCG, 2022.
Michael Quellmalz, Robert Beinert, and Gabriele Steidl. Sliced optimal transport the sphere. Inverse Problems, 1550010 August 2023. ISSN 1361-6420. doi: 10.1088/1361 (1000//1361--12/accffaaa URL http: 7/dx. doi. 20g/10.10 11111/366161010
Julien Rabin, Julie Delon, and Yann Gousseau. A statistical approach to the matching of local features. SIAM Journal on Imaging Sciences, ((3))93,-958, 2009.
Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot. Wasserstein barycenter and its application to texture mixing. In Scale Space and Variational Methods in Computer Vision: Third International Conference, SSVM 2011, Ein-Gedi, Israel, May 29 June 2, 2011, Revised Selected Papers 3, pp. 435 3454466. Springer, 2012.
Mark Rowland, Jiri Hron, Yunhao Tang, Krzysztof Choromanski, Tamas Sarlos, and Adrian Weller. Or- thogonal estimation of wasserstein distances, 2019.
Tim Salimans, Han Zhang, Alec Radford, and Dimitris Metaxas. Improving gans using optimal transport. In International Conference Learning Representations, 2018.
Filippo Santambrogio. Optimal transport for applied mathermaticians. Birkäuser, NY, 55(58-63):94, 2015. Sobol. The distribution of points in a cube and the accurate evaluation of integrals (in russian) zh. Vychisl. Mat. Mater. Phys, 7:784-802. 1967.
Kenneth B. Stolarsky. Sums of distances between points on a sphere. ii Proceedings of the American Mathematical Society, (1(2):575555 582, 1973. ISSN 00029939, 10886826.
Eloi Tanguy, Rémi Flamary, and Julie Delon. Reconstructing discrete measures from projections. conse- quences on the empirical sliced wasserstein distance. Comptes Rendus Mathématique, 2023.
J. G. van der Corput. Verteilngssfunktionen.. I. Proc. Akad, Wet. Amsterdam, 38:813-82 1935. ISSN 0370-0348.
L.J. P. van der Maaten and G.E. Hinton. Visualizing Data Using t-SNE. JMLR, 2008.
C. Villani. Optimal transport: old and new. Springer Verlag, 2008.
Jiqing Wu, Zhiwu Huang, Dinesh Acharya, Wen Li, Janine Thoma, Danda Pani Paudel, and Luc Van Gool. Sliced wasserstein generative models. In Proceodings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 3713 3722, 2019.
28

---

Published in Transactions on Machine Learning Research (06/2025)
Xu Xianliang and Huang Zhongyi. Central limit theorem for the sliced 1-wasserstein distance and the max-sliced 1-wasserstein distance, 2022.
29


Reviews:
[{"review_id": "icMjAML5rH", "reviewer": "Reviewer_u77K", "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": "Strengths:\n\n- The work is very thorough, covering both theoretical aspects and practical implementations.  \n- A clear view of sampling methods is provided. This will help readers understand trade-offs (e.g., between convergence rate and computational overhead).  \n- The experimental section is broad, spanning synthetic Gaussian mixtures, persistence diagrams, 3D point clouds, and MNIST.  \n\nWeaknesses: \n- Some parts of the paper are quite dense and technical, which may pose challenges for practitioners who are less familiar with advanced numerical integration and optimal transport theory.\n- Some sections (especially the theoretical derivations) could benefit from additional intuition or summary to guide readers through the complex analysis.  \n- The experimental comparisons, while broad, could be more detailed in discussing the limitations or failure cases of some methods, especially in very high-dimensional settings.  \n- Implementation details (e.g., parameter choices for optimization routines or sensitivity to hyperparameters) are mentioned but might require deeper discussion for reproducibility.", "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 paper provides a comprehensive survey of various sampling strategies used in the estimation of the Sliced Wasserstein distance. It revisits key theoretical properties of the SW distance—such as its continuity, piecewise smoothness, and Lipschitz behavior—and uses these properties to motivate a taxonomy of sampling methods. In particular, the authors detail three broad classes of approaches: random samplings, discrepancy-based approaches, and variance reduction techniques. Some numerical experiments are shown to accompany the methodologies presented.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The authors should add more intuitive explanations alongside the technical derivations—especially regarding the behavior of the discrepancy measures and the curse of dimensionality—to make the material more accessible.\n- It would also be beneficial to include a more in-depth discussion of the limitations observed in the experiments, such as scalability issues and the effect of hyperparameter choices on convergence.\n- The authors should provide a more detailed description of the experimental setup (e.g., parameter settings and implementation choices) or direct readers to a well-documented repository to ensure reproducibility.\n- A dedicated section discussing the trade-offs between theoretical guarantees and actual computational overhead in practice.", "broader_impact_concerns": "Not applicable.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "NEGjqnvNpr", "reviewer": "Reviewer_BLhD", "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": "Strengths:\n\n- The comparative study is thorough where most existing methods for sampling from a sphere are presented and experimented.\n- The theoretical footings of each method are presented and discussed.\n\nWeaknesses:\n\n- This work is purely a comparative study and does not introduce any technical contributions.\n- The insights are largely empirical.", "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 conducts a comparative study of various sampling methods from the unit sphere for computing the sliced Wasserstein distance between two distributions. The authors examine three main classes of sampling strategies: (i) random sampling (uniform and orthonormal), (ii) discrepancy-based sampling (e.g., Riesz points, Sobol, Halton), and (iii) variance reduction techniques (e.g., spherical harmonics control variate).\n\nThe main contribution of the paper is the comprehensive comparison between these methods, providing practical guidelines on selecting the most suitable strategy based on the problem's dimensionality and computational constraints.\n\nAdditionally, the paper establishes the regularity of the (1D) Wasserstein distance between the projections of two discrete measures along a given direction $\\theta$.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- In section 2.2., it lacks discussions on how to prove proposition 1 when $\\mu$ and $\\nu$ are continuous.\n- Right after Proposition 3, it is unclear to me if $Y_M$ is an unbiased/consistent estimator for the sliced Wasserstein distance. It said \"As a consequence, this method inherits the properties of classical Monte Carlo, with a rate of convergence independent of the dimension of the input measures\", referring to the variance inequality. However, it did not say anything about the bias, so it is not enough to conclude.\n- Three methods in Section 3.2.1.4. would need proper references.\n- In the experiments, it is unclear to me why some methods are missing from certain dimensions. For example, in Figure 4.1, the legends of the subplots are not at all consistent.\n- The Koksma-Hlawka inequality is referenced a couple times in the discussion, yet no formal definition is provided.\n\nTypo:\n\n- In the caption of Figure 2.1 (b), $SW_2^2(_m u,\\nu)$ should be $SW_2^2(\\mu,\\nu)$ .", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "dqTcEF3lW4", "reviewer": "Reviewer_b6ih", "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": "This work provides a nice overview of different strategies to approximate the Sliced-Wasserstein distance with Monte-Carlo estimators. In addition to compare known methods (such as uniform sampling, orthonormal sampling, and some of the QMC methods), it also compares other methods which have not been tested to compute SW to the best of my knowledge (e.g. Fibonacci point sets, or minimizers of Riesz energy).\n\nI believe this is a valuable contribution to the community, as it provides a useful summary for practitioners on selecting an estimate Sliced-Wasserstein distance.\n\n**Strengths:**\n- Analysis of the regularity of the integrand\n- Clear description of different strategies to approximate integrals with respect to the uniform measure on the sphere, with known results on the theoretical side\n- Comparison of the methods in several settings\n- Guidelines on which method to choose given the setting\n\n\n**Weaknesses:**\n- As acknowledged in the paper, some methods lack theoretical guarantees on the sphere\n- One might say that MNIST is not a very high dimensional setting, and it could be interesting to see what happens on an even higher dimensional dataset.", "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 studies the approximation of the Sliced-Wasserstein distance by comparing different approximations of the integral with respect to the uniform measure on the sphere. It begins by providing a theoretical analysis of the regularity of the integrand of the Sliced-Wasserstein distance. Then, it describes different strategies to estimate the integral, and provides complexity and known theoretical results for the different strategies. Finally, the authors provide a comparison of the described strategies in different settings, and make recommendations on which method to use given the setting.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As there are a lot of results, it would be nice in addition to Table 1 to also summarize in a Table the known results of the different method (i.e. the rates, the computational complexities...).\n\nFor the timings reported in Figure 4.2, I am wondering if there are all optimal. I had a quick look at the code, and it seems that the uniform sampling normalization is done with a loop, while it could be done in parallel (which might speed things up a bit?). I also think that the orthonormal sampling could be done in parallel by first sampling a $(M\\times d\\times d)$ matrix with independent Gaussian components, and then taking the QR decomposition, see e.g. [1] or [2]. I don't know if it would do a big difference though.\n\nI am also wondering if sampling uniform points on the sphere by e.g., minimizing the spherical Sliced-Wasserstein distance using the closed-form in Proposition 1 of [2], could provide good results?\n\n\nTypos:\n- Page 2: \"Indeed the higher the dimension the harder it is to generate samplings with Q.M.C. and R.Q.M.C.\napproaches.\" -> to generated samples?\n- Equation 11: \"$\\mathbb{E}_{\\theta\\sim \\mathcal{U}(\\mathbb{S}^{d-1})]}$\"\n- In the proof of Proposition 1, it is written \"On each $V_l$, $f_{|\\mathbb{S}^{d-1}}=f_{\\sigma_l \\mathbb{S}^{d-1}}$. I guess it should be $=f_{\\sigma_l}$?\n- Table 1: \"Incremental computation\": I am not sure to understand what it means. I would say that incremental computations are when we can't compute points in parallel.\n- Section 3.1.2 title: \"Orthonornal\"\n- Page 13: \"s-Riesz energy for $s=-1$\", but is only defined for $s\\ge 0$?\n- 2 citations of the same paper \"Sliced Wasserstein Auto-encoders\" (with publication dates 2018 and 2019)\n\n\n[1] Lin, T., Zheng, Z., Chen, E., Cuturi, M., & Jordan, M. I. (2021, March). On projection robust optimal transport: Sample complexity and model misspecification. In International Conference on Artificial Intelligence and Statistics (pp. 262-270). PMLR.\n\n[2] Bonet, C., Berg, P., Courty, N., Septier, F., Drumetz, L., & Pham, M. T. (2022). Spherical sliced-wasserstein. arXiv preprint arXiv:2206.08780.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents a comprehensive survey and comparative analysis of sampling strategies for approximating the Sliced Wasserstein (SW) distance. It begins with a theoretical review of the SW distance’s continuity, smoothness, and Lipschitz properties, which motivate a taxonomy of estimation methods. The authors categorize sampling strategies into three main classes: (i) random sampling (uniform, orthonormal), (ii) discrepancy-based methods (e.g., Sobol, Halton, Riesz), and (iii) variance reduction techniques (e.g., spherical harmonics control variate). The paper includes complexity analysis, theoretical insights, and practical guidelines based on numerical comparisons across different settings and dimensionalities.

Overall, this paper provides a valuable summary of important sampling techniques for the Sliced Wasserstein distance and is a meaningful contribution to the machine learning community. All reviewers agree on accepting the paper, and their concerns have been adequately addressed. Therefore, I also recommend acceptance.
Decision: Accept (poster)