Gromov-Wasserstein-like Distances in the Gaussian Mixture Models Space
Authors: Antoine Salmona
Agnes Desolneux
Julie Delon
Venue: Accepted by TMLR ()

Abstract: The Gromov-Wasserstein (GW) distance is frequently used in machine learning to compare distributions across distinct metric spaces. Despite its utility, it remains computationally intensive, especially for large-scale problems. Recently, a novel Wasserstein distance specifically tailored for Gaussian mixture models and known as  $ MW_2 $ (mixture Wasserstein) has been introduced by several authors. In scenarios where data exhibit clustering, this approach simplifies to a small-scale discrete optimal transport problem, which complexity depends solely on the number of Gaussian components in the GMMs. This paper aims to extend $ MW_2 $  by introducing new Gromov-type distances. These distances are designed to be isometry-invariant in Euclidean spaces and are applicable for comparing GMMs across different dimensional spaces. Our first contribution is the Mixture Gromov Wasserstein distance ($MGW_2$), which can be viewed as a ’Gromovized’ version of $ MW_2 $ . This new distance has a straightforward discrete formulation, making it highly efficient for estimating distances between GMMs in practical applications. To facilitate the derivation of a transport plan between GMMs, we present a second distance, the Embedded Wasserstein distance ($ EW_2 $). This distance turns out to be closely related to several recent alternatives to Gromov-Wasserstein. We show that  can be adapted to derive a distance as well as optimal transportation plans between GMMs. We demonstrate the efficiency of these newly proposed distances on medium to large-scale problems, including shape matching and hyperspectral image color transfer.

Published in Transactions on Machine Learning Research (09/2024)
Gromov-Wassersteen-iikk Distances in the Gaussian Mixture Models Space
Antoine Salmona
ENS Poris Saclay, CNRS, Centre Borelli UMR 9010
Julie Delon
Université de Paris, CNRS, MAP5 UMR 8145 and Institut Universitarre de Fronce
Agnès Desolneux
ENS Poris Saclay, CNRS, Centre Borelli UMR 9010
Reviewed on OpenReview: https: opwneeveee.eeeeeeoooeeee net/ forum? id- JT4JJT
Abstract
The Gromov-Wassorster (GW) distance is frequently used in machine learning to compare distributions across distinct metric spaces. Despite its utility, it remains computationally intensive, especially for large-scale problems. Recently, novel Wasserstein distance specif- ically tailored for Gaussian mixture models (written GMMs in the paper for the sake of brevity) and known as MW2 (mixture Wasserstein) has been introduced by several authors. In scenarios where data exhibit clustering, this approach simplifies to a small-scale discrete optimal transport problem, which complexity depends solely on the number of Gaussian components in the GMMs. This paper aims to incorporate invariance properties into MW2. This is done by introducing new Gromov-type distances, designed to be isometyy-invariant in Euclidean spaces and applicable for comparing GMMMs across different dimensional spaces. Our first contribution is the Mixture Gromov Wasserstein distance (MGW2), which can be viewed as a 'Gromovized' version of MW2. This new distance has a straightforward discrete formulation, making it highly efficient for estimating distances between GMMs in practical applications. To facilitate the derivation of a transport plan between GMMs, we present second distance, the Embedded Wasserstein distance (FW2). This distance turns out to be closely related to several recent alternatives to Gromov- Wasserstein. We show that EW2 can be adapted to derive distance well as optimal transportation plans between GMMs.. We demonstrate the efficiency of these newly proposed distances on medium to large-scale problems, including shape matching and hyperspectral image color transfer.
1 Introduction
The goal of optimal transport (OT) theory is to design meaningful ways to compare probability distributions. It provides very useful mathematical tools for diverse imaging sciences and machine learning lasks including generative modeling (Arjovsky et al., 2017; Genevay et al., 2018; Tolstikhin al., 2018), domain adaptation (Courty et al., 2016), image processing (Rabin et al., 2012; 2014), and embedding learning (Courty al., 2018: Xu et al., 2018). For two probability distributions i and v, respectively two Polish (i.e complete, separable, metrizable) spaces X' and }, and given lower semi-continuous function e: X x ) - R called cost, optimal transport in its most classic form aims at solving the following optimization problem,
inf / c(x, (x.
(1)
xell(w.))

---

Published in Transactions on Machine Learning Research (09/2024)
where I(,, is the set of probability measures on X' x )/ with marginals j and v. When ) is equal to X, the choice of cost cp(x; v) = dx(a, v)"; with p21 1 and dx the metric of the space X, induces distance between probability distributions with finite p-th moments, called the Wasserstein distance Wp. In the discrete setting, Problem (1) becomes
inf EO.wo. b))l( R,
where a = (a1,.... am)" and b = (bi...., 0n)) are respoctively in the R" and R" simplexes Am and An. ITa, b) = {w E Amxn : wl,, = a and wº' = bj and C is a non-negalive matrix of size m x n, called cost. matrix.
Optimal transport is known to be computationally challenging. Between discrete distributions, its computation involves solving linear program that rapidly becomes costly as soon as the number of points is moderately large. Between two sets of n points, its computation complexity is O(r"log(n)) (Seguy et al., 2017), which compromises its usability for settings with more than a few tens of thousand of points. To lighten OT computational cost, a large number of works have developped efficient computational tools. In particular, Cuturi (2013) proposes to solve an entropic regularized OT problem using the Sinkhorn-Knopp algorithm (Sinkhorn and Knopp, 1967), reducing the cost of the problem to O(n²). Over the last past years, a large body of works have focused speeding up the Saalhhorn-Koopp algorithm, building mostly on diverse low-rank approximations (Solomon et al., 2015; Altschuler et al., 2018; 2019; Forrow et al., 2019; Scetbon and Cuturi, 2020; Scetbon et al., 2021). These approaches have helped to reduce the computational cost of the problem from cubic (for the non-regularized problem) to linear complexity. Another lype of commonly used solvers are building on sliced mechanisms (Rabin et al., 2012; Kolouri et al., 2019). These solvers average Wasserstein distances between several one dimensional projections of the high-dimensional distributions, leveraging the fact. that the OT problem between onc-dimentional distributions can be solved using simple sorting algorithm. Alternatively, Delon and Desolneux (2020) have proposed an OT distance between Gaussion minture models (GMM), called Mixture Wasserstein (MW), where the admissible couplings are themselves constrained to be GMMs. They demonstrated that this specific continuous OT problem could be equivalently reformulated into a discrete version (which had been also proposed independently by Chen et al. (2018)): for two GMM with respectively Ko and K1 components, solving this formulation boils down to solve a small scale Ko x K₃ discrete OT problem. This distance can be applied to real data by first fitting GMMs on each distribution, making it particularly suited for scenarios where clustering structure already exists in the data. The main advantage of this approach is that its computational cost arises almost exclusively from fitting the GMMs to the
the data, since the complexity of the composite OT problem depends neither the dimension number of points, but solely on the number of components in the GMMs. This approach offers scalable and computationally efficient OT distance, which has been used for instance for texture synthesis (Leclaire et al., 2023), evaluating generative models (Luzi et al., 2023), Gaussian Mixture reduction (Zhang and Chen, 2020) or approximate Bayesian computation (Forbes et al., 2021).
One weakness of the classical optimal transport approach lies in the fact that it implicitly assumes that the spaces X' and y are comparable, i.c. that there exists relevant cost function c: X' - R. to compare them. Yet, this assumption is not always verified. For instance, if X = Rd and y Rdd with d * d', the definition of meaningful cost function c: Rdd x Rdd + R__ is not straightforward. Furthermore, some applications such as shape matching require having an OT distance that is invariant to a given family of transformations, such as translations or rotations, or more generally to isometries?. Even if the two distributions involved in these applications do live in the same ground space, it is not straight.foraat to design a cost. function such that the resulting OT distance will be invariant to these familics of transformations. To overcome those limitations, several non-convex variants of Problem (1) have been proposed (Cohen and Guibasm, 1999; Pele and Taskar, 2013; Alvarez-Melis et al., 2019; Cai and Lim, 2022). Among these, the Maanov-Wassesstein (GW) (Mémoli, 2011) distance is perhaps the most frequently utilized, recently gaining significant attention for the versatility it provides. Indeed, it only requires modeling topological aspects of the distributions within each domain to compare them without having to specify first a subset of invariances nor to design a relevant. cost function between the spaces the distributions lie on. The GW problem between two measures μ and v
The simplex Am is the subset of mm of (#1..... Fmy YT such that for all S S m, 24 2 0, and En 34 =1 We say that de X. + y is an isometry if for all (a) E X2, (((()) = (x)
2

---

Published in Transactions on Machine Learning Research (09/2024)
living respectively on X' and y aims at solving
inf Jexo xEll(
Lxyx xxxy
where ex : X' x X' - R and cy 3/x x y' - lK are two cost functions. Since this optimization problem only requires to define cost functions in each respective space, it remains very versatile and can be defined with very little assumptions on the spaces X' and y This approach has been applied to shape matching (Mémoli, 2009), or more generally to correspondence problems (Solomon et al., 2016), word embedding (Alvarez-Ma) and Jaakkola, 2018); graph classification (Vayer et al., 2019a), graph prediction (Brogat-Motte et al., 2022), and generative modeling (Bunne et al., 2019).
Computationally speaking, the Gromov- Wasserstein problem is known to be much more costly to solve than the classic linear OT problem. Indeed, the problem is non convex, quadratic with respect to it and known to be NP-hard. One possible approach to solve GW consists in linearizing the cost and to solve iteratively several classic OT problems. Entropic regularization of GW has also been proposed in (Peyré et al., 2016; Solomon et al., 2016) and results in a still non convex problem which can be solved by a projected gradient algorithm, where cach projection is itsolf an entropic lincar optimal transport problem. In recent years, several practical approximations of GW have been proposed in the literalure to reduce its compulational complexity and solve it efficiently, either through quantization of input measures (Chowdhury et al., 2021), recursive clustering approches (Xu et al., 2019; Blumberg et al., 2020), using minibatch scheme (Fatras et al., 2021). Specifically to the Euclidean setting, Vayer et al. (2019b) has introduced solver buiding on a sliced mechanism, and leveraging the observation that the GW problem seems most of the time easy to solve between one-dimensional distributions. More recently, Scetbon et al. (2022) have shown that the low-rank approximations used to speed-up the Sinkhorn-Knopp algorithm were particularly suited for the regularized GW problem, resulting in much more computationally efficient solver. In this work, we propose to build the ideas of Delon and Desolneux (2020) in order to construct OT distances between GMMs that are invariant to isometries and that stay relevant between GMMs of different dimensions. These distances share similarities with the one defined in (Chowdhury et al., 2021), since they rely on a form of quantization of the original data through the GMM representation. One of these distances is a "Gromovization" of the Mixture Wasserstein distance, that we call MGW. We will see that the structured representation of MGW makes it very robust in practice, and permits to design an efficient and scalable solver using a fixed small number of Gaussian components, while keeping competitive precision and running times (when compared to the state-of-tho-art methods described above) when the number of points of the underlying data increases.
Contributions of the paper. In this paper, we introduce two Maanoo--asssesttein type OT distances between GMMs that are designed to be invariant (at least)) to isemetries. More precisely, we introduce in Section a natural Gromov version of the distance introduced by Chen et al. (2018) and Delon and Desolneux (2020), that we call MGW for Mizture Gromov Wasserstein. This distance can be used for applications which only require to evaluate how far the distributions are from each other, without having to identify correspondences between points. However, this formulation does not directly allow to derive an optimal transportation plan between the points. To design a way to define such transportation plan, we define in Section 4 another distance that we call EW for Embedded Wasserstein. This latter turns out to be closely related to the Grommov-aasserseiin distance and coincides with the OT distance introduced by Alvarez-Melis et al. (2019). We show that EW can be adapted to derive a distance and optimal transportation plans between GMMs and we then define heuristic transportation plan for MGW by analogy with EW. Finally, in Section 5, we illustrate the pratical use of our distances on medium-to-large scale problems such as shape matching and hyperspectral image color transfer and we compare the performance of our methods with other recent GW based approaches, both on assessing distances between clouds on points and drawing correspondences between points. All the proofs are postponed to the appendix.
Notation
We define in the following some of the notation that will be used in the paper.
3

---

Published in Transactions on Machine Learning Research (09/2024)
(x,x")a stands for the Euclidean inner-product in R between x and x'. We will use the notation
(r,=²) when the dimension is clear and unambiguous.
|| stands for the Euclidean norm of I E IR. We will use the notation ||||| when the dimension is clear and unambiguous.
Ir(M) denotes the Irace a matrix M.
|| stands for the Frobenius norm of a matrix M, ic. ||M|| ((r(MTM).
|| stands for the nuclear norm of a matrix M, i.e. |[M]] tr((MT M)
The notation (MV) denotes the vector singular values of the matrix M.
Ida is the identity matrix size d d.
For any x E IRd diag(r) denotes the matrix of size x d with diagonal vector x.
Ia stands for any matrix of size d the form diag((=l)
Suppose d 2 d'. For any matrix M of size d x d, we denote by M(d') the submatrix of size d' d' containing the d' first rows and the d' first columns of A.
Let r r and 8 S d'. For any matrix M size x s, we denote by Mid,d") the matrix of size of the form (or M 0) When d= d', will write M14,
We use the notation Sd for the of symmetric matrices of size d x d, S the set of symmetric positive semi-definits matrices, and $4 the set of symmetric positive definite matrices. 1 = si)si denotes the matrix of ones with d' rows and d columns. 13/5d
The notation X ~ means that X is random variable with probability distribution μ.. If is a positive measure on X' and Φ: X' + ) is a mapping, Φ# stands for the push-forward measure of μ by 6, ie, the measure on y such that for any measurable set A of V, Φ* (A)). If μ a positive measure on X' suрр(μ) denotes its support, i.c. the subset of X defined as supp(u) = {a E X for all open set N₂ such that x Nz. p(Nx) > 00.
X and Y are random vectors on IRd and RRd", we use the notation Cov(X,Y) for the matrix of size d x d' of the form [XX EY E[Y))*].
For any positive measure μ. denote by its associated centered measure, i.e. the measure such that if X ~ μ, we have X Ex=y [X] ~ j..
For any m E Rd and any Σ E S we denote by N(m, Σ) the Gaussian measure of mean m and covariance matrix Σ.
For x E X, 8R denotes the Dirac distribution at x.
2
Background Mixture-Wasserstee and Gromov-Wasserstein-tyye distances
We recall in this section the definitions and some important properties of the different OT distances used throughout the paper. For any Polish space X, we write P(X') the probability measures on X. For d d2 and p 2 1, the Wasserstein space (R.)) is defined as the set of probability measures P Rd with finite moment of order p. i.c. such that
[ [=||ddp(x) +00,
with ||-|| being the Euclidean norm Rd
2.1 Mixture- Wassers-Wn distance between GMMs
We present here the distance introduced in Delon and Desolneux (2020), as well as some results that will be useful in the rest of the paper. We denote GMMx (R!) the set of Gaussian mixtures on IRd with less than K components, ie. the set of measures in P(RR) which can be written
K'
where K' < K, (aj...., )T is in Age, and {uu he is family of pairwise distinct Gaussian distributions, each of mean my E Rd and covariance matrix Eg s4. Again, to avoid degeneracy issues where locations
4

---

Published in Transactions on Machine Learning Research (09/2024)
with no mass are accounted for, we will assume that the elements of are all positive. The set of all finite Gaussian mixture distributions on Ied is then written
GMM (R U GMM)(RR).
K20
Note that the condition that the Gaussian components are pairwise distinct ensures the identifiability of the elements of GMM,(R (I²) (Yakowitz and Spragins, 1968), in the sense that two GMMs μ Σx a&pi and V Et bir are equal if and only if K L, and we can reorder the indices such that for all k, as by and μ4 Vk- It can been shown that MM (Re)) is dense in W3 (IR") for the metric Wp, meaning that any measure in W,(R.) can be approximated with any precision for the distance Wp by finite Gaussian mixture distribution. Let f GMMK (R.) and V E GMM)( The Mixturo-Wasse Wasserstein distance of order 2 is defined as
MW2(u, v) (
inf (TEILR.))GGMMM (B24) [
RaxRs ||= ul1 (T)
(MW2)
As for W2 with (2 2W2 defines a metric on GM (R2) (Delon and Desolneux, 2020). In general, the transportation plan solution of the W2 problem is not a Gaussian mixture, thus by restricting the set of admissible couplings, we most the time have (p, v/) > W) It can be shown that the difference between MW2(H,, v) and W2(n, v) is upper-bounded by a term that only depends on the weights and the covariances matrices of the components of the two mixtures. An important property of M MW2 is that it can be written in equivalent form, which had already been introduced in Chen et al. (2018): if μ Σ ax/4: and v = Σ bi then
K
MW2(H,
inf Σ Ww. (2 R.
wEll(a.))
where a = (ai,..., ax) )T, b= (h,,..., br) From a computational point of view, this latter formulation reduces the problem to a simple small-scale discrete optimal transport problem since the W2 distance between Gaussian distributions has a closed form: indcod, recall that if μe N(m,, Σ&) and vs N(m, Ex), then
W (2 v1) = lime- mil1 + tr Ex (E Ex 2 2 (<
(3)
Finally, the respective solutions n² and wn of Problems (MW2) and (2) are linked by the following relationship, for all vx, E Rd Rd,
n" (x 3)
Σ k.
whi (z) y=
(4)
where Tolel the optimal W2 transport map between pa and in and Pps the probability density function of Pk.
2.2 Cromov-Wasserste distance
The Gromov-Wasserster Wasscrstein problem (Mémoli, 2011) can be definod as the following: given two network measure spaces, i.e. triplets of the form (x, where X' is a Polish space, cx X' X R is measurable function and μ E P(X), it aims at finding
inf f / (FEll(..)) xxy Jxxy (cx(xxx') (/ Pd s)dx(s', .'))
GW cx,
()) (w)
(5)
with 2 1. The fundamental metric propertics of GW, have been studied in depth in (Mémoli, 2011; Sturm, 2012; Chowdhury and Mémoli, 2019). When ex and cy are powers of the metrics dx and dy of the base spaces X and V, GW, induces a metric over the space metric measure spaces (i.e. the triplets (X, (X,dx,H)) H)) quotiented by the strong isomorphisms (Sturm, 2012), where one says that two metric mcasure spaces (X' dx, p)) and (Y, dy, p) are strongly isomorphic if there exists an isometric bijection Φ: supp(μ) - supp(v) that transports y into When ex and by are not powers of the metrics of the base spaces, GWp still defines
5

---

Published in Transactions on Machine Learning Research (09/2024)
metric, but this time over the space of network measure spaces quotiented by the weak isomorphisms (Chowdhury and Mémoli, 2019), which are spaces isomorphic for the costs ex and ey relatively to a third space, see (Chowdhury and Mémoli, 2019) for details. Note that in both cases, the metric property of GW, stricly holds only when it takes finite values, and so it is natural to restrict it to the following space
Mp = {(X,, cx, Jxxx Ex (d((r)dp((z)) < +00]
Note that when X' and y are fixed as well as ex and cy, it is natural to see GWp as a distance between the two measures ls and v rather than a distance between the two network measure spaces (X, ex, (4) and (V,y,v). Therefore, we will denote in that case with a slight abuse of notations GW p(u, v) instead of ((X, v)). Finally, in the discrete setting, given a (a)...., 0m) and b = (bi,..., b,,)" being respectively in A. and Ax, and given two non-negative cost matrices C and C of respective size m x m and n x n, the Gromov-Wasserster Wasserstein distance can be written
bW)(a,b) inf Σ 10 - CWPllurruua. well(a,b) (J.kdl
2.3 Other invariant distances
Sturm (2006) has introduced another distance between metric measures spaces which takes the following form Dp((X, dx p), (V, dy, v)) inf (.
(6)
2,0.6
where (X, po) and (V,dy, 2) are two metric measure spaces as defined in Section 2.2, E is third Polish space, and where th: X + 2 and Φ: y - Z are two isometric mappings. More recently, Alvarez-Meltt et al (2019) have roduced another family of invariant OT distances in the Euclidean setting which can also be used to compare distributions on spaces of different dimensions. Initially, Alvarez-Melis et al. (2019) have introduced this OT distance in the sctting where μ and v are both living in the same Euclidean space Rd. Yet, generalizes well to set sings where j and v are living in spaces of different dimensions. Between two measures a and v on Rdd and Rd this reads as
( inf inf / (FETT(E))) hEW xdx
1/2
IW2 IW,He(,v))
1| h(v) h(y) Md))Pam(a.)))
(IW₂)
where H. is a class of mappings from Rd to Rd encoding the invariance. This is a non-convez optimization problem in т and & that becomes convex in æ if h is fixed and becomes also convex in ha if a is fixed and H is a convex set. When d is equal to d' and both measures are centered, notably shown that when v is such that Ey
Mellss et al. (2019) have Ida and when H = H, == {P E Raxd ||] S vd), Problem (IW) is equivalent to the Gromov- Wasserstein problem (5) of order with inner-product costs. Indeed, it can be shown that both problems are equivalent in that case
1/ my (I( *EI[(H,P)
sup
(F-COV)
where for any matrix A of size d x d, ||/||r denotes the Frobenius norm, i.c. Vtr(AT A). Another interesting case is when H H2 == 0(R))= == {P E Rddd prp Ida] is the set of orthogonal matrices of size d x d. In that case, Problem (IW2) is equivalent to
// xy dh dz w) (ell(,.))
sup
T
(+-COV)
gdx
where for any matrix A of size d x d, |A]|+ is the nuclear norm of A, i.e. |||||- both Problems (F-COV) and (=-COV) are non-conuez. These results have been shown by Alvarcz-Melis et al. (2019) in the case where y and v are discrete but can easily be extended to continuous distributions. Observe that problem (+-COV consists in maximizing the sum of the singular values of the cross-covariance matrix I dz v). whereas the Problem (F-COV) consists maximizing the sum of the squared singular values of the cross-covariance matrix. In general, these two problems are not equivalent despite being structurally similar, the example of Figure I illustrates it.
tr((ATA))), =lAA)== Note that

---

Published in Transactions on Machine Learning Research (09/2024)
3.4
0.4
3.1
0.1
-2.2
-0.2
3.3
0.3
-04
-.44
-0.2
0.0
0.2
04
0.8
-04
-0.2
0.0
0.4
5.6
Figure 1: Transport plans between two discrete centered distributions on RR2 composed of three points. Left: optimal coupling given by the maximization of Problem (F-COV). Right: optimal coupling given by the maximization of Problem (+-COV).
3 Gromov-Wasserstein distance between mixture of Gaussians
In this section, we define a Gromov-Woan Wasserstein type distance between Gaussian mixture distributions. This distance is a natural 'Gromovization' Problem (2). Indeed, as has already been observed the literature (Chen et al., 2018; Lambert al., 2022), any Gaussian mixture in dimension can be identified with a probability distribution Rd x $4. i.c. the product space of mcans and covariance matrices. Equivalently, discrete probability distribution on the space of Gaussian distributions
finite Gaussian mixture can be seen N(Rd), which has been proven to be complete metric space when endowed with W2 (Takatsu, 2010) and is furthermore separable since it is a subspace of (2 which is itself a separable metric space when endowed with W2 (Bolley, 2008). Since the theory of optimal transport still applics on measures over non- non-Euclidean spaces (Villani, 2008), follows that Problem (2) can formally be thought as a simple OT problem between two discrete mcasures in (R.).. Thus, one can define dircctly its Gromov version.
Definition 1. Let μ
Ee a&pa and v = E, bier be two Gaussian mixtures respectively on IR and IR we
defme
(G,)( inf Σ well(a.b) k.j.Kl
W2 -
(MGW2)
Unlike MW2, there is no straight/frward equivalent continuous formulation of this latter problem. In particular, it not clear whether Problem (MGW2) is equivalent not to the continuous GW problem between and seen continuous measures on Rd and æd -whree the set of admissible couplings is restricted to Gaussian mixture distributions. Thanks to the identifiability property of the set of finite Gaussian mixtures, we have that each E GMMx (Rd)) associated with a unique discrete distribution f E (Rd) and MGW, between a and coincides with GW2 with squared W2 costs between their associated discrete measures and in (R)).. Finally, note that we have defined MGW, only between GMMs with finite number components because there is general no identifiability property for infinite Gaussian mixtures. As outcome, for given infinite GMM μ on Rdd there might be more than associated continuous measure i (Rd).. For instance, the standard Normal distribution N(O, 1) can naturally be identified in P(V(R) with the Dirac distribution at N(0, 1), but also with the Normal distribution N(O, 1/2) over the parametrized line (N(0, 1/2) e (RR) : 0 E R}, or with N(0, 1) over the parametrized line {ōg E WBB) : 0 R).
(R(Rd) includes the degenerate Gaussian distributions, as for instance the Dirac distributions.

---

Published in Transactions on Machine Learning Research (09/2024)
3.1 Metric properties
Here study the metric properly of 2GW2 that mainly arises from the Gromov- Warn-rstesrrnnn structure of Problem (MGW2). Indeed, the following result is a direct consequence of the theory developed by Sturm (2012).
Proposition 2. In the following, μ Ex Orpa and v == Exbir bi are two GMMs respectively in GMMx and (MMM,(R*)
(i) 2GW2 is non-negative and symmetric.
(8) 2GW2 satisfies the triangle inequality, i.e. for any & GMM,,
MGW2(u,V) S MGW2(u, + MGW2(E,v).
(tt) MGW2(1,2)=0 if and only there erists bijection Φ: {pr - {vih such that Ex (xo(u)) and & is an isomelry for W2, Lee for all k and i smaller than K, W2(((us), M(u)) W2( 14)).
Sketch of proof. The proof of these results mainly consists in applying the theory of Sturm (2012), using the facts that (TR)) is complete (Takatsu, 2010), separable (Bolley, 2008), and metrizable with the Wasserstein distance. See Appendix B.1 for the full proof.
MGW2 defines thus a pscudometric on the set of all finite Gaussian mixtures of arbitrary dimensions, i.c. the set,
GMMx G. (RR),
d21
that is invariant to the mappings é that transform a finite Gaussian mixture Σx=1 4x44i into another finite Gaussian mixture of the form Σ ag4i such that for all k and i smaller than K, W2(M), v4) W2(Hu, pc). A question that arises is: are all these mappings & between M (R=") and P(Rd) always associated with mappings T: Rdd Fd that are isometries for the Euclidean norm and such that Tμh coincides with 6(u) for every μ E GMM (R") Y We can already state the following converse result.
Proposition 3. Let d d', and let T: Rdd + IRd be mapping that is an isometry for the Euclidean norm. Then the mapping Φr: GMMxx + P(Rd) defined as 6T(μ)= = TTH for all E GM (R."), such that for any of the form _ (1())) is in GMM (Re)) and is of the form ΣK D 0k/k, with {uk) k being such that, for all k and smaller than K, W2(",, = (2(Pk, and so MGW,(n, Tæ O.
Sketch of proof. The proof of this result mainly consists in showing that for any μ E GM (R @T(H) is in GMM& (R:)) because T is necessarily affine, as direct consequence of the Mazur-Ulam theorem (Mazur and Ulam, 1932) which implies that any isometry from Rdd to Rd (endowed with the Euclidean norm) is necessarily affine. See Appendix B.2 for the full proof.
Hence, T:: Rdd + Rd is an isometry for the Euclidean norm, then MGW₂ invariant to the mapping Φr : GMM (R ( given for all μ E GM Ao (IE by 6T(4) = Typ. Yet, in general, there exist mappings Φ (2 + W2(Rd) that are isometrics for W₂ and that are not induced by any mapping T: IRd - Rdd that is an isometry for the Fuclidean norm. This has been proven by Kloeckner (2010) the general case when considering isometrics defined all over W2(R,,, but remains truc when restricting isometries defined over subspaces of (RR) as the following example suggests.
Example 4. lel N++ (R) be the set of one-dimensional Gaussian distributions with strictly posilive mean. Let Φ: N++ + (B) be the mapping that swaps the mean and the standard deviation, i.e. such that for any 7 = N(m,, o²) with my > O and o, 0, Φ(7) = N(a m.,). Then & is isometry for Wz. Observe indeed that for y and S in N4 (B), we have
((()))))((())))))))))))))))))(() =aa, - =c) + (my mc)² W2(7,
8

---

Published in Transactions on Machine Learning Research (09/2024)
Thus @ is an isometry for W2. yet is not induced by any isometry of R. Hence there exist mappings from GMM to G.M (Rd) that satisfy the conditions above but which are not induced by isometries for the Euclidean norm from Rdd to Rd.
3.2 2GW, in practice
Using MGW2 on discrete data distributions. Most applications of optimal transport involve discrete data that can be thought as samples drawn from underlying distributions, which are not GMMs in general. those applications, we aim to evaluate an OT distance between two distributions of the form f = (1/M) E, Szg and & = (1/N) E5 By, where {a:) and {ug}y are families of respectively M and N vectors of Rd and IR." Though f and can be thought as mixtures degenerate Gaussian distributions, evaluating directly (GW2(P, P) is not particularly interesting since we have in that case MGW2( f)) = GW2(II-II I1-11 D). However, we can design pseudometric MGWK, K,2 between and & by fitting two GMMs f and with K and & and then sctting MGW D) MGW2(P, v). The approximation of and D by μ and can be done by maximizing the log-likelihood of the GMMs with the EM algorithm (Dempster et al., 1977). Note that if K is chosen too small, the approximations of f and & will be of bad quality and we are likely to observe undesirable behaviors, as for instance having MGWK,2 D) = 0 despite fa and not being equal up to an isometry. Thus, the choice of K must be a compromise between the quality of the approximation given by the GMM and the computational cost. To illustraae the pratical use of MGW2 on simple toy example, we draw 150 samples from the spiral dataset provided in the scikit-learn toolbox' (Pedregosa et al., 2011) and apply rotations with various angles this dataset. We then fit independently GMMs with 20 components 2GW2 between the two obtained GMMs. We also compute GW2 with inner-product as cost functions, MW2 using also 20 Gaussian components and W2. The results can be found in Figure 2. As expected, MGW2 is as 2W2 which is not the case of MW2 and W2.
components
the initial and the larget rotated datasets and we compute
Spiral datasets
Evolution of OT distances
MGW;
GW2
25
MW2
W:
20
15
10
5.
0
44
x2
3
Figure 2: Left first column: spiral datasets (in blue and red) composed of 150 points R22 The red dataset corresponds to points sampled from the distribution of the blue dataset rotated from π. Left second columm: The two corresponding learned GMMs with 20 components via EM algorithm (each color corresponds to a Gaussian component of the GMMs). Right: evolution of MGW2, GW2, MW2, and W2 between the initial distribution (in blue) and the rotated ones in function of the angle of rotation. Experiments are averaged over 10 runs and the colored bands correspond to +/- the standard deviation. This experiment is inspired from Vayer et al. (2019b).
4The package accessible here: littps:/ wwwrikiiit-eeeeeooo wew.iiiteeeaaaare/sstaeeeeee /stable//
9

---

Published in Transactions on Machine Learning Research (09/2024)
Difficulty of designing transportation plan. The 2GW2 problem can be used on discrete data to provide an optimal coupling w" between the Gaussian components of the two Gaussian mixtures and v that approximate the discrete data distributions á and D. However, some applications require an coupling x between the points that compose f and is not straightforward to derive such a transportation plan æ from the plan wn that minimizes the 2GW, problem. naive heuristic approach would be define æ similar way to (4), transporting the Gaussian components using restrictedGGW, transport maps (Delon al., 2022a) instead of W2 transport maps. Yet this approach introduces too many degrees of freedom it consists in transporting the Gaussian components independently, without taking into account the global structure of the mixture, see Figure 3 for an illustrative example. Since selecting the solution that preserves the global structure of the mixture among all the candidates seems to be a difficult combinatorial problem, a better solution to design such plan would be to derive explicitely the isometric transformation that has been implicitely applied to one of the two measures when solving the MGW2 problem. This is the idea behind the embedded Wasserstein distance that we introduce the following section.
Figure 3: Left: two discrete distributions (in gradient of colors) and (in blue) that have been drawn from two GMMs. The colors have been added to f in order to visualize the couplings between and D. Middle and right: two possible solutions of transport of f obtained by plugging the discrete plan that minimizes M MGW, in (4), using restricted-GWy transport maps (Delon et al., 2022a) to transport the Gaussian components. Observe that the middle solution preserves the global structure of the mixture, in the sense that points that are close to each other but associated with different Gaussian components remain close when tranported. This is not the case for the right solution.
4 Embedded Wasserstein distance
In this section, we define an alternative distance to Gromov- Wasserstein also invariant to isometries which specifies the isometric transformation applied to one of the measure when computing the distance.
Definition 5. Let P(Rd) and P(R))). For 2 and 2 1, let denote Isom, (X") the set of all isometries for the Euclidean norm from R to IR. We define
FW2(J,, inf { inf de (B.)
We(n, d4v),
inf (BT')
W2(
(EW₂)
with the convention that the infonimum over an emply set is equal to +00.
Observe that if > d', the set Isoma(Rd") is empty and so FW2(u, v) = inf inf selsom (Rd) W2(n, d4v). In contrast, if d < d', (Boms(RR) is emply and so EW²(u,)) inf (R=" W2(Wpp v). When d d= d', the two infinimums are equivalent. In all what follows, we will suppose without loss of generality that d26 d'
4.1 Properties of EW2
We present here three properties of EW2 which motivate its use between Gaussian mixture models. First, we start by showing that FW2 defines pseudometric that is invariant to isometries.
10

---

Published in Transactions on Machine Learning Research (09/2024)
Proposition 6. 2W2 defines pseudometric on U621 (2 such that for any μ W2(R") and e W2(R"), DW if and only if there exists an isometry Φ: Rd + R such that = 0HHl.
Sketch of proof. Non-negativity and symmetry are straightforward. The triangular inequality can be proved observing first that the infinimum in ø is always achicved, then remarking that EW₂ remains unchanged when one of the two measures is immersed a third Euclidean space of greater dimension than d and d'. This makes 2 closely related to the distance between metric measure spaces introduced by Sturm (2006) presented in Section 2.3. See Appendix B.4 for the full proof.
Now we show that EW2 is equivalent to the OT distance introduced by Alvarez-Melis al. (2019) described in Section 2.3 for a particular choice of transformation space H. In all what follows, denote Va (IR) the Stiefel manifold (James, 1976), i.e. the set of rectangular othogonal matrices of size d x d' such that prp Idg. More precisely, we show the following result.
Proposition 7. Let j E W2(R1.) and v E W2(TR") and let suppose without loss of generality 2 d'. Then,
(W,(D,
inf =Ell(u.)) ae (R"), bERo Jeexse || Py dT),
inf
Moreover this lotter problem is equivalent to
zy d dmleew)
sup *EI((M)))
(*-COV)
Sketch of proof. Equation (7) a consoquence of (Delon al., 2022b, Lemma 3.3) and of the Mazur- Ulam theorem (Mazur and Ulam, 1932), which implies that any isometry from Rd to Id (endowed with the Euclidean norm) is necessarily affine. The equivalence with Problem (=-COV) is then roughly consequence of MAlvarez-Melis et al., 2019, Lemin 4.2) which implies that Problem (7) is achieved in P at P* Ua (ajd.4") VT where Uz E ((R4) and V* ((Rd") are the left and right orthogonal matrices associated with the Singular Value Decomposition (SVD) of w dr( v). See Appendix B.3 for the full proof.
Since Problem (*-COV) is in general not equivalent to Problem (F-COV), the EW2 problem is in general not equivalent to the 2W2 problem with inner-product costs. However, the following result shows that between Gaussian distributions, the two problems share some common solutions.
Proposition 8. Suppose without loss of generality that d'. Let μ N(0, Eo) and v= N(O, En) be two centered Gaussian measures Rd and Rdd, Let Po, Do and Pı, Dj be the respective dingonalizations of Eo (= P.DDP and En (= (= that sort the eigenvalues in non-incceessing order. We suppose that μ is not degenerale, i.e. Eo is non-singular. Then the problem
EW2(u,v)
inf Walu,P441 (R.)
admils solutions of the form (x", P*) with p of the form p Potddd" 1 pf and n" (Ida, T) with T being any affine map such that for all Rd,
T(x) = (f= [d' PHz
In other terms, the solutions of Problem (5) with inner-product costs exhibited in Delon et al. (2022a) are also solutions of Problem (EW2). Parthermore,
(u, - tr(D:) - (0) DF)
Sketch of proof. The proof of this result is inspired from the proof of Equation (3) given by Givens et al. (1984) that is based on Lagrangian analysis. The main difference with the proof of Equation (3) lies the introduction of an additional variable with constraint P V (Rd). See Appendix B.5 for the full proof.
11

---

Published in Transactions on Machine Learning Research (09/2024)
Note that it is not clear if the two problems strictly equivalent only share some common solutions because it is not clear, to the best of our knowledge, that the solutions exhibited above are the only solutions of the 2W2 problem with inner-products costs, see Delon et al. (2022a) for more details. To complete this section, we emphasize that 2 different from the distance proposed in Cai and Lim (2022), that we call here 2W2 for Projection Wasserstein discrepancy. Details about this difference can be found Appendix C.
4.2 Embedded Wasserstein distance between GMMs
Similarly to Delon and Desolneux (2020), one can define an OT distance derived from FW2 when ja and v are GMMs by restricting the set of admissible couplings to themselves GMMss
Definition 9. Let μ E GMMM (R4) and v EE GMML(R")) and suppose that d 2 d'. We define
MEW2(H,V) v) = inf inf ddEloomy (R.))
MW2(H, 041), inf "ETsom,(("")
MWI(V+H..))
(8)
As before, one can reformulate this latter problem by observing that the isomorphic mappings for the Euclidean norm are necessarily of the form Px+ b with P (RR) and E Rd. Similarly to EW2, one can show that the infinimum in Φ is always achieved and that MEW2 satisfies all the properties of a pseudometric G.MMa by simply replacing W2 by MW2 in the proof of Proposition 6. Supposing without loss of generality that d 2 d' and using the equivalent discrete formulation (2) of the MW2 problem, get that for a = Ee 0al4 and E, brn, the problem is equivalent to
inf inf Σ Wa Pavi), k.
(MEW))
(R.) wEl((a.))
where for any k S K and S L, dhe and of are the Gaussian components respectively associated to the centered GMMs and v. Note that ue and of are not necessarily themselves centered.
4.2.1 Numerical solver
This time, is not possible to derive analytically the closed form of the optimal P for Problem (MEW)). However, one can still solve the problem nunerically using an alternate minimization scheme. Indeed, Problem (MEW2) is not convex in P and w, but is convex in if P is fixed and is furthermore simple small-scale discrete OT problem in that case, which motivates the use of an alternating optimization scheme for solving this problem. However, Problem (MEW)) is not convex in P for fixed w because the feasible set, i.e. the Stiefel manifold Va (Rd), is not convex. For a fixed w, the minimization in P can be done by projected gradient descent (Calamai and Moré, 1987), i.c. for a given iterate p(i) and a given w, the next itcrate p(i+1} given by
pfi+1] ky (n (lw(P1t)) (Pat-)))((P()))
-
OP
where (Rd) is the projection mapping the Stiefel manifold, where 1) > 0 and where for all matrices P of size d' x d, (P) = Ex, W2,:W P=)(). Observe that we have, as a byproduct of Proposition 7, that for all P of size d' x d, the projection Ky (Rd) written
Ky (Ra) = Ia Ia
[d.'' VT
where Up e O(Rd) and Vp e ((R") are respectively the left and right orthogonal matrices associated with the SVD of Indood, this projection can be written
By (Ra arg min ||P P1B (EV)
arg min [[[]]]] + || 2tr(p" (p'((p)) P)]
(R
Since for all PE V (R.), |||1| = d', we get that the problem is equivalent to sup (Rd) tr(PTT P) which is maximized when P Up Uplaldee Iald,"d VF, see the sketch of proof of Proposition 7. Finally, this yields to
12

---

Published in Transactions on Machine Learning Research (09/2024)
Algorithm Mixture Embedded Wasserstein solver
Require: μ == ==E* akpk: = Ei bun, p(o) E V (R3)) > O.
1: while not converged do
2: [C]ka + W (1k, v)) for 1. l....,, K; I,...,,
wlit SSoIv-OT(aa))) C)
Solve a classic OT problem. Do projected gradient descent on P.
while not converged do
A p(i-1) cn (p(-
(A) 7: plg &- 1 VT
6:
8: end while
end while
10: return wi, P
Algorithm 1. Note that more involved optimization procedures using the specific structure of the Stiefel manifold could probably be used here (Boumal, 2023).
When f and y are only composed of non-degenerate Gaussian components, one can compute (P) either by using automatic differentiation (Baydin et al., 2018) by using the following technical result, whose proof is postponed to Appendix B.
Lemma 10. Let for any S K, pi N(mOR Σok) with mos Fd and Eor $4 and for any 1 S L, i = N(mun Eu) with mu E IRd' and Eu E = For any w in the K x L simplex, let Jw: Rdxd' - R be the functional defined, for all matrix P of size x d', by
Jo(P) = Σ
W. (2 P4v))
kd
Then for avey full-runk matrix P of size d x d', we have
(J)(P) = = wall Pmmmm] mo - Eos PE (E Eos PE
0P
k.1
Initialization procedure. Since the problem is non-convex, the solution to which Algorithm 1 converges strongly depends the initialization of P. It is therefore crucial to design a good initialization procedure. To do so, propose to use the annealing scheme introduced by Alvarez-Melis et al. (2019). More precisely, we propose to set the initial P as the solution of the following iterative procedure. First we solve an entropic-regulariee W22 problem between the two discrete measures ntt = En a and ve = = bi0mm with a large value of regularization Eg in order to obtain a coupling (1),, Then we set
p(l) Ny (Rd) (EE, Wk. monn!t)
We then solve another entropii-roguulriizdd W2 problem, this time between dº' and p(U,, using smaller value of regularization é = x Eo with a E (0, 1). We obtain thus a new coupling ( and we can then denive p(2) as previously. We repeat this procedure N₂ times until the regularization term EN, becomes small enough. This boils down Algorithm 2.
In practice, set in all our experiments a 0.95 and en 1 in Alvarez-Melis al. (2019). Furthermore observed that in most cases, setting N₂ 10 was sufficient to obtain good initialization of P for Algorithm 1.
4.2.2 Transportation plans and transportation maps
Since (MEW2) has continous equivalent formulation (8), one can derive from any optimal solution (w*, PP) of the former, optimal solution (Tr,, "d)) of the latter. More precisely, we have the hand for all
13

---

Published in Transactions on Machine Learning Research (09/2024)
Algorithm Annealed initialization procedure for mixture embedded Wasserstein
Require: a, b, (mos Jx, {mu)f, 6o 0, (0, p(O) Ia
for i = 1....,,. Nu do
[C]x, + l|mos mull
2:
3:
wlil + =-OT(a, (oT(a,b b, E(1)))
> Solve a regularized OT problem.
At Whi (3]
moxmii
5:
U,E, V + (VD(A)
p(() U (d,d") VT
-8 end for return P
Annealing scheme.
y e Rd 00(1)) (² P* bb, where 6' = [Exxa(X) - P By and on the other hand for all (I, v) E Hd x Rd
x (1, Σ
E WhepPee (1))0 po R.
(9)
where T the optimal W₂ transport map between 1 and p (where we recall that p4 and 1/ are the Gaussian components of the centered GMMs) and s²: Idd + gd is defined for all E Rd as (x) (- (b). As in Delon and Desolneux (2020), it is possible to define a unique assignment of each by setting for all x eell Rd.
(o (2lV T Teaan (x) = (X [Y = x] En axpre (x)pm) (x)
Note that Tmean is not a Monge map since n- is not of the form (Ida,T) T) not equal y and Tmean is not necessarily the gradient of a convex function. When using MEW, obtain an assigmment between two sets { and {y;}, of respectively M and N vectors of IRd and Rd one compute Tmean (2;) for each F3; and then determine which yj is the closest of T meeau (24) using nearest-neighhor algorithm (Fix and Hodges, 1951).
In particular, 4 #4 is
4.2.3 Improving the MGW, method
Inspired by the MEW2 method presented above, we propose in this section to improve the MGW, method by: (i) proposing an annealed scheme similarly to Algorithm 2 in order to reduce the chances of converging to sub-optimal local minima, (ii) designing a transportation plan for MGW, similarly to (9).
only guaranteed to converge towards local minimum when solving it with classic non-regularized GW solver (Peyré et al., 2016). Furthermore, the convergence towards a particular minimum depends strongly on the initialization of the coupling w. Since the discrete GW problem in MGW, is of very small scale and so not costly in itself, we propose, by anology with 2EW,, to use a similar annealing scheme as in Algorithm 2 to reduce the chance of converging to a sub-optimal local minimum. More precisely, this gives the following algorithm.
Annealing scheme. Since Problem (MGW2) is non-convex,
As previously, we set in our experiments that, in toy cases where we know what the global minimum is, that No. = 10 seemed to be sufficient number of iterations to prevent the algorithm from converging towards a sub-optimal minimum.
0.95 and 6o
in Alvarez-Melis et al. (2019) and observed
Designing a transportation plan. Still by analogy with MEW2, one can design transportation plan for 2GW, by defining a matrix PMOW, E V (Rd) and a vector bMow, E Rd, and then replacing v0 o Tw, in (9) by WMGG, Tw,, where for all x E Rd, (MMG,,()) = Pi (x - bMGw.)) More precisely, this can be done the following way. Given two GMMs = Σ. 4x4% and == E,bn bu respectively in GMMx and GMM, (R.) and given the optimal discrete plan wt solution of Problem (MGW2), one can define the matrix
14

---

Published in Transactions on Machine Learning Research (09/2024)
Algorithm Annealed mixture Gromor-Wasserste solver
Require: P Σx akpl: = Σ' bug, (0, 1), ED, w() = ab? 1: [C" t W2 (4;) for k = 1, ,...., K l.... K
2: [C,,, f (2):)))) for I = ,,...,Ly, ,,,,,,,,,,,,,, = I.....L
for n m==1... I,..,, Nut do
w(n) + E-G C C, C", 1n-1,
Solve a regularized GW problem. Anncaling scheme.
En f- aεn-1
6: end for
7: return SOLVE-GW(o, C", C
C (Nu))
b Solve the non-regularized GW problem.
PMOw, as the solution of the following problem
inf (R))
E Wi P#1/),
(10)
where Ue and 1/ the Gaussian component of the centered GMMs j and v, then set bMG Ex~ PMrw, Bylly As above, this problem can be solved numerically by performing a projected gradient descent on P, using either automatic differentiation or Lemma 10. This is also a non-convex optimization problem since (RR( is non-convex and so the solution given by the projected gradient descent depends on the initialization. We propose thus to initialize with the projection on the Stiefel manifold of the discrete cross-covariance matrix between the means of the Gaussian components, i.e.
p d0o MGW, (R4) (En, (..mii)mmmm
Finally, using PMGW, one can define continous plan TMWW, associated with the discrete optimal plan wf solution of the 2GW, problem similarly to (9). We can therefore use MGW2 to transport distributions, We can also, as for M EW2, use MGW₂ to obtain an assignment between two sets
using as previously T meean-
of points.
5 Experiments
In what follows, we use MGW, and MEW, solve Gromov- Wasserstein related tasks on various datasets. More precisely, we apply first the methods on simple toy low-dimensional GMMs. Then, we show that both methods be used to solve relatively efficiently GW related tasks on real datasets in moderate to large scale settings involving sometimes several tens of thousands of points, both for evaluating distances between clouds of points and drawing correspondences between points. In all our experiments, we use the numerical solvers provided by the Python Optimal Transport (POT) package (Flamary et al., 2021) that implements solvers for the ren-ragularized and regularized classic OT and GW problems. Code is available here
5.1 Low dimensional GMMs
In Figure 4, we use again the example of Figure and we derive an optimal transport plan for the MGW, problem as described in Section 4.2.3. We also show the plan obtained by solving the EW, problem. One that with both solutions, the global structure of the distribution is preserved in the sense that points that are closed to each other but in two different Gaussian components have been sent to points that are also close to each other but in different Gaussian components.
5.2 Distances between clouds of points
In this section, we illustrate the usability of our methods to assess distances between clouds of points. First, we reproduce an experiment originally conducted in Rustamov et al. (2013) and presented in Solomon et. al.
TThe package accessible here: https://ppthhonotggtuubbb
15

---

Published in Transactions on Machine Learning Research (09/2024)
data
MGW2
MEW2
Figure 4: Left: two discrete distributions (in gradient of colors) and (in blue) that have been drawn from two GMMs. The colors have been added to in order to visualize the couplings between ft and D. Middle: transport of obtained by solving the MGW2 problem, then deriving PMGW, E V22(R2) by solving Problem (10). Right: transport of f obtained by solving the MEW, problem.
(2016) with the of entropic-regularrizdd GW, that aims to recover the cyclical nature of horse's gallop. Then, we perform comparison between runtimes of MGW, and other methods existing in the literalure that provide a GW-type distance between point clouds.
Galloping horse sequence. Here we repoduce the experiment of the galloping horse, that has been originally conducted in Rustamov et al. (2013) and presented in Solomon al. (2016) with the use of entropie-regnllarized GW. In this experiment, compute a matrix of pairwise distances (either for MGW2 or M MEW2) between 45 meshes representing galloping horse. Then, we conduct Multi-Dimensional Scaling (MDS) (Borg and Groenen, 2005) which roughly can be thought as generalization of PCA of the 45 x 45 matrix of pairwise distances between meshes, in order to plot. each mesh as a 2-dimensional point. Figure 5 shows these 2-dimensional embeddings of the sequence. As observed in Solomon et al. (2016), the interesting part here that those points are positioned in a cyclical fashion, which means that the original set of pairwise distances seem to respect the periodic aspect of the sequence (both for MGW, and M MWW). Each mesh is composed of approximately 9000 vertices and the average time to compute implementation of the entropie-regularized GW solver is around 30 minutes which makes the computation of the full pairwise distance matrix impractical, as mentioned in Solomon et al. (2016). In constrast, when using our methods with GMMs with K = 20 components, it took us only approximately 10 minutes to compute the full distance matrix using M 2W,, and around one hour using M 2EW, these times including the fitting with EM of all the GMMs.
distance when using the POT
MGW2
M MEW2
Data
10 20 10
40
20
40
Medkh
Medi
Figure 5: MDS on the galloping horse animation using the MGW2 distance (left), and the MEW2 distance (middle). Each point corresponds to a given mesh and the meshes are colored in function of their number in the sequence. Right: examples among the 45 meshes that composes the sequence. The computations of both distances have been done by first fitting GMMs with 20 components on each mesh independently.
16

---

Published in Transactions on Machine Learning Research (09/2024)
Local minima. To highlight the importance of using an annealing scheme when deriving MGW2 MEW2. we have reconducted the previous experiment but this time without the annealing schemes described in Algorithm 3 and Algorithm 2. In Figure 6, we plot the evolutions of the values of MGW2 and M EEW between one given fixed mesh and all the others. In both cases, the anncaling scheme scems to be useful to prevent the solver to converge towards sub-optimal mininima. However, if the MGW, solver seems to often converge to the same optimum regardless the use of the annealing scheme, this is not the case of MEW, which, without the annealing initialization procedure (Algorithm 2), converges most of the time to a sub-optimal minimum, much that the periodical aspect doesn't even appear in that case. This experiment also emphasizes the fact that when solving GW problem with classic non-regularized or entropic solvers from Peyré et al. (2016), we are not at all guaranteed to converge towards a global minimum and, more critically, we have in general no ways to know if the solution we converged to is actually optimal or sub-optinal.
MGW2
MEW;
withour conealing
WithO
0.91
with armo.iing
0.05
with anzzaling
0.04
0.05
0.02
0.02
0.01
0.01
0.00
16
20 30
db
10
20
30
40
Mesh id
Meshi
Figure 6: Left: Evolution of MGW3 between the second mesh and all the others, using an annealing scheme (Algorithm 3) in blue, and without the annealing scheme in orange. Right: Evolution of M MEW? between the first mesh and all the others, with the annealing initialization procedure (Algorithm 2) in blue, and without in orange. The computation of both distances has been done by first fitting GMMs with 20 components on each mesh independently.
Runtimes comparison. We perform a comparison between runtimes of MGW2, sliced GW (SGW) (Vayer et al., 2019b), low-rank GW (Seebone et al., 2022), minibatch GW (mbGW) (Fatras et al., 2021), entropic-regularized GW (erGW) (Peyré et al., 2016) and quantized GW (qGW) (Chowdhury et al., 2021) between two 2D random discrete distributions with varying number of points from n 10 to n 10% We use the codes provided by the authors their dedicated Github repositories. Note that MEW, is not included in this comparison as we observed in the previous experiment that this latter method was significantly slower than 2GW,, For MGW,, we use GMMs with respectively K = {10, 20, 50) components. For SGW, we use the implementation on CPU with L = (50, 200} projections. IrGW has a parameter r corresponding to the rank of the coupling matrix. We choose here respectively r = n/100 (this choice is advised by the authors of (Scetbon al., 2022) for IrGW to be good approximation of erGW) and 100 (which yields lincar computational time). For mbGW, we use batches of size batches (these values are ad by the authors of (Fatras et al., 2021)). For erGW, implementations of the method, the first one from POT and the second from Scetbon et al. (2022)7, both with regularization parameter E = 0.1. Finally, for qGW, we use a proportion p = 0.1 of the points as partition block representatives and then we take a Voronoi partition with respoct to these representatives. Note that this latter method only provides a coupling but we reinject it in the GW objective. Results can be found in Figure 7. We can observe that MGW, has similar runtimes as SGW (CPU version) and seems even a bit faster in large scale settings. Several algorithms fail to converge when the number of points is too large. The limits observed are: 10 for both implementations of erGW, 2 x 10* for IGG with r n/100 1/100, and 4 x 1 for 1GGW with 100. In the same way, considering our computational ressources, using gGW to compute a distance between the two point clouds with more than 3 x 10t points was impossible because the two full pair-to-pair distance matrices are becoming too heavy in terms of memory. Note that it is still possible to
50 with k n/ 10
use two different
httpe:///ghhabbcmm//meeeeeeeett http://wwwweeeeeeooaaaaaaooooooe /bloh hwtt://ahahmmbbcaaeaaaaaaaaaaaea Fast://mnnna//aaaaaaaaaa Gromor GromooWassspp Wass.py
17

---

Published in Transactions on Machine Learning Research (09/2024)
compute coupling using qGW afterwards, (Chowdhury et al., 2021) for details, but it is more possible to evaluate the GW objective, which necessarily requires to access the full pair-to-pair distance matrices.
NOW; (K-10)
SGW (L-200)
dGW (p4)
300%; 05-00
NGW 10100
VEWW (RorEE 2017)
10
N0W (5-50)
5OW 100)
cOW
RACW
Figure Runtimes comparison between MGW,, SGW (Vayer et al., 2019b) (CPU), IrGW (Scetbon al., 2022), mbGW (Fatras et al., 2021), erGW (Peyré et al., 2016) and gGW (Chowdhury et al., 2021) between two 2D random discrete distributions with varying number of points from 100 to 10t in log-log scale. The time includes the computation of the pair-to-pair distance matrices and the fitting of the GMMs for MGW, (using scikit-learn).
5.3 Drawing correspondences between points
In this section, we illustrate the usability of our methods to establish correspondences between clouds of points on two shape matching applications.
Source
Target (MGW2)
Target (MEW;)
Figure 8: Shape matching between shapes and their distorted versions. We plot the output of the Tmean transport map applied to the source shape the left for respectively MGW, (middle) and MEW, (right).. GMMs with 20 components have been fitted independently on the source shape and on its noisy version. Colors have been added to visualize where the points have been transported.
Quality of the transport map. We reproduce here an experiment from Chowdhury et al. (2021). The goal is to match 3D meshes from the CAPOD dataset (Papadakis, 2014) with copies of themselves whose
18

---

Published in Transactions on Machine Learning Research (09/2024)
vertices are permuted and perturbed randomly. To do so, fit GMM on each mesh well as on its noisy version, then we derive a discrete coupling w between the Gaussian components using 2GW2 or M 2. Finally, we transport the points using Tmean (see Section 4.2.2). Two examples of matching using GMMs with 20 components can be found in Figure 8. Observe that both methods MGW2 and MEW, seem to be able to recover relatively well the correct matching of the points on these examples. To complete this experiment with quantitative results, we compute, as in Chowdhury et al. (2021), the distortion at cach point x as the distance from its ground truth copy x and its matched point y with a given coupling л. The distortion score of the coupling is then the mean squared distortion. We report class average distortion scores for 2GW2 and MEW2 with GMMs with K = {10, 20, 50} components as well as computation times in Table 1. Since the points at the output of T mean are not exactly corresponding to the points in the noisy version of the shapes, we introduce an additional nearest neighbors step in order to reproject the points onlo the noisy shape. We also report in Table class average distortion scores and computation times for qGW using a proportion p = {0.01, 0.1, 0.2, 0.5) of the points as partition block representatives and using Voronoi partition with respect. to these representatives. In terms of distortions scores, we observe that MGW2 yields similar results that qGW for most classes. We remark that the results we obtained on qGW are consistent with the ones reported in (Chowdhury et al., 2021, Table 1). Note that Chowdhury et al. (2021) have shown that qGW significantly outperforms MREC (Blumberg et al., 2020) and mbGW (Fatras et al., 2021) on this task (which is why we don't 't include them in the comparison). In terms of running time, note that the situation is a little bit different from Figure since the MGW2 method includes here an additional step of deriving a coupling between the points given the coupling between the Gaussian components whereas the qGW method doesn't 't require to assess the GW objective after computing the coupling anymore. We observe that both methods MGW₂ and MEW2 are significantly slower than qGW in that setting. Yet the gap of running time between qGW and our methods seems to reduce as the nunber of points increases. Indeed, using 2GW2 or 2EW, seems to largely decorrelate computation time from the number of points (since the number of Gaussian component stays fixed), without observable deterioration in terms of registration accuracy. These methods therefore make a lot of sense as the number of points increases.
Method
Param
Humans
Planes
Spiders 2664
Cars
Dogs
Trees
Vases
1926
2144
5220
8937
10433
15828
MCW;
10
0.04 (17.5)
0.40 (17.6)
0.03 (17.7)
0.12 (17.8)
0.13 (18.8)
0.10 (18.8)
0.34 (19.1)
20
0.15 (69.4)
0.36 (69.5)
0.04 (69.9)
0.17 (70.8)
0.20 (71.9)
0.10 (71.9)
0.28 (73.1)
50
0.18 (431)
0.10 (428)
0.007 (431)
0.12 (435)
0.20 (437)
0.04 (438)
0.20 (441)
MEW,
0.09 (17.1)
0.37 (22.9)
0.02 (16.3)
0,23 (17.6)
0.20 (24.5)
0.11 (18.3)
0.29 (23.2)
20
0.21 (77.2)
0.39 (66.6)
0.02 (64.1)
0.25 (78.0)
0.20 (85.6)
0.13 (67.1)
0.30 (76.8)
50
0.16 (555)
0.17 (421)
0.009 (397)
0.20 (423)
0.21 (462)
0.08 (465)
0.19 (486)
vGW
0.01
0.25 (0.59)
0.46 (0.78)
0.05 (1.08)
0.24 (3.88)
0.28 (11.3)
0.13 (17.4)
0.28 (32.4)
0.1
0.16 (1.04)
0.10 (1.33)
0,02 (1.84)
0.21 (5.80)
0.02 (16.5)
0.05 (26.8)
0.18 (54.9)
0.2
0.11 (1.65)
0.08 (2.12)
0.01 (2.86)
0.15 (9.25)
0.008 (28.2)
0.04 (53.6)
0.21 (123)
0.5
0.10 (4.39)
0.07 (5.77) 0.007 (7.73)
0.16 (34.9)
0.007 (104)
0.15 (165)
0.22 (418)
Distortion scores (lower is better) and runtimes (in parentheses) for MGW,, MEW,, and qGW. The average number of points in each shape class is provided under the shape class name. Results are listed for several parameter choices of each method. Results have been averaged on 10 runs of the experiment.
Table 1:
Matching human shaped meshes. To demonstrate the usability of our methods in larger scale settings, we use the SHREC'19 dataset (Melzi et al., 2019) that contains human shaped meshes that can sometimes be composed of more than 300000 vertices. Our goal is to draw correspondences between the shapes using only the information of the vertices (the dataset also includes edges). To do so, we first fit independently GMMs with 20 components on each mesh and we derive directly couplings at the scale of the Gaussian components that represent the different parts of the bodies. In such large scale settings, the main bottleneck of the methods in terms of computational time is clearly the fitting of the GMMs that can take at worst 2 minutes for the meshes composed of the highest number of vcrtices. The results are displayed on Figure 9. Observe that in most cases, both methods seem to be able to match correctly the colored parts. Yet in the last row, MEW, matches a leg at the left in red to an arm at the right. This probably implics that the method has been trapped in a local minimum despite the annealing initialization procedure. Finally, note that we presented here cases where the methods performed relatively well, but there are cases where MGW2
19

---

Published in Transactions on Machine Learning Research (09/2024)
2EW, fail to find correct correspondences and exhibit behaviors similar to MEW, in the last row, which suggests that the methods converge sometimes to sub-optimal minima despite the annealing schemes.
Source
Target (MGW2)
Target (MEW2)
Figure 9: Shape matching between human-shaped meshes using MGW2 (middle) and M MEW, (right). Each shape on the left column is matched with the shapes on the same row. GMMs with 20 components have been fitted independently on each shape and the points colored in green and purple correspond to Gaussian components that are matched together when solving MGW2 or MEW,, From left right and top to bottom, the meshes are composed respectively of 84912, 30300, 75000, 273624, 360678, and 360357 vertices.
5.4 Application to hyperspectral image color transfer
The goal here is to reproduce the experiment color transfer conducted in Delon and Desolneux (2020), but this time using a hyperspectral image, i.e an image with more than color channels. More precisely, we aim create an RGB image from an hyperspectral image using the colors of another RGB image To do so, we consider images empirical distributions in the color spaces and solve
Maamoo-Wasserstein problem between the distributions Í NØ Σn du, and & = N Σi deg where M and N are the number of pixels in respectively the hyperspectral image and the RGB image we use as color palette, and {ughk 1M and {ur} are the values at each pixel, i.e for here all I, v E R and for all k, E Rd with >> 3. We thus fit two GMMs j and on respectively P and & and we use MGW2 or MEW,, to derive a mapping Tmeaa: Rd R as described in Section 4.2.2. We apply this process to a hyperspectral image of 512 512 pixels with 15 channels that are displayed in Figure 10 top left. We use color palettes two paintings by Gauguin and Renoir, displayed in Figure 10 top right, that are respectively Manhana no atua (top) and Le déjeuner des canotiers (bottom). These two images are composed of 1024 x 768 pixels. The resulting images Tmean(u) are displayed in Figure 10 bottom (Gauguin at the left and Renoir at the right). For this experiment, we observed that setting the number of Gaussian components to K = 15 was a good compromise between capturing the complexity of the color distributions and obtaining a relatively regular mapping Thean- This experiment shows that MGW, and MEW, can be used large scale scttings: observe indeed that the color distributions and & are composed respectively of approximatively 300000 and 800000 points, which makes the problem intractable with entropic-GW solvers such as Peyré al. (2016) Solomon et al. (2016).
20

---

In term of computation time, the fitting of the two GMMs for the hyperspectral image takes aproximatively one minute against 20 seconds for the GMM for the RGB image. The projected gradient descent becomes rather slow in that setting, which makes it preferable to few updates of P at each step of Algorithm 1 for the computation of 2.W2 Finally, for both methods, it takes around minutes to compute the whole RGB image T mean (u) Published in Transactions on Machine Learning Research (09/2024)
MGW
MEW2
Color transfers between a hyperspectral image with 15 channels (top left) and two paintings by Gauguin and Renoir (vop right, middle right). Bottom line: the obtained RGB images using MGW, and M EW2. For this experiment, we used GMMs with 15 components. Image taken by Francesca Ramacciotti (Alma Mater Studiorum University of Bologna) and Laure Cazals (supported by the European Commission in the framework of the GoGreen project (GA no. 101060768)).
Figure 10:
6 Conclusion and perspectives
In this paper, we have introduced two new OT distances on the set of Gaussian mixture models, MGW, and M 2W2, and have shown that they both can be used solve relatively efficiently Gromov- Wasserstein related problems on Euclidean spaces, especially in moderate-to- sosaatee-t-llaae large scale settings involving several tens of thousands of points. These OT distances also by design particularly suited to settings where there already exists kind of clustcring structure in the data. This being said, if MEW, remains an efficient alternative to the entropic GW solvers proposed by Peyré et al. (2016) and Solomon et al. (2016), we observed that the method was actually slower and perhaps harder to tune than MGW2 for a slighty lower quality of results, and so we believe that MGW, is a better choice in practice. This latter distance part of the family of Gromov- Wasserstein type OT distances that reduce the size of the GW problem, which also includes notably qGW (Chowdhury al., 2021), MREC (Blumberg et al., 2020) and scalable GW (Xu et al., 2019). To the
21

---

Published in Transactions on Machine Learning Research (09/2024)
best of our knowledge, such method specifically based on Gaussian mixture clustering had already been proposed in the literature. Furthermore, our method differs from these three other approaches in the fact that we only need here to solve numerically one single GW problem at the scale of the clusters, using not only the centroid position information but also order 2 statistics. Note in particular that our method has strong similarities with qGW, which uses Voronoi quantizations of the spaces instead of GMMs However, an important difference lies in the fact that the MGW; method primarily provides a distance between clouds of points before providing heuristic coupling between them, while the qGW method directly derives heuristic coupling which can then be reinjected into the objective function to derive distance. For both methods, the optional additional step of computing coupling for 2GW2 and computing a distance for gGW increases the computation time. Consequently, 2GW2 seemns to be a more appropriate choice than qGW for tasks that require only a distance ween clouds of points, see the computation times of Figure 7, while qGW seemss to be a more appropriate choice for tasks that require only a coupling between points. Still, we have shown in our experiments that MGW2 can also provide couplings, with equivalent performance to qGW in terms of accuracy, although it is significantly slower in the small-to-medium scale setting of Table 1. However, we believe that one advantage of 2GW2 over qGW for deriving couplings is that MGW, appears to better decorrelate the number of clusters needed to achieve good accuracy from the nunber of points, which may become an important feature in larger scale settings.
Perspectives for future work MGW, could be easily extended to other type of mixtures as soon as we have an identifiability property between the mixtures and the probability distributions on the space of the distributions that compose the mixtures. If in the Euclidean setting GMMs seem to be versatile enough to represent large classes of concrete and applicd problems, an interesting extension of our work could be consider mixture of distributions on non-Euclidean spaces.
Computationally speaking, the main bottleneck of the method probably comes from the fitting of the GMMs with the costly in large scale settings or as soon as the dimension increases. If the EM algorithm remains invariably the classical algorithm for learning GMMs, some recent approaches (Hosseini and Sra, 2020; Sembach et al., 2022; Pasande et al., 2022) have proposed alternative algorithuss that seem to outperform it. These approaches are based on Riemannian stochastic optimization, leveraging the rich Riemannian structure of the set of positive definite matrices. Another interesting alternative that has been shown to outperform the FM algorithm has been proposed by (Kolouri et al., 2018) and is based on the minimization of the Integrating this in our method could result thus in an approach fully-based optimal transport.
(EM) algorithm (Dempster et al., 1977) which can become relatively
distance.
Another possible limitation of our work lies in the fact that the MGW, solver converges sometimes to sub-optimal local minima. If the anncaled procedure introduced in Section 4.2.3 seems to reduce this issuc, we generally have no guarantee that the solution we have converged to is optimal. This not specific to our method and comes from the gradient descent structure of the classic GW solvers. Still, when solving the GW problem between GMMs rather than solving it directly between the points, it is likely that we increase the probability of converging towards sub-optimal local minimum because inevitably introduce symmetries by simplifying the problem and so we probably increase in the mean time the number of local minima in the GW objective. In the Euclidean setting, the recent work of Ryner et al. (2023) proposes an algorithm for solving the GW problem that is guaranteed to converge toward a global minimum, leveraging the low-rank structure of the cost matrices when the cost functions are the squared Fuclidean distances. A future perspective of work could be therefore to study if a similar idca could be applied for solving the MGW2 problem.
Acknowledgements
This research was funded, in part, by the Agence nationale de la recherche (ANR), through the SOCOT project (ANR-22-CFA40-0017), and the PEPR PDE-AI project (ANR-23-PEIA-0004))
22

---

Published in Transactions on Machine Learning Research (09/2024)
References
Altschuler, J., Bach, F., Rudi, A., and Niles-Weed, J. (2019). Massively scalable Sinkhorn distances via the Nyström method. Advances in neural information processing systems, 32.
Altschuler, J., Bach, F., Rudi, A., and Weed, J. (2018). Approximating the quadratic transportation metric in ncar-lincar time. arXiv preprint Xiv:1810010044
Alvarez-Melis, D. and Jaakkola, T. (2018). Proceedings of the 2018 Conference Empirical Methods in Natural Language Processing, pages 1881-1890.
alignment of word embedding spaces. In
Alvarez-Melis, D., Jegelka, S., and Jaakkola, T. S. (2019). Towards optimal transport with global invariances. In The 22nd International Conference on Artificial Intelligence and Stalislics, pages 1870 1879. PMLR.
Arjovsky, M., Chintala, S., and Bottou, L. (2017). Wasserstein generative adversarial networks. In International Conference on Machine Learning, pages 214-223. PMLR.
Baker, J. (1971). Isometrics normed spaces. The American Mathematical Monthly, 78(6):655-66) 658.
Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M. (2018). Automatic diffcrentiation in machine learning; a survey, Journal of Marchine Learning Research, 18:1 43.
Blumberg, A. J., Carriere, M., Mandell, M. A., Rabadan, R., and Villar, S. (2020). MREC: fast and versatile framework for aligning and matching point clouds with applications to single cell molecular data. stal, 1050:20.
Bolley, F. (2008). Separability and completeness for the Wasserstein distance. Lecture Noles in Mothemalics- Springer- Verlog-, 1934:371.
Borg, I. and Groenen, P. J. (2005). Modern multidimenssional scaling: theory and applications. Springer Science & Business Media.
Boumal, N. (2023). An introduction to oplimizalion on smooth manifolds. Cambridge University Press.
Brogat-Motte, L., Flamary, R., Brouard, C., Rousu, J., and d'Alché Buc, F. (2022). Learning to predict graphs with fused Gromov Kromov-Wassergee Wasserstein barycenters. In International Conference on Machine Learning, pages 2321 2335. PMLR.
Bunne, C., Alvarez-Melis, D., Krause, A., and Jegelka, S. (2019). Learning generative models across incomparable spaces. In International conference on machine learning, pages 851-861. PMLR.
Cai, Y. and Lim, L.-H. (2022). Distances between probability distributions of different dimensions. IEFE Transactions on Information Theory, 68(6):4020 4031.
Calamai, P. H. and Moré, J. J. (1987). Projected gradient methods for linearly constrained problems. Mathematical programming, 39(1):98-116.
Chen, Y., Georgiou, T. T., and Tannenbaum, A. (2018). Optimal transport for Gaussian mixture models. IEEE Access, 7:6269 6278.
Chowdhury, S. and Mémoli, F. (2019). The Gromov-Wasserrtenn distance between networks and stable network invariants. Information and Inference: A Journal of the IMA, 8(4):757-787.
Chowdhury, S., Miller, D., and Needham, T. (2021). Quantized Gromov-Wasserster In Machine Learning and Knowledge Discovery in Databases. Research Track: European Conference, ECML PKDD 2021, Bilboo, Spain, Seplember 13-17, 17, 2021, Proceedings, Parl III 21, pages 811-827. Springer.
Cohen, S. and Guibasm, L. (1999). The Earth mover's distance under transformation sets. In Proceedings of the Seventh IEEE International Conference on Computer Vision, volume 2, pages 1076-1083. IEEE.
Courant, R. (1920). Über die eigenwerte bei den differentialgleichnngnn der mathematischen physik. Mathe- malische Zeilschrift, 7(1-4):1-57.
23

---

Published in Transactions on Machine Learning Research (09/2024)
Courty, N., Flamary, R.,; and Ducoffe, M. (2018). Learning Wasserstein embeddings. In ICLR 2018-6th International Conference on Learning Representations, pages 1-13.
Courty, N., Flamary, R., Tuia, D., and Rakotomamonjy, A. (2016). Optimal transport for domain adaptation. In Transactions on Pattern Analysis and Machine Intelligence, volume 39, pages 1853-1865. IEFF.
Cuturi, M. (2013). Sinkhorn distances: lightspeed computation optimal transport. Advances in neural
information processing systems, 26.
Delon, J. and Desolneux, A. (2020). A Wasserstein-type distance in the space of Gaussian mixture models. SIAM Journal on Imaging Sciences, 13(2):936-97) 970.
Delon, J., Desolneux, A., and Salmona, A. (2022a). Gromov distributions. Journal of Applied Probability, 59(4):1178-1198.
Gromov-wasserstel wasserstein distances between gaussian
Delon, J., Gozlan, N., and Saint-Dizier, A. (2022b). Generalized Wasserstein barycenters between probability measures living different subspaces. Annals of Applied Probability.
Dempster, A. P., Laird, M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the royal statistical society: series B (methodological), 39(1):1-22.
Fatras, K., Zine, Y., Majewski, S., Flamary, R., Gribonval, R., and Courty, N. (2021). Minibatch optimal transport distances; analysis and applications. arXiv e-prints, pages arXiv-2101.
Fischer, E. (1905). Über quadratische formen mit reellen koellizienten. Monatshefte für Mothemaläk und Physik, 16:234-22 249..
Fix, E. and Hodges, J. (1951). Discriminatory analysis: nonparametric discrimination: consistency properties. report. 4. T. USAF School of Aviation Medicine.
Flamary, R., Courty, N., Gramfort, A., Alaya, M. Z., Boisbunon, A., Chambon, S., Chapel, L.., Corenflos, A., Fatras, K., Fournier, N., et al. (2021). POT: Python Optimal Transport. The Journal of Machine Learning Research, 2211:3357333578
Forbes, F., Nguyen, H. D., Nguyen, T. T., and Arbel, J. (2021). Approximate bayesian computation with surrogate posteriors. Preprint hal-03139256. (Cited on pages 5, 65, 69, 71, 101, 103, 251, 265, and 278.).
Forrow, A., Hütter, J.-C., Nitzan, M., Rigollet, P., Schiebinger, G., and Weed, J. (2019). Statistical optimal transport via factored couplings. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1454-2463 PMLR.
Genevay, A., Peyre, G., and Cuturi, M. (2018). Learning generative models with Sinkhorn divergences, In International Conference on Artificial Intelligence and Statistics, volume 84, pages 1608 1617. PMLR.
Givens, C. R., Shortt,, R. M., et al. (1984). A class of Wasserstein metrics for probability distributions. Tn Michigan Mathematical Journal, volume 31, pages 231 240. the University of Michigan.
Hosseini, R. and Sra, S. (2020). An alternative EM for Gaussian mixture models: batch and stochastic Riemannian optimization. Mathematical programming, 181(1):187 223.
James, I. M. (1976). The topology of Stiefel manifolds, volume 24. Cambridge University Press.
Kloeckner, B. (2010). A geometric study of Wasserstein spaces: Fuclidean spaces. Annali della Scuola Normale Superiore di Pisa-Clon Classe di Scienze, ((2):297-323.
Kolouri, S., Nadjahi, K., Simsekli, U., Badeau, R., and Rohde, G. (2019). Generalized sliced Wasserstein distances. Advances in neural information processing systems, 32.
Kolouri, S., Rohde, G. K., and Hoffmann, H. (2018). Sliced Wasscrstein distance for learning Gaussian mixture models. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1427-3436.
24

---

Published in Transactions on Machine Learning Research (09/2024)
Lambert, M., Chewi, S., Bach, F., Bonnabel, S., and Rigollet, P. (2022). Variational inference via Wasserstein gradient flows. In Advances in Neural Information Processing Systems.
Leclaire, A., Delon, J., and Desolneux, A. (2023). Optimal transport between GMMs for texture synthesis. In Scale Space and Variational Methods in Computer Vision: 9th International Conference, SSVM 2023.
Luzi, L., Marrero, C. O., Wynar, N., Baraniuk, R. G., and Henry, M. J. (2023). Evaluating generative networks using Gaussian mixtures of image features. In Proceedings of the IEEE/CVF Winter Conference Applications of Computer Vision, pages 279-288.
Magmus, J. R. and Neudecker, H. (2019). Matriz differential calculus with applications in statistics and econometries. John Wiley & Sons.
Mazur, S. and Ulam, S. (1932). Sur les transformations isométriques d'espaces vectoriels normés. CR Acad. Sci. Paris, 194
Melzi, S., Marin, R., Rodolå, E., Castellani, U., Ren, J., Poulenard, A., Wonka, P., and Ovsjanikov, M. (2019). SHREC 2019: matching humans with different connectivity. In Eurographics Workshop on 3D Object Retrieval, volume 7, page 3. The Furographics Association.
Mémoli, F. (2009). Spectral Gromov Wasserstein distances for shape matching. In 2009 IEEE 12th International Conference on Computer Vision Workshops, TCCV Workshops, pages 256-263. IEFE
Mémoli, F. (2011). Gromov-Wasserster Wasserstein distances and the metric approach to object matching. In Foundalions of Computational Mathematics, volume 11, pages 417 487. Springer.
Papadakis, P. (2014). The canonically posed 3d objects dataset. In Eurographics Workshop on 3D Object Retrieval, pages 33-36.
Pasande, M., Hosseini, R., and Araabi, B. N. (2022). Stochastic first-order learning for large-scale flexibly tied Gaussian mixture model. Xivi preprint 22Xi:2212.055/00
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al. (2011). Scikit-learn: machine learning in python. the Journal of machine Learning research,
Pele, O. and Taskar, B. (2013). The tangent Earth mover's distance. In Geometric Science of Information: First International Conference, GSI 2013, Paris, France, August 28-30, 2013. Proceedings, pages 397-404. Springer.
Petersen, K. B., Pedersen, M. S., et al. (2008). The matrix cookbook. Technical University of Denmark, 7(15):510.
Peyré, G., Cuturi, M., and Solomon, J. (2016). Gromov-Wassersee Wasserstein averaging of kernel and distance matrices. In International conference on machine learning, pages 2664 2672. PMLR.
Rabin, J., Ferradans, S., and Papadakis, N. (2014). Adaptive color transfer with relaxed optimal transport. In 2014 IEEE international conference on image processing (ICIP), pages 4852-4856. IEEE.
Rabin, J., Peyré, G., Delon, J., and Bernot, M. (2012). Wasserstein barycenter and its application to texture mixing. In Scale Space and Variational Methods in Computer Vision: Third International Conference, SSVM 2011, Din-Gedi, Israel, May 29 June 2, 2011, Revised Selecled Papers 3, pages 435 446. Springer.
Rustamov, R. M., Ovsjanikov, M., Azencot, O., Ben-Chen, M., Chazal, F., and Guibas, L. (2013). Map-based exploration of intrinsic shape differences and variability. ACM Transactions on Graphics (TOG), 32(4):1-12.
Ryner, M., Kronqvist, J., and Karlsson, J. (2023). Globally solving the Wasses-Wassersten problem for point clouds in low dimensional euclidean spaces. arXiv preprint arXiv:2307.09051.
Scetbon, M. and Cuturi, M. (2020). Linear time Sinkhorn divergences using positive features. Advances in Neural Information Processing Systems, 18133688-34888
25

---

Published in Transactions on Machine Learning Research (09/2024)
Scetbon, M., Cuturi, M., and Peyré, G. (2021). Low-rank Sinkhorn factorization. In International Conference on Machine Learning, pages 2344-9354. PMLR.
Scetbon, M., Peyré, G., and Cuturi, M. (2022). Linear-time Gromov Stomov-Wasseersterr distances using low rank couplings and costs. In International Conference on Machine Learning, pages 19347-19365. PMTR.
Seguy, V., Damodaran, B. B., Flamary, R., Courty, N., Rolet, A., and Blondel, M. (2017). Large-scale optimal transport and mapping estimation. arXiv preprint
Sembach, L.., Burgard, J. P., and Schulz, V. (2022). A Riemannian Newton trust-region method for fitting Gaussian mixture models. Statistics and Computing, 32(1):8.
Sinkhorn, R. and Knopp, P. (1967). Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics, 21(2):343-38) 348.
Solomon, J., De Goes, F., Peyré, G., Cuturi, M., Butscher, A., Nguyen, A., Du, T., and Guibas, L. (2015). Convolutional Wasserstein distances: efficient optimal transportation on geometric domains. ACM Transactions on Graphics (ToG), 34(4):1-11.
Solomon, J., Peyré, G., Kim, V. G., and Sra, S. (2016). Entropic metric alignment. for correspondence problems. ACM Transactions on Graphies (ToG), 35(4):1-133
Sturm, K-T. (2006). On the geometry of metric measure spaces. i. Acto Math, 196:65-14 131.
Sturm, K.-T. (2012). The space of spaces: curvature bounds and gradient flows on the space of metric measure spaces. arXiv preprint arXiv:1208.0434
Takatsu, A. (2010). On Wasserstein geometry of Gaussian measures. In Probabilistic approach to geometry, pages 463-472. Mathematical Society of Japan.
Tolstikhin, I., Bousquet, O., Gelly, S., and Schölkopf, B. (2018). Wasserstein auto-encoders. In 6th International Conference on Learning Representations (ICLR 2018). OpenReview. net.
Vayer, T., Courty, N., Tavenard, R., and Flamary, R. (2019a). Optimal transport for structured data with application on graphs. In International Conference on Machine Learning, pages 6275 6284. PMLR.
Vayer, T., Flamary, R., Courty, N., Tavenard, R., and Chapel, L. (2019b). Sliced Gromov-Wassersteie Wasserstein. Advances in Neural Information Processing Systems, 32.
Villani, C. (2008). Oplämal transport: old and new, volume 338. Springer Science & Business Media.
Xu, H., Luo, D., and Carin, L. (2019). Scalable Gromov-Wasse Wasserstein Icarning for graph partitioning and matching. Advances in neural information processing systems, 32.
Xu, H., Wang, W., Liu, W., and Carin, L. (2018). Distilled Wasserstein learning for word embedding and topic modeling. Advances in Neural Information Processing Systems, 31.
Yakowitz, S. J. and Spragins, J. D. (1968). On the identifiability of finite mixtures. The Annals of Mathemationl Statistics, (9:1):209 214.
Zhang, Q and Chen, J. (2020). A unitied framework for Gaussian mixture reduction with composite transportation distance. aiv preprint arXiv:2002.08;10
26

---

Published in Transactions on Machine Learning Research (09/2024)
Organization of the supplementary
The supplementary is organized as follows. First, in Appendix A, we show six technical results that will be used throughout the proofs of the paper. In Appendix B, we give the full proofs of the technical results of the paper. Finally, in Appendix C, we give more details on the difference between EW2 and the OT distance introduced by Cai and Lim (2022).
A Technical lemmas
Before lurning to the proofs of the theoretical results, we state here six technical lemmas that will be used throughout the proofs of the results of the paper.
A.1 property of couplings between measures living in different dimensions
First we start by recalling the following result (Delon al., 2022b, Lemma 3.3).
Lemma A1 (Delon et al., 2022b). Let μ W2 (TR)) and E W2(R"") with d not necessarily greater than d', and let T; Rdd + Rd be measurable map. Then x E [1 Tœv) if and only if there is some E I1(, v) such that n' (Ida, TT. In particular, if there ewist a, b 2 0 such that [T(y)ll S a + bll for all y E I2 then
inf / rEll(ae)) X
= inf S | 2 dx,,)
| T(y)[
rEll(u RR
A.2 Isometries in Euclidean spaces
We show the following result, that states that any isometry T: Rdd - Rd for the Euclidean norm is affine and of the form, for all y E Rd T(y) Py b, where b Ed and P is in the Stiefel manifold V (Rd). Lemma A2. Suppose d 2 d', Then o: Rd + Rd is an isometry for the Euclidean norm if and only if there exist Va (R4) and Idd such that for all y E Rdd", 6 is of the form
<ly) = Py + b.
Proof. First observe that for P E (Rd) and Rdd, y + Py + b an isometry since we have, for any y and y' in Idd
Ply - P IPy+b- Pab b12 y0)12 - (1-
(y- 0y)=(y - -') = 11 s/11².
The converse is a consequence of the Mazur-Ula theorem (Mazur and Ulam, 1932) that states in the version of Baker (1971) that an isometry from a real normed space to a strictly convez normed space, iee. a normed space where the unit ball is a stricly convex set, is necessarily affine. Since it is easy to show that the unit ball {a e Ed ||=||| < 1} strictly convex set, we get that for all x E Id Φ of the form y + Py + with P being matrix of size d x d', and E Rdd Moreover we have for all y, D Rdd
( + Py py112 = [| 11122 (y v' P(y-v') v')
Since & is an isometry, it follows that 11 51122 = (u -/- PT - y²) and so pT P Id which concludes
the proof.
A.3 Centering of measures
Lemma A3. Let μ E (2R) and v E (2R")) with d not necessarily greater than d'. Let n and i denote the centered measures associated to μ and and let & be any subset of matrices of size d x d', Then,
inf inf I *ell(n.)) Pep, bERd Zxx3
|| - Py - bl = inf inf
|| - PTll2 (2,
xell(p, PeP JRA
27

---

Published in Transactions on Machine Learning Research (09/2024)
Proof. Denoting mo Ex [X], m1 BY M8], I mu, and ÿ y m1, have for any I1(,, v),
Jev ||0 Py- - b11 fll (1r -
JB JRMXRU 11 Pÿj - b +- mp Pmm (1 (2)
= 11mo - - Pmull 112 Ixexxe 1| Pill (11 (x,
since (ā - Pÿj ma - - Pmi)dT((E,)) u) = Thus it follows,
inf inl
| Py - 2yy-b[2ace b112 (1,
nell(u..) PeR, DER* RdX R"
inf PeB
inf bERU
||ma
Pm1 - b112 + inf /
|| - Pyll y)
TEll(aP)) /xxx"
Observe now that for any P E B. |Imo- Pm1 bbll 00 =
Pm1, which concludes the proof.
A.4 A matrix linear program
Lemma A4. Let K be matrix of size d x d' with Singular Value Decomposition (SVD) K = UKERVE and let B be any compact set of matrices of size d x d'. Then,
trttrpp K) max tr (E(EE)) PEP
sup
PEB
where Σp
with (P)) E Rd denoting the vector of singular values of P. Furthermore the supremum is achieved at P of the form,
P-UKIPV
Proof. Note that this lemma can be proven with a proof similar to the one of Alvarez-Melis et al. (2019, Lemma. 4.2), using the min-max theorem for singular values. Here we offer an alternative proof based on Lagragian analysis. First observe that the supremum is achieved as a direct consequence of the Weierstrass theorem because & compact and the mapping P + tr(P" K) is continuous. For a given E P, let UpE be the SVD of P. The problem can be rewritten
max maXt(VVZ/UUUREEE PEB
Now, let us denote U= Up and V VTVK. Observe that U is in 0(R") and V in O(R"). Using the cyclical permutation of the trace operator, the problem becomes
max (AATUE))ER PEB
Now, for given fixed Σp, we determine which U and V maximize (EETUEEVVT) This problem reads
tr(ETUERVV)
UEO(R"), VeO(R")
The Lagrangian of this problem reads
C(U, V Co Ci) -t - tr(Ci(VTV -lde))
where Co E gd and C E gd are the Lagrange multipliers respectively associated with the constraints U E ((R)) and V E O(R"). The first order condition gives
ΣPVEX 2UC ΣtUEK = 2VC1,
28

---

Published in Transactions on Machine Learning Research (09/2024)
equivalently
{
O'ERVER 2Co DUUXKV = =2V,v
Since Co and Ci are symmetric matrices (because they are associated with symmetric constraints), we get. that both left-hand terms are symmetric. This gives the following conditions
{
UTEPVER = EKVVEEL ΣΛEΣKE VE
Now, observe that when multiplying the first condition at right by U'Ep and multiplying the second condition at left by Σx V we get by combining the two conditions
U EPEPUER {
UTEPVERIS RK = -
or equivalently,
Dp D
U'ERVDK = -D
where Dp = (P)) and Dx = diag(o (K)). Multiplying the first condition at left by V U and the second condition at right by VEŞU, this yields to
{ VD Dr VVUT Dp UE UT U- UUEVDRR
It follows that VD DTDP Dp and U UT D are symmetric matrices and VD VT commutes with Dp and UT DJU commutes with D/g Thus we can deduce that U and V are permutation matrices. Since the singular values are ordered non-increasing order, we deduce that the problem maximized when U Ida and = Idy. This implies that Up UK and Vp VK, which concludes the proof.
Note that Lemma A4 especially useful when the constraint of belonging to the Z can be expressed as a constraint on the singular values. Observe that this is the case of V since for all P E Va (IR",, we have pT P = Ida and equivalent condition of belonging in Var (TR" is that dPP) 1
A.5 Some properties of symmetric matrices
Here we state technical results on symmetric matrices that will be useful in the proofs of the results on Gaussian distributions.
Lemma A5. Let A E St. We denote A and Xa its largest and smallest eigenvalues. For all Rd such that | 1, have
(i) is eigenvector of A associated to A if and only if xT Ax d1
(W) x is an eigenvector of A associated to da if and only if xf Ar = Xa-
Proof. Let x E Rd such |||||| = 1. Since A is symmetric, there exists O e O(R)) and A diag(()u) ssksa) such that HT Ax TT TTOAO x. Denoting the vector O we get thus
A TT
A=
2
Hence it follows that
Ad|l=|/2 2 x Ar S 21112112
with equality if and only if is eigenvector associated with d or Ad-
29

---

Published in Transactions on Machine Learning Research (09/2024)
Lemma A6. Suppose that d d'. Let E be positive semi-defintte (PSD) matrix of size d d' of the form
Eo KT
K Σ
with Eg $4., E1 S and K being rectangular matriz of size d x d', Let S En K2-20' K be the Schur complement of E. Then there erists S d' and B, E (R() such that
K K= L
where U, (RR))) and A, is diogonal positive matric of size such that
E =
Proof. For a given Schur complement S Σ, K=E-t K; we have KTE Ei S. Since En Sd we can deduce that - K E S.. and so that Σ₁ - s We note r the rank of IT K. One can observe that
r d' dd.
where the left- hand side inequality follows from the fact that rk(AB) s min(rk(A), Then, En S can be diagonalized
UA = U,A Σ S KT x =
(11)
with A = diag(A) (Ai: 12) A3 = diag Af. (......,),), X²); and U, E VRRRR") such that U = (Ur Ur-r). From (11), we can deduce that
E KU, KUA Id,,,
where A, is the unique PSD square-root of A.. Let us B, KUUA, such that B, E VR").. It follows that
KU, = Eg B,Ar
Moreover, since U2,,KT 1 =
and No E *))) it follows that KU = 0 and so
= KUUT KUUU," B,A.U.
which concludes the proof.
B Proofs of the theoretical results
B.1 Proof of Proposition
Proof of Proposition Takatsu (2010) has shown that the space of Gaussian distributions ARRR)) is a complete metric space when endowed with W2. Moreover, (V(R)) is separable since it is a subspace of (RR)) which is itself a separable metric space when endowed with W2 (Bolley, 2008). Thus, A M(RR)) is Polish and can directly apply the Gromov- Wasserstein theory developped in Sturm (2012). Let (N(R"), Wa, Ø) and (N(R") W2, D) be two metric measure spaces respectively in MI,. Let us define
D(ū,D) = inf xell(p.o) (NTEE)) J W3(7.7)
/
Applying Sturm (2012, Corollary 9.3), we get that D defines metric over the space of metric measure spaces of the form (WRR)), W2, (i) quotiented by the strong isomorphisms, and thus we get directly that is symmetric, non-negative, satisfies the triangle inequality and D(ū, = 0 and only if there exists bijection Φ: supp(ū) supp(Ø) such that dy R. where for any y and 7 in supp(ū), W2(0(7), o(')) Wa(15
30

---

Published in Transactions on Machine Learning Research (09/2024)
Now observe that if μa Σx asfie and v -= =, bin respectively in GMMx (R*) and GMM,, (RM)(RR) and p Ex an and = Σ byrn are respectively in P(V((R))) and P(N(R()) )) we have
a (2 < t00, SNcresentte) W2(())))4(())77()) = Σ W
h
and
J ) WG(G..))e/()ddic))
Σ byby (v)))) +00
(W(R"), W2, ūi and (N(RR") ), W2, P) are both in M4. Purthermore, we have MGW2(u, D(ū, Hence 2GW2 inherits the metric properties of D, which concludes the proof.
7
B.2 Proof of Proposition
Proof of Proposition 3. First recall that the push-foward measure Tœn with on IRd and T: Rd - Rd is defined as the measure on Idd such that for every Borel set of Rd, Tau(A) = p(T Equivalently, for any measurable map h: Rd R, we have
J (())d(Tw))()>> [n (h ())))dp())
Now observe that for any finite GMM I on Red of the form p = Ex akliks have
Se (h o (
= far E ampu(v)
Σk ak Saw (h (y))dum())
Σ Jzn h(T))d(TTTR)))))
Sax harh( (EK au an(ThMe)()))
and T is of the form Σ ag pR) with Type Gaussian since T necessarily affine a consequence of Lemma A2. Thus, Tæμ is in G. (RR). This proves that ΦT takes its values only in GMM (R²) and that dr aspa) is of the form Σk=1 agk.. Now observe that, for every and smaller than K,
inf ((rr(()) = me Jaaxed || v112
Using two times successively Lemma Al using the fact that T is an isometry an so for any E Rd [T(y)]| = Illll, it follows
inf I
inf xen(T To xdd || x' d'1pd (1) TEll( .4) Red xRØ ll o V'lRdryy,. v/) = W2(mx.)))
Thus, MGW2(H, Tw) a direct consequence of Proposition 2, which concludes the proof.
B.3 Proof of Proposition 7
We prove Proposition 7 before proving Proposition 6 because we will use the former in the proof of the latter.
Proof of Proposition 7. Since suppose d 2 d', we have
v) EW2 (W2 dElsom, (Rd)
inf
W3 (1, dgv).
31

---

Published in Transactions on Machine Learning Research (09/2024)
Let Φ Isoma (R*) for the Euclidean norm. Using Lemma A2, get that there exists P (e() and E Rd such that for all E TR, Φ(y) = Py b. Moreover, we have, using Lemma A1,
J
inf uK, dE (B4) xEN(H RUC
inf
||x: vll2d
{ inf inf
|| ((v)]
@Elsomy (R4) xell(u))) xBd xRd
inf inf =EM af (R*),, heRd JR xR
| Py - bl| dv
which proves Equation (7). Now we show the equivalence with Problem (+-COV). Using Lemma A3, Problem
(7) can be rewritten
EW2(k))= = inf (RV,((R)) well(..)) J ||x - Pyll'dTTT,,, y)
inf
I Jeaxee (11xo) 2 + || - 2(x, Pu)) dr(r,y)
inf PEV. (Rd) rEll(a,))
inf
Since for all P E V (IR"), Pyll doesn't depend on P. we get that the problem is equivalent to
I
sup PeV. (R4) rel((.))) RaxRd
sup
(x, (z,))(z.))
Now observe that for all a E D),,
1 (o, Py)dr((,
[
JRAxR" tr(zy" pT )dm(x,)) = J tr( tr(PT my (dr(,,
xe""
where used the cyclical permutation property of the trace operator. Finally using the lincarity of the trace, get that the problem is equivalent to
(pT pT
xy (g'dd(x,)))
sup PEV. (R') (R') sEl((.))
sup
JRAXR"
or equivalently,
(P.f xyy d dd(
sup PEY (R*) xEll(.e))
sup
P.
Now, using Lemma A4 and using the fact that if P (R() o(P) = la, we get that the problem reduces to
zy dn dm(z, w
sup =EI((E,))
and this is achieved for P = U, ld (dd') VT where U₂ E O(R)) and V4 E C(R") are respectively the left and right orthogonal matrices of the SVD of Jaexa zy dr((,,)) y), which concludes the proof.
B.4 Proof of Proposition 6
Before turning to the proof of Proposition 6, we will prove two useful results. First, we show that the EW2 problem always achieved at an optimal couple (r",
Lemma B1. Let μ E (2RR)) and v E (2(R') and let suppose d 2 d'. Then there exists an optimal isometry o* IRd + Rdd such that EW2(p, = W2(u, div)).
Proof. Using Lemma A3 and Lemma A1, we have that
inf / (w)( = (Rd) rEIl(a.)) J
|| Py dy y)
32

---

Published in Transactions on Machine Learning Research (09/2024)
inf (2, P=D) P
where a and i/ the centered measures associated with a and v. Let us denote J: P W2( P*V) and let us show that J is continuous. For any Po and P1 in (RR",, we have,
((P)) - J(R)) = (7, Po D) - (2 P1æD)1 W2(Po (14
where we used the triangular inequality property of 22. Furthermore,
[ ||= y1l-
inf Pi
inf I xEll(P.P) 32 xR
|Pox Pivl1
where we used Lemma Al twice. Now observe that the coupling (Ida Ide ) )
is in 11( 2), it follows
/ Je"xee xR"" || - dT v) S [ || - dio))
inf rET((..))
Finally, for any pdd we have
|Poz P1z1| ||0| sup
[I(P P1)=11² S IPo 111B311=113
||===1
and so it follows that
P1133 / |,.||dd
(J(P)) J(P))
Since v is in WIR is in (I and so forr | < too. It follows that (J(Pp) - J(P)) - 0 when P - 0 and so is continuous. Moreover, since Var (Rd) is compact (James, 1976), J has a minimum on V (Rd) as a result of the classic Weierstrass theorem that states that any real-valued continous function defined on a compact set achieves its infinimum. Thus, there exists P such that EW2(p, 2/) (2 P1 and setting b* = [x] P*Er~u[Y] and (x)x = P"x + b* for all x E IR, we get that there exists o² E Isoma (Rd) such that EW2(H.1) = W2(u, v4v), which concludes the proof.
Now we show the following results, which imply that EW2 remains unchanged when one of the two measures is immersed in a third Euclidean space of greater dimension than d and d'.
Lemma B2. Let μ E (Rd) and v E (P²)) with d not necessarily greater than d'. Let r 2 max(d, d'} and let (IR". Then, EW2(p, v) EW2 (W2(uu,
Proof. First, using Lemma A2, we get that there exists Pi E Va(R") and b E Rr such that for all x E Rd ((x)= P4. Since d', we have, denoting μ. dau and the centered measures respectively associated with μt thn and u, and using successively Lemma A3 and Lemma A1,
EW2 inf *Ell(Wan.e) PEV (R"), bER" I J ||z Py dl)
inf
inf inf I *ET((NMM)) PeVx( Ja 1|2- -Pyll
inf inf I 1 - PyllddT(2, v) *Ell(f,,) PeVa(( RR
J/ || [ |Puul - sup sup tr(P" Kn) rEl(P,)) (Rr)
J ()z(((dj((a))= + [ [1ol)dd(())- -2 sup sup tr(P" P KT),
xell(a,)) (R))
33

---

Published in Transactions on Machine Learning Research (09/2024)
where Kx Spexner zy (x v). Using the equivalence with Problem (+-COV), we get
sup tr(p" P,K,)= sup nell(P.))
sup xel((.D)) (R")
Now observe that , has the same singular values as K since KT P P K KTK,. Thus || ||Ka| and so WW2(wuu, v) EW,( p), which concludes the proof.
Observe that Lemma B2 highlights that EW2 shares close connections with the distance between metric measure spaces introduced in Sturm (2006) and defined in Equation (6). However it is not clear whether the two distances are strictly equivalent or not because the infinimum in 2 in Equation (6) also includes non-Eute non-Euulidean spaces. However, if we restrict the problem only to Euclidean spaces Z, then Lemma B2 directly implics that the two distances are equivalent. Now we are ready to prove Proposition 6.
Proof of Proposition 6. First observe that non non- negatiiiati is straightforward. Furthermore, observe also that if d + d', symmetry is also straightfoward. Now suppose d d' and observe that that the set V (Rd) coincides with the set of orthogonal matrices O(B²). Thus have
inf inf L =Ell(a..) PEO(B"), bERR RURR"
inf W2(n,oyu) (B1))
||0 Py -1,,
inf inf I rell(p.)) PEO(R"), beRd Re inf W2( velsoma(k")
||TT x y
and so EW2 is also symmetric that case. Before turning to the proof of the two other points, we recall that the infinimum in always achieved, see Lemma B1.
(i) Now prove the triangle inequality. Let > max(d, d', d"), do E Isoma((R)) and for & E W2(R"),, let Φ1 E arg min penoselsmmndd (Rr) W2(0o#11, O. We have, using first Lemma B2, then using the triangle inequality property of W2,
EW2(u, v) EWa(Pu+1,1)
inf dElsom, (IT)) inf aelsomy (2r)
[W2(Po $1 Wa(dinE:
W2(0osp
inf dElsom, (R²) (2(
S EW2(do#u EW
We conclude then by applying Lemma B2 on both terms.
(ii) Suppose without any loss of generality that d 2 d' and suppose EW2(u, v) 0. Since the infinimum in & is achieved, there exists & E Isomar (TR) such that W2(,, 2wv) = 0 and so = Φgv. v. The reverse implication is obvious.
Finally, observe that if μ and have finite order moments, then EW2 necessarily takes finite values, and EW2 defines pseudometric on L4z1 (IR+).
B.5 Proof of Proposition
Proof of Proposition 8. As seen above, Problem (EW2) is equivalent to
sup =El((.)) PeY, (RR) (P, KP)
sup
34

---

Published in Transactions on Machine Learning Research (09/2024)
where K zy" dx(z, v). As in Delon al. (2022a), use the necessary condition for to be in I( v) that is that the covariance matrix En of the law a is PSD matrix, or equivalently that the Schur complement of Dn, i.e. E1 - K2 Kx is also PSD matrix. This gives the following inequality:
sup TE) (R.)
sup (P, S
max KT Keed (
max (P, K)
The rest of the proof is inspired from the proof of the closed-form of the W2 between two Gaussians provided by Givens et al. (1984). We want to solve the following constrained optimization problem
min E E Eo KES
2tr K).
(3")
Using Lemma A6, we can write ()(P as function of B,. This gives the following equivalent constrained optimization problem
min B B,= p=id
-2tr(P 20 B,A.
The Lagrangian of this latter problem reads as
L(Br, P.Co, C1) = BAAU))
+ Id,))
where Co E 8" and C1 E gd' are the Lagrange multipliers respectively associated with the constraints B B, Id, and p P Idg. The first order condition gives
E B,Co {
B,A,A = PC1.
Since Vo, P, Ur, and Ar are full rank, 2 PU,A, is of rank and so Co is also of rank Thus get that
B, = PU,
and so
By B, B, Id,
-o A. pT PT ETE
Thus,
Co= = P
On the other hand, by reinjecting the expression of B, the other first order condition we get
Dd By multiplying this equation by itself we get
Cr.
PT PU, = =G.
Since Of is symmetric we get that po commutes with U, and so Σ₁ S. Moreover, as before we have
tr(pT K) = - - - =tr((E, - (PTS
Using the Courant-Fischer min-max theorem (Courant, 1920; Fischer, 1905) to characterize the eigenvalues of Σj - S, see (Givens et al., 1984, Proposition 7) for details, get that tr(P" K) is maximized when S 0 and so the problem is equivalent to the following problem
tr tr D p)
max (20)
35

---

Published in Transactions on Machine Learning Research (09/2024)
where D2) is any diagonalization of Σn and Do, PT For all R we have ad adllvll 2 - y² E' S on a
where a1...., are the cigenvalues of Eg ordered in non-increasing order. Thus, denoting A...... Ad the eigenvalues of p2 EoP, we get that for all k =d'
ad S A an
Since want to maximize tr( Di.r), set the largest eigenvalue A1 of ET P to 0n We denote In E IRd the eigenvector associated. We have P En Py we get that | || is an eigenvector of No associated with 01 Let As and 3 be any other eigenvalue and its associated cigenvector in the orthonormal basis in which p E is diagonal. We have y y1 = 0 and so vi pa Pt/1 = 0. Thus Pys is orthogonal to Py1- Since | = 1, we get that Pyk also an eigenvector of Eg and so it exists i - such that A u pT En Pyk as- Thus, we conclude that the eigenvalues of the optimal pd P are the d' largest eigenvalues of Eg- Morcover, tr(ßî Do.p) is clearly maximized when Do.r and D2 have their cigenvalues sorted the same order. We conclude then that setting Do = D and Dj = D1, where Do and D₁ are the diagonal matrices associated with the diagonalizations that sort the eigenvalues in non-increasing, maximizes the problem and so it follows that
and || = 1 l| = so using Lemma. A5,
2tr(p" K) =2t(DD
(tr ( D
max E1 KTT Kes (R*)
Finally, observe that when setting K of the form
K = =o( (or's
we have
||||| = =)) (d))
DF).
Morcover, observe that this is the solution of Equation (F-COV) exhibited in (Delon al., 2022a, Lemma 3.2). Thus K* cleary in the feasible set and is optimal. By reinjecting the optimal value in the expression of FW, we get
EW ( (r())) - tr(D:) (t D
Furthermore, using the results of Delon et al. (2022a), we get directly that the optimal plans "t are of the form (Ida, T p with T linear the form
[d',4]
T To (T
pT
Finally, observe that K admits as SVD Po(D TT. For given fixed Tæx we get using Lemma A4, that the optimal pa associated with K* is p== =P, ild.'d pf, which concludes the proof. ]
B.6 Proof of Lemma 10
Proof of Lemma 10. First note that in this proof, we denote dxddd the set of matrices of size d x d' that we distinguish from the set ddd' vector with d x d' coordinates. We set g: P E dd E pT E'' PxR and h: Q E sd + Q² such that for all matrix P of size d x d', we have
S(P) tr(h(g(P)) (P)))
36

---

Published in Transactions on Machine Learning Research (09/2024)
For any matrix A E Raxdd, denote vec(A) Iddd' the vector obtained by stacking the columns of A. Observe that, see (Magnus and Neudecker, 2019) for details, for any function 6z IRdxd' -> RR the Jacobian matrix J J[ø] of & can be defined as, for all P E Rdxd"
J(O[P)) (vPc Övec(P)
Moreover, observe that since f: Raxd R. J[f[[P] E Raid' and
Of(P) vec (J=
aP
where vec-1 is the inverse vector operator, i.e. such that for any A E Raxd', vee vec-((vec(A)) = A. Applying the chain rule to derive f, we have
= = (|
First, we compute (P)) It follows, using formula provided by Petersen et al. (2008) and Magnus and Neudecker (2019),
a pf En =tOP" + p" EPPPE
and so
Ovec(E, Io = / @x m Dvec((( + (2) &K Σ pT Ov)e))((P)
-E; Eo ox Ka) Ovec P) + (E) @K Ej pT Eo) (v)ec((PP) -Ta² Kaa) Øk Σ pm Eg) OveeeP)
where 0g denotes the Kronecker product and for any r, K, is the commutation matrix of size r X T, see (Magmus and Neudecker, 2019) for details. Thus,
J [[g (Lara + (2² 0K i 2o)
Now we compute (Q) Observe that we have for any Q E S
Q!Q==Q
Thus it follows, denoting s: Q+QQQQ
as(Q) + +++as(Q) = =QQ
This latter equation is a Sylvester equation with variable Os(Q), which is equivalent to the following linear system:
(Qà @K Q (ovec(s(Q)) ivec(Q),
where @k stands for the Kronecker sum. If Q is non-degenerate, Q3 Ox Q is also non-degenerate and so in that case
((M(()) = (Q* Ok 1T
Finally, it is casy to see that for R E dd' xd' we have
(R) vee" veoT((d.a))
37

---

Published in Transactions on Machine Learning Research (09/2024)
Thus, denoting A= E pT Pop and observing that is symmetric and full-rank when P is full-rank (since we supposed that Eg and E, are full rank), it follows that for all full-rank matrix P of size d d',
(IA((P) = (Σ} 0x En + (A² 0x ~vo-dv))(
where used that g and (A@x A) were symmetric. Observe now that (A @x AB) ! voc(X), where X E Rd''d' xd' is the unique solution of the following Sylvester equation
A XA* = Idd
Since is symmetric, one can set A QDQ" where Q ((Rx) and D is a diagonal matrix of size d'. The Sylvester equation can be rewritten
D Ida
where Y = Q" XQ. Since A is full-rank, D is invertible and it is easy to see that the unique solution of this latter equation is Y (1/2)D-* and so X (1/2)A and thus
(At ØK veo 2 vee
Moreover, since A is symmetric, we have a (A vee(AA b)
and so it follows that
('(((P) = ( @k Eg E3 vec (A
= voc(En PI A ト
which concludes the proof.
C More details on Projection Wasserstein discrepancy
In this section, we give more details on the difference between EW22 and the OT distance introduced in Cai and Lim (2022) that we call here projection Wasserstein discrepancy. We recall that for E W2(RR)) and v e (2 with 2 d', this OT distance is defined as
(12 = inf W2(0#1, (ETu(R<")
(PW2)
where (TR is the set of all affine mapping from Idd to Rd of the form spx) = pT - b) with P P (R()) and b E Idd One key results of Cai and Lim (2022) to show that PW2 has the following equivalent formulation
PW*(u, v) inf W2(,,
(R
where (RR) is the subset of (R(R defined as WV(R") (E (IR)) there exists d(z) = pT b) with Ve (R.)) and Rd such that d*t m)
Observe that this latter formulation is structurally different of EW2 since for any isometry Φ: Rd + Rd the distribution ΦΦv v is necossarily dogenerate, whereas this not the case for the distribution ç. The difference between EW2 and PW₂ is illustrated Figure Cl.
To highlight even more the difference between EW2 and PW,, we derive an equivalent problem of Problem (PW2). Observe that in that case, the mapping I in (PW2) is not an isometry since it is not injective. As result, the term that previously depended only on the marginal in the developpement, of the square of the Fuclidean distance will now depend on P. More precisely, this gives the following result.
Proposition C1. Let E (2R)))) and u E W2(RR") and let suppose d 2 d'. Problem (PW2) is equivalent to
inf inf rell(p.s) P
(tr(PT E2P) - -2r((2T (tr(PTKa))) KT))
(12)
where Ex -gee xx xxTddxddaaa du((), K2 Seixpe zy dh (x y), and where µ and are the centered measures associated with ps and
38

---

Published in Transactions on Machine Learning Research (09/2024)
x H pT b)
(x-b)
PW2(H,
y Py b
EW2(H,
Figure C1: Link between PW₂ and 2W2 for two distributions μr and v respectively on R² and R. In PW2, P is projected into R by a mapping of the form x + pa - b). In EW2, v is transformed into a degenerate measure (lying on the purple line) on R² with an isometric mapping the form + Py + b.
Proof of Proposition C1. First observe that using Lemma A3, we can consider without any loss generality that μ and v are centered and omit b. Using Lemma A1, it follows
f xz2
(W](H, = inf inf (Bd) af eMI(P (MN)) Jee
||2 ddl (2. v)
inf inf (B") wEll(H.)) 3Ax3N
||P" 1--ylidd(2a.)) 8)
(L [1PTTx//Pddla + S. 11vl/'duuy) - sup
/ =Ell(a.) Rd
inf (R))((")
and the problem is equivalent to
(f IIPT x
inf
- sup I TEM((.)) gd xR
(PTxx) udd)(((()))))
P (Be)) (B:r)
which is itself equivalent to (12), which concludes the proof.
Observe that Problem (12) can be interpreted regularization in P of the EW2 problem since have seen above that this latter equivalent to the following problem
tr(P" Kx).
sup TE (R4)
sup
It can also be interpreted as W2 problem between v and measure jí' which has different second-order moment than j.
39


Reviews:
[{"review_id": "pcTUfePeSy", "reviewer": "Reviewer_oCzP", "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 - The paper is written clearly, and with both mathematical rigor and clarify.  The new approaches are well described, and I found accessible despite the formality.  \n - Combining GMM with Wasserstein and Gromov-Wasserstein is a very natural approach, and so this study is well-warranted.  \n - The empirical study seemed fine, and the benefits of methods were made clear and convincing.  \n\nWeaknesses:\n - The MEW algorithm is not guaranteed to converge to optimum.  Some examples demonstrate the issues.  This is an instance of a known challenging task in this space of shape matching. \n - The MGW seems to do well, but not as well as the quantized-GW approach in terms of time/accuracy trade-off.", "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 introduces two new shape/embedded-distribution matching formulations, and algorithms.  The idea is to leverage the Gromov-Wasserstein framework which finds the optimal transportation plan between the pairwise distances among points within each input.  The Wasserstein part is known to be slow, and the Gromov aspect even slower; so this paper explores first approximating the distributions with a Gaussian Mixture Model of size at most K.  They introduce two specific variants:\n  - MGW is as described above and there is a closed form way to compute the distance that depends only on K.  The authors show it is a psuedometric, and through a few examples shows it approximates the GW pretty well.  \n  - MEW replaces the Gromov aspect (which learns over all pairwise distances) with \"embedded\" aspect that tries to learn the best affine-transformation between compared shapes.  The optimization is non-convex, and the paper offers an annealing approach that seems to work reasonably well.  As a benefit it produces a transportation plan.  \n\nA point of comparison is the quantized Wasserstein that use a Voronoi quantization into K points instead of a Gaussian Mixture Model.  This is also known to converge to the full Wasserstein in the limit, and seems to usually do better in their experiments.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- More clearly state the advantages of MGW (and MEW) over qGW.  \n  - The algorithm labels in the legends of the figures did not quite match the description in the text.  Please fix.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "L5l5UGn3uD", "reviewer": "Reviewer_BZ1F", "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  - The proposed distance is novel.\n  - Further, in practice, the proposed methods appear to be on par or better than existing methods in terms of time and performance.\n\n*Weaknesses*\n  - This distance is limited to GMMs. To compute the distance for discrete samples, one needs to fit a GMM to the data. In turn, this yields an additional source of error that the current paper does not consider.\n  - Further, fitting the GMM turns out to be theoretically delicate, as it has been shown that the optimization landscape of maximum likelihood is not nicely behaved, and the EM algorithm can get stuck at spurious critical points.", "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": "Computing the classical Wasserstein distance between probability distributions has two practical limitations: it is computationally expensive and prevents users from comparing distributions supported on two different probability spaces. This paper presents a way to bypass these drawbacks in the setting where the two probability distributions are (nearly) Gaussian Mixture Models (GMM). The paper has the following contributions:\n\n- (Definition MGW2 distance) The authors combine ideas due to Chen et al. (2018) and Delon and Desolneux (2020) on the Wasserstein distance between GMM and to Memoli (2011) on the Gromov-Wasserstein distance to define a Gromov-Wasserstein distance for GMM, dubed MWG2. The main benefit of this distance is that it reduces to a simpler-looking (yet still hard-to-solve) problem (see Definition 1).\n- The authors also investigate some of the basic properties of MGW2 (see section 3.1), but most of these follow directly from existing results.\n- (Definition MEW2) The definition of MGW2 does not provide a direct way to extract a transportation map between distributions. So, the authors propose a second distance MEW2 (see Definition 8), which bypasses this issue.\n- MGW2 and MEW2 yield non-convex problems, so the authors propose two-stage iterative methods to compute each one of these distances. In the first stage, they initialize using an annealing procedure inspired by a method by Alvarez-Melis et al. (2019) and then refine the solution using alternating minimization.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- One of the critical goals of this work seems to be having an invariant distance (to certain transformations) between probability measures; how the abstract and the introduction are written makes it sound like a nice byproduct. Consider rewriting to emphasize the fact that you aimed for invariance.\n- In the abstract, introduce the abbreviation GMMs (currently, it is used but not introduced).\n- The title has a grammatical mistake it should be Wasserstein (note the s before the t) instead of Wassertein.\n- P2 First line: footnotes should go after punctuation.\n- P2 Change: \"is invariant to important families of transformations\" to \"is invariant to a given family of transformations''\n- P3 Notations: I would call this section Notation (without the s); the same goes for the first line of this section.\n- P3 and P4 This is mostly personal preference, but consider changing the format of the notation section to a couple of paragraphs inline instead of a list. The list occupies too much space and does not match the style of the rest of the paper.\n- P5 Add a reference for the fact that $MW_2$ defines a metric on $GGM_\\infty(R^d)$.\n- P8 Before Proposition 2: 'developped' -> developed.\n- The P8 Question before Proposition 3 is hard to parse; consider rewriting.\n- P8 The Euclidian norm is only one; it is defined in any dimension, but only one. Change the statement to say it is an isometry for the Euclidian norm.\n- P8 Consider adding an Example environment for the example at the end of this page.\n- P9 is hard to follow for the uninitiated; consider rewriting. Explicitly introduce $\\pi*$ and $\\omega*$. Don't call a matrix $Id$ if it is not the identity. Consider dropping the notation $A^{[d', d]}$ and simply write it explicitly in the definition of $T_{GGW_2}^{k, l}$.\n- All of Page 9 feels unnecessary; it could be compressed to a paragraph at best. The authors are explaining how not to do things, which feels distracting.\n- Page 11 Section 4: \"explicit'' is NOT a verb but an adjective or noun.\n- Page 11 Drop the use of \"any'' in the sentence \"we will suppose without any loss of generality that...'' (it's redundant).\n- Page 13 Problems involving the Stiefel Manifold can often be solved via Riemannian optimization methods; consider trying https://www.manopt.org/, see https://www.nicolasboumal.net/book/IntroOptimManifolds_Boumal_2022.pdf for a reference.\n- Page 14 What is the function $\\kappa$ in the definition of $P^{\\{1\\}}$?\n- Page 15 Use \\colon (the command) as supposed to \":'' when defining functions.\n- Page 17 Figure 5 is impossible to parse. Are you somehow comparing the error between the two methods? It seems that you are just plotting labels.\n- For the experiments, use a color-blind-friendly palette.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Lmt9HKbRiK", "reviewer": "Reviewer_QjXg", "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\n1) The authors provide detailed proofs and the paper is well-written and easy to follow\n\n2) These new distances exhibit certain innovation and practicality in handling Gaussian Mixture Models (GMMs) and demonstrate high efficiency under large-scale problems.\n\n3) The experiments of this paper are well-developed and comprehensive, fully demonstrating the superiority of these two methods.\n\n\n\nWeakness:\n\n1) The innovation of this paper may be lacking. Can you highlight the difficulty in your proof?\n\n2) Although the algorithms in this paper are primarily designed for large-scale problems, I am still uncertain about their performance on smaller-scale problems. Since the paper focuses on large-scale settings, even if the performance on smaller scales is not good, it does not affect my overall opinion of the paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper investigates the optimal transport problem and first proposes a novel distance function to extend MW2. This distance, named $MGW_2$, is efficient for estimating distances between GMMs. Furthermore, the authors develop another distance function, the Embedded Wasserstein distance ($MEW_2$), which is also adapted to derive a distance as well as optimal transportation plans between GMMs. These OT distances are also by design particularly suited to settings where there already exists a kind of clustering structure in the data Finally, this paper conducts experiments to validate the efficiency of their methods. The results demonstrate the effectiveness of their method. \nAlthough $MEW_2$ has some shortcomings in practical applications, $MGW_2$ remains effective.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I advise the authors to divide the section 'Discussion' into conclusion and future work.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The main topic is a clear interest:  can Wasserstein computation be improved by first approximating with Gaussian Mixture Modeling.  The answer is mathematically involved, but ultimately the improvements are incredibly modest.  It is worth having this worked out.
Decision: Accept (poster)