lemmatheorem lemma corollarytheorem corollary propositiontheorem proposition
definitiontheorem definition remarktheorem remark lemmaBlemma lemmaB
lemmaClemma lemmaC propositionCproposition propositionC

# Gromov-Wassertein-like Distances in the Gaussian Mixture Models Space

Antoine Salmona1, Julie Delon2, Agnès Desolneux1

(1 ENS Paris-Saclay, CNRS, Centre Borelli UMR 9010
2 Université de Paris, CNRS, MAP5 UMR 8145 and Institut Universitaire de
France
)

###### Abstract

In this paper, we introduce two Gromov-Wasserstein-type distances on the set
of Gaussian mixture models. The first one takes the form of a Gromov-
Wasserstein distance between two discrete distributions on the space of
Gaussian measures. This distance can be used as an alternative to Gromov-
Wasserstein for applications which only require to evaluate how far the
distributions are from each other but does not allow to derive directly an
optimal transportation plan between clouds of points. To design a way to
define such a transportation plan, we introduce another distance between
measures living in incomparable spaces that turns out to be closely related to
Gromov-Wasserstein. When restricting the set of admissible transportation
couplings to be themselves Gaussian mixture models in this latter, this
defines another distance between Gaussian mixture models that can be used as
another alternative to Gromov-Wasserstein and which allows to derive an
optimal assignment between points. Finally, we design a transportation plan
associated with the first distance by analogy with the second, and we
illustrate their practical uses on medium-to-large scale problems such as
shape matching and hyperspectral image color transfer.

#####

Keywords— optimal transport, Wasserstein distance, Gromov-Wasserstein
distance, Gaussian mixture models

##  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 tasks including generative
modeling Arjovsky et al., (2017); Genevay et al., (2018); Tolstikhin et al.,
(2018), domain adaptation Courty et al., (2016), image processing Rabin et
al., (2012, 2014), and embedding learning Courty et al., (2018); Xu et al.,
(2018). For two probability distributions $\mu$ and $\nu$, respectively on two
Polish (i.e complete, separable, metrizable) spaces $\mathcal{X}$ and
$\mathcal{Y}$ and given a lower semi-continuous function
$c:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}_{+}$ called _cost_ ,
optimal transport in its most classic form, aims at solving the following
optimization problem,

$\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}c(x,y)\mathrm{d}\pi(x,y)\;,$
(1)

where $\Pi(\mu,\nu)$ is the set of measures on $\mathcal{X}\times\mathcal{Y}$
with marginals $\mu$ and $\nu$. When $\mathcal{Y}$ is equal to $\mathcal{X}$,
the choice of cost $c_{p}(x,y)=d_{\mathcal{X}}(x,y)^{p}$, with $p\geq 1$ and
$d_{\mathcal{X}}$ being the metric of the space $\mathcal{X}$, induces a
distance between probability distributions with finite $p$-th moments called
the Wasserstein distance $W_{p}$. In the discrete setting, Problem (1) becomes

$\inf_{\omega\in\Pi(a,b)}\sum_{k,l}C_{k,l}\omega_{k,l}\;,$ (2)

where $a=(a_{1},\dots,a_{m})^{T}$ and $b=(b_{1},\dots,b_{n})^{T}$ are
respectively in the $\mathbb{R}^{m}$ and $\mathbb{R}^{n}$ simplexes
$\Delta_{m}$ and $\Delta_{n}$111The simplex $\Delta_{m}$ is the subset of
$\mathbb{R}^{m}$ of $x=(x_{1},\dots,x_{m})^{T}$ such that for all $k\leq m$,
$x_{k}\geq 0$, and $\sum_{k=1}^{m}x_{k}=1$. ,
$\Pi(a,b)=\\{\omega\in\Delta_{m\times n}\,:\;\omega\mathbbm{1}_{n}=a\text{ and
}\omega^{T}\mathbbm{1}_{m}=b\\}$ and $C$ is a non-negative matrix of size
$m\times n$ also called cost.

One weakness of this classical optimal transport approach lies in the fact
that it implicitly assumes that the spaces $\mathcal{X}$ and $\mathcal{Y}$ are
_comparable_ , i.e. that there exists a relevant cost function
$c:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}_{+}$ to compare them.
Yet, this assumption is not always verified. For instance, if
$\mathcal{X}=\mathbb{R}^{d}$ and $\mathcal{Y}=\mathbb{R}^{d^{\prime}}$ with
$d\neq{d^{\prime}}$, the definition of a meaningful cost function
$c:\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}_{+}$ is
not straightforward. Furthermore, some applications such as shape matching
require having an OT distance that is invariant to important families of
transformations, such as translations or rotations or more generally to
_isometries_ 222We say that $\phi:\mathcal{X}\rightarrow\mathcal{Y}$ is an
isometry if for all $(x,x^{\prime})\in\mathcal{X}^{2}$,
$d_{\mathcal{Y}}(\phi(x),\phi(x^{\prime}))=d_{\mathcal{X}}(x,x^{\prime})$..
Even if the two distributions involved in these applications do live in the
same ground space, it is not straightforward to design a cost function such
that the resulting OT distance will be invariant to these families 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). The most commonly
used among them is perhaps the Gromov-Wasserstein (GW) distance Mémoli, (2011)
which has recently received high interest thanks to the flexibility this
approach offers. 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. This approach has been applied to shape
matching Mémoli, (2009) or more generally to correspondence problems Solomon
et al., (2016), word embedding Alvarez-Melis and Jaakkola, (2018), graph
classification Vayer et al., 2019a , graph prediction Brogat-Motte et al.,
(2022), and generative modeling Bunne et al., (2019).

Optimal transport is known to be computationally challenging. Between discrete
distributions, its computation involves solving a linear program that rapidly
becomes costly as soon as the number of points is moderately large. Between
two sets of $n$ points, its computation is done in $O(n^{3}\mathrm{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 order
to solve OT problems. In particular, Cuturi, (2013) proposes to solve a
entropic regularized OT problem using the Sinkhorn-Knopp algorithm Sinkhorn
and Knopp, (1967), reducing the cost of the problem to $O(n^{2})$. Over the
last past years, a large body of works have focused on speeding up the
Sinkhorn-Knopp 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 type of commonly
used solvers are building on sliced mechanisms Rabin et al., (2012); Kolouri
et al., (2019) that leverage the fact that the OT problem between one-
dimentional distributions can be solved using a simple sorting algorithm.
These solvers roughly consist in computing infinitely many linear projections
of the high-dimensional distributions to one-dimensional distributions and
then computing the average of the Wasserstein distances between these one-
dimensional representations. Alternatively, Delon and Desolneux, (2020) have
proposed an OT distance not relying on direct comparison of histograms of
points. First, a _Gaussian mixture model_ (GMM) is fitted on each
distribution, then the two obtained GMMs are compared using a restricted
version of the $W_{2}$ distance where the admissible transportation couplings
$\pi$ must themselves be GMMs. This OT problem has an equivalent formulation,
that had already been proposed by Chen et al., (2018), that allows to solve it
by merely calculating $W_{2}$ distances between pairwise Gaussian components,
which can be done analytically, and then solving a small-scale discrete OT
problem. The main benefit of this latter approach is that the complexity of
the composite OT problem does not depend on the dimension nor on the number of
points but only on the number of components in the GMMs, implying that the
computational cost of this approach comes almost entirely from the fitting of
the GMMs. Although this method probably doesn’t compete with the fastest
recent refinements of the Sinkhorn-Knopp algorithm in terms of pure
computational cost, it provides a relatively scalable and computationally
effective OT distance that is particularly suited when there already exists a
kind of clustering structure in the data. This approach has been used for
texture synthesis Leclaire et al., (2022), evaluating generative models Luzi
et al., (2023), or Gaussian Mixture reduction Zhang and Chen, (2020).

Computationally speaking, the Gromov-Wasserstein problem is known to be much
more costly to solve than the classic linear OT problem. Indeed, solving
numerically the non-regularized GW problem involves solving a classic OT
problem at each iteration. As for linear OT, entropic-regularized solvers have
also been proposed Peyré et al., (2016); Solomon et al., (2016) and involve in
that case solving a regularized linear OT problem at each iteration. 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 a much more computationally efficient
solver. Alternatively, other works have proposed efficient solvers by reducing
the size of the GW problem, either through quantization of input measures
Chowdhury et al., (2021), or by recursive clustering approches Xu et al.,
(2019); Blumberg et al., (2020). Specifically to the Euclidean setting, Vayer
et al., 2019b  has introduced a 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. In this work, we aim to construct
an OT distance between GMMs that is invariant to isometries and that stays
relevant between GMMs of different dimensions, in order to design a relatively
efficient and scalable solver for Gromov-Wasserstein related problems,
especially when there already exists a kind of clustering structure in the
data.

##### Contributions of the paper.

In this paper, we introduce two Gromov-Wasserstein type OT distances between
GMMs. More precisely, we introduce in Section 3 a natural Gromov version of
the distance introduced by Chen et al., (2018) and Delon and Desolneux,
(2020), that we call $MGW_{2}$ for _Mixture 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 allow to derive directly an
optimal transportation plan between the points. To design a way to define such
an optimal transportation plan, we define in Section 4 another distance that
we call $EW_{2}$ for _Embedded Wasserstein_. This latter turns out to be
closely related to the Gromov-Wasserstein distance and to coincides with the
OT distance introduced by Alvarez-Melis et al., (2019) as well as with the
distance between metric measure spaces introduced by Sturm, (2006). We show
that when restricting the admissible transportation couplings $\pi$ to be
themselves GMMs, this distance can be used as another alternative to Gromov-
Wasserstein, not as computationally effective as $MGW_{2}$ but that allows to
derive directly an optimal transportation plan between the points. We also
define a transportation plan for $MGW_{2}$ by analogy with $EW_{2}$. 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. All the proofs are postponed to the appendix.

## Notations

We define in the following some of the notations that will be used in the
paper.

  * •

$\langle x,x^{\prime}\rangle_{d}$ stands for the Euclidean inner-product in
$\mathbb{R}^{d}$ between $x$ and $x^{\prime}$. We will denote $\langle
x,x^{\prime}\rangle$ when there is no ambiguity about the dimension.

  * •

$\|x\|_{\mathbb{R}^{d}}$ stands for the Euclidean norm of
$x\in\mathbb{R}^{d}$. We will denote $\|x\|$ when there is no ambiguity about
the dimension.

  * •

the notation $\text{tr}(M)$ denotes the trace of a matrix $M$.

  * •

$\|M\|_{\mathcal{F}}$ stands for the Frobenius norm of a matrix $M$, i.e.
$\|M\|_{\mathcal{F}}=\sqrt{\text{tr}(M^{T}M)}$.

  * •

$\|M\|_{*}$ stands for the nuclear norm of a matrix $M$, i.e.
$\|M\|_{*}=\mathrm{tr}((M^{T}M)^{\frac{1}{2}})$.

  * •

the notation $\bm{\sigma}(M)$ denotes the vector of singular values of the
matrix $M$.

  * •

$\operatorname{Id}_{d}$ is the identity matrix of size $d\times d$.

  * •

$\widetilde{\operatorname{Id}}_{d}$ stands for any matrix of size $d\times d$
of the form $\textup{diag}((\pm 1)_{1\leq i\leq d})$

  * •

Suppose $d\geq{d^{\prime}}$. For any matrix $M$ of size $d\times d$, we denote
$M^{({d^{\prime}})}$ the submatrix of size ${d^{\prime}}\times{d^{\prime}}$
containing the ${d^{\prime}}$ first rows and the ${d^{\prime}}$ first columns
of $A$.

  * •

Let $r\leq d$ and $s\leq{d^{\prime}}$. For any matrix $M$ of size $r\times s$,
we denote $M^{[d,{d^{\prime}}]}$ the matrix of size of the form
$\begin{pmatrix}M&0\\\ 0&0\end{pmatrix}$. When $d={d^{\prime}}$, we will
denote $M^{[d]}$.

  * •

For any $x\in\mathbb{R}^{d}$, $\mathrm{diag}(x)$ denotes the matrix of size
$d\times d$ with diagonal vector $x$.

  * •

We denote $\mathbb{S}^{d}$ the set of symmetric matrices of size $d\times d$,
$\mathbb{S}_{+}^{d}$ the set of symmetric positive semi-definite matrices, and
$\mathbb{S}_{++}^{d}$ the set of symmetric positive definite matrices.

  * •

$\mathbbm{1}_{{d^{\prime}},d}=(1)_{\begin{subarray}{c}1\leq
i\leq{d^{\prime}}\\\ 1\leq j\leq d\end{subarray}}$ denotes the matrix of ones
with ${d^{\prime}}$ rows and $d$ columns.

  * •

The notation $X\sim\mu$ means that $X$ is a random variable with probability
distribution $\mu$.

  * •

If $\mu$ is a positive measure on $\mathcal{X}$ and
$\phi:\mathcal{X}\rightarrow\mathcal{Y}$ is a mapping $\phi_{\\#}\mu$ stands
for the push-forward measure of $\mu$ by $T$, i.e. the measure on
$\mathcal{Y}$ such that for any measurable set $\mathsf{A}$ of $\mathcal{Y}$,
$\phi_{\\#}\mu(\mathsf{A})=\mu(\phi^{-1}(\mathsf{A}))$.

  * •

If $\mu$ is a positive measure on $\mathcal{X}$ , $\mathrm{supp}(\mu)$ denotes
its support, i.e. the subset of $\mathcal{X}$ defined as
$\mathrm{supp}(\mu)=\\{x\in\mathcal{X}\ |\ \text{there exits }N_{x}\text{ open
such that }x\in N_{x}\text{ and }\mu(N_{x})>0\\}$.

  * •

If $X$ and $Y$ are random vectors on $\mathbb{R}^{d}$ and
$\mathbb{R}^{d^{\prime}}$, we denote $\text{Cov}(X,Y)$ the matrix of size
$d\times{d^{\prime}}$ of the form
$\mathbb{E}\left[(X-\mathbb{E}[X])(Y-\mathbb{E}[Y])^{T}\right]$.

  * •

For any positive measure $\mu$, we denote $\bar{\mu}$ its associated centered
measure, i.e. the measure such that if $X\sim\mu$, we have
$X-\mathbb{E}_{X\sim\mu}[X]\sim\bar{\mu}$.

  * •

For any $m\in\mathbb{R}^{d}$ and any $\Sigma\in\mathbb{S}^{d}_{+}$, we denote
$\mathrm{N}(m,\Sigma)$ the Gaussian measure of mean $m$ and covariance matrix
$\Sigma$.

  * •

For $x\in\mathcal{X}$, $\delta_{x}$ denotes the Dirac distribution at $x$.

##  2 Background : Mixture-Wasserstein and Gromov-Wasserstein-type 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
$\mathcal{X}$, we write $\mathcal{P}(\mathcal{X})$ the set probability
measures on $\mathcal{X}$. For $d\geq 1$ and $p\geq 1$, the Wasserstein space
$\mathcal{W}_{p}(\mathbb{R}^{d})$ is defined as the set of probability
measures $\mu$ on $\mathbb{R}^{d}$ with finite moment of order $p$, i.e. such
that

$\int_{\mathbb{R}^{m}}\|x\|^{p}\mathrm{d}\mu(x)<+\infty\;,$ (3)

with $\|.\|$ being the Euclidean norm on $\mathbb{R}^{d}$.

###  2.1 Mixture-Wasserstein distance between GMMs

We present here the distance introduced in Delon and Desolneux, (2020), as
well as various results of this latter paper. We denote
$GMM_{K}(\mathbb{R}^{d})$ the set of Gaussian mixtures on $\mathbb{R}^{d}$
with less than $K$ components, i.e. the set of measures in
$\mathcal{P}(\mathbb{R}^{d})$ which can be written

$\mu=\textstyle{\sum\limits_{k=1}^{K^{\prime}}a_{k}\mu_{k}\;,}$ (4)

where $K^{\prime}\leq K$, $a=(a_{1},\dots,a_{K^{\prime}})^{T}$ is in
$\Delta_{K^{\prime}}$, and $\\{\mu_{k}\\}_{k}$ is a family of pairwise
distinct Gaussian distributions, each of mean $m_{k}\in\mathbb{R}^{d}$ and
covariance matrix $\Sigma_{k}\in\mathbb{S}^{d}_{+}$. Again, to avoid
degeneracy issues where locations with no mass are accounted for, we will
assume that the elements of $a$ are all positive. The set of all finite
Gaussian mixture distributions on $\mathbb{R}^{d}$ is then written

$\textstyle{GMM_{\infty}(\mathbb{R}^{d})=\underset{K\geq
0}{\bigcup}GMM_{K}(\mathbb{R}^{d})\;.}$ (5)

Note that the condition that the Gaussian components are pairwise distinct
ensures the identifiability of the elements of $GMM_{\infty}(\mathbb{R}^{d})$
Yakowitz and Spragins, (1968), in the sense that two GMMs
$\mu=\sum^{K}_{k}a_{k}\mu_{k}$ and $\nu=\sum^{L}_{l}b_{l}\nu_{l}$ are equal if
and only if $K=L$, and we can reorder the indices such that for all $k$,
$a_{k}=b_{k}$ and $\mu_{k}=\nu_{k}$. It can been shown that
$GMM_{\infty}(\mathbb{R}^{d})$ is dense in $\mathcal{W}_{p}(\mathbb{R}^{d})$
for the metric $W_{p}$, meaning that any measure in
$\mathcal{W}_{p}(\mathbb{R}^{d})$ can be approximated with any precision for
the distance $W_{p}$ by a finite Gaussian mixture distribution. Let $\mu\in
GMM_{K}(\mathbb{R}^{d})$ and $\nu\in GMM_{L}(\mathbb{R}^{d})$. The Mixture-
Wasserstein distance of order $2$ is defined as

$MW_{2}(\mu,\nu)=\left(\inf_{\pi\in\Pi(\mu,\nu)\cap
GMM_{\infty}(\mathbb{R}^{2d})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x-y\|^{2}\mathrm{d}\pi(x,y)\right)^{\frac{1}{2}}\;.$
($MW_{2}$)

As for $W_{2}$ with $\mathcal{W}_{2}(\mathbb{R}^{d})$, $MW_{2}$ defines a
metric on $GMM_{\infty}(\mathbb{R}^{d})$. In general, the transportation plan
solution of the $W_{2}$ problem is not a Gaussian mixture, thus by restricting
the set of admissible couplings, we most of the time have
$MW_{2}(\mu,\nu)>W_{2}(\mu,\nu)$. It can be shown that the difference between
$MW_{2}(\mu,\nu)$ and $W_{2}(\mu,\nu)$ is upper-bounded by a term that only
depends on the weights and the covariances matrices of the components of the
two mixtures. Finally, $MW_{2}$ can be written in an equivalent form, which
had already been introduced in Chen et al., (2018): if
$\mu=\sum_{k}^{K}a_{k}\mu_{k}$ and $\nu=\sum_{l}^{L}b_{l}\nu_{l}$, then

$MW_{2}^{2}(\mu,\nu)=\inf_{\omega\in\Pi(a,b)}\sum_{k,l}\omega_{k,l}W_{2}^{2}(\mu_{k},\nu_{l})\;,$
(6)

where $a=(a_{1},\dots,a_{K})^{T}$, $b=(b_{1},\dots,b_{L})^{T}$. From a
computational point of view, this latter formulation reduces the problem to a
simple small-scale discrete optimal transport problem since the $W_{2}$
distance between Gaussian distributions has a closed form: indeed, recall that
if $\mu_{k}=\mathrm{N}(m_{k},\Sigma_{k})$ and
$\nu_{l}=\mathrm{N}(m_{l},\Sigma_{l})$, then

$W_{2}^{2}(\mu_{k},\nu_{l})=\|m_{k}-m_{l}\|^{2}+\mathrm{tr}\left(\Sigma_{k}+\Sigma_{l}-2\left(\Sigma_{l}^{\frac{1}{2}}\Sigma_{k}\Sigma_{l}^{\frac{1}{2}}\right)^{\frac{1}{2}}\right)\;.$
(7)

###  2.2 Gromov-Wasserstein distance

The Gromov-Wasserstein problem Mémoli, (2011) can be defined as the following:
given two Polish spaces $\mathcal{X}$ and $\mathcal{Y}$, two measurable
functions
$c_{\mathcal{X}}:\mathcal{X}\times\mathcal{X}\shortrightarrow\mathbb{R}$ and
$c_{\mathcal{Y}}:\mathcal{Y}\times\mathcal{Y}\shortrightarrow\mathbb{R}$, and
two probability measures $\mu\in\mathcal{P}(\mathcal{X})$ and
$\nu\in\mathcal{P}(\mathcal{Y})$ with $p\geq 1$, it aims at finding

$GW_{p}(c_{\mathcal{X}},c_{\mathcal{Y}},\mu,\nu)=\left(\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathcal{X}\times\mathcal{Y}}\int_{\mathcal{X}\times\mathcal{Y}}|c_{\mathcal{X}}(x,x^{\prime})-c_{\mathcal{Y}}(y,y^{\prime})|^{p}\mathrm{d}\pi(x,y)\mathrm{d}\pi(x^{\prime},y^{\prime})\right)^{\frac{1}{p}}\;.$
($GW_{p}$)

The fundamental metric properties of $GW_{p}$ have been studied in depth in
Mémoli, (2011); Sturm, (2012); Chowdhury and Mémoli, (2019). When
$c_{\mathcal{X}}$ and $c_{\mathcal{Y}}$ are powers of the metrics of the base
spaces $\mathcal{X}$ and $\mathcal{Y}$, $GW_{p}$ induces a metric over the
space of _metric measure spaces_ (i.e. the triplets
$(\mathcal{X},d_{\mathcal{X}},\mu)$) quotiented by the _strong isomorphisms_
Sturm, (2012), where one say that two metric measures spaces
$(\mathcal{X},d_{\mathcal{X}},\mu)$ and $(\mathcal{Y},d_{\mathcal{Y}},\mu)$
are strongly isomorphic if there exists an isometric bijection
$\phi:\mathrm{supp}(\mu)\rightarrow\mathrm{supp}(\nu)$ that transports $\mu$
into $\nu$. When $c_{\mathcal{X}}$ and $c_{\mathcal{Y}}$ are not powers of the
metrics of the base spaces, $GW_{p}$ still defines a metric, but this time
over the space of _network measure spaces_ (i.e. the triplets
$(\mathcal{X},c_{\mathcal{X}},\mu$)), quotiented by the _weak isomorphisms_
Chowdhury and Mémoli, (2019), which are spaces isomorphic for the costs
$c_{\mathcal{X}}$ and $c_{\mathcal{Y}}$ relatively to a third space, see
Chowdhury and Mémoli, (2019) for details. Note that in both cases, the metric
property of $GW_{p}$ stricly holds only when it takes finite values, and so it
is natural to restrict it to the following space

$\mathbb{M}_{p}=\textstyle{\\{(\mathcal{X},c_{\mathcal{X}},\mu)\,:\;\int_{\mathcal{X}\times\mathcal{X}}c^{p}_{\mathcal{X}}(x,x^{\prime})\mathrm{d}\mu(x)\mathrm{d}\mu(x^{\prime})<+\infty\\}}\;.$
(8)

Finally, in the discrete setting, given $a=(a_{1},\dots,a_{m})^{T}$ and
$b=(b_{1},\dots,b_{n})^{T}$ being respectively in $\Delta_{m}$ and
$\Delta_{n}$, and given two non-negative cost matrices $C^{x}$ and $C^{y}$ of
respective size $m\times m$ and $n\times n$, the Gromov-Wasserstein distance
can be written as

$GW_{p}(C^{x},C^{y},\alpha,\beta)=\inf_{\omega\in\Pi(a,b)}\sum_{i,j,k,l}|C^{x}_{i,k}-C^{y}_{j,l}|^{p}\omega_{i,j}\omega_{k,l}\;.$
(9)

###  2.3 Other invariant distances

In the Euclidean setting, Alvarez-Melis et al., (2019) have introduced another
family of invariant OT distances 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 setting where $\mu$ and
$\nu$ are both living in the same Euclidean space $\mathbb{R}^{d}$. Yet, it
generalizes well to settings where $\mu$ and $\nu$ are living in spaces of
different dimensions. Between two measures $\mu$ and $\nu$ on
$\mathbb{R}^{d}$, this reads as

$IW_{2}(\mathcal{H},\mu,\nu)=\left(\inf_{\pi\in\Pi(\mu,\nu)}\inf_{h\in\mathcal{H}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-h(y)\|^{2}\mathrm{d}\pi(x,y)\right)^{\frac{1}{2}}\;,$
($IW_{2}$)

where $\mathcal{H}$ is a class of mappings from $\mathbb{R}^{d^{\prime}}$ to
$\mathbb{R}^{d}$ encoding the invariance. This is a _non-convex_ optimization
problem in $\pi$ and $h$ that becomes convex in $\pi$ if $h$ is fixed and
becomes also convex in $h$ if $\pi$ is fixed and $\mathcal{H}$ is a convex
set. When $d$ is equal to ${d^{\prime}}$ ant both measures are centered,
Alvarez-Melis et al., (2019) have notably shown that when $\nu$ is such that
$\mathbb{E}_{Y\sim\nu}[YY^{T}]=\operatorname{Id}_{d}$ and when
$\mathcal{H}=\mathcal{H}_{1}:=\\{P\in\mathbb{R}^{d\times
d}\,:\;\|P\|_{\mathcal{F}}\leq\sqrt{d}\\}$, Problem ($IW_{2}$) is equivalent
to the Gromov-Wasserstein problem ($GW_{p}$) of order $2$ with inner-product
costs. Indeed, it can be shown that both problems are equivalent in that case
to

$\sup\limits_{\pi\in\Pi(\mu,\nu)}\left\|\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}xy^{T}\mathrm{d}\pi(x,y)\right\|_{\mathcal{F}}\;,$
($\mathcal{F}$-COV)

where for any matrix $A$ of size $d\times d$, $\|A\|_{\mathcal{F}}$ denotes
the Frobenius norm, i.e. $\sqrt{\mathrm{tr}(A^{T}A)}$. Another interesting
case is when
$\mathcal{H}=\mathcal{H}_{2}:=\mathbb{O}(\mathbb{R}^{d}):=\\{P\in\mathbb{R}^{d\times
d}\,:\;P^{T}P=\operatorname{Id}_{d}\\}$ is the set of orthogonal matrices of
size $d\times d$. In that case, Problem ($IW_{2}$) is equivalent to

$\sup\limits_{\pi\in\Pi(\mu,\nu)}\left\|\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}xy^{T}\mathrm{d}\pi(x,y)\right\|_{*}\;,$
($*$-COV)

where for any matrix $A$ of size $d\times d$, $\|A\|_{*}$ is the nuclear norm
of $A$, i.e. $\|A\|_{*}=\mathrm{tr}((A^{T}A)^{\frac{1}{2}})$. Note that both
Problems ($\mathcal{F}$-COV) and ($*$-COV) are _non-convex_. These results
have been shown by Alvarez-Melis et al., (2019) in the case where $\mu$ and
$\nu$ 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 $\int xy^{T}\mathrm{d}\pi(x,y)$, whereas
the Problem ($\mathcal{F}$-COV) consists in 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, as the example of
Figure 1 illustrates it.

|
---|---
Figure 1: Transport plans between two discrete centered distributions on
$\mathbb{R}^{2}$ composed of three points. Left: optimal coupling given by the
maximization of Problem ($\mathcal{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-Wasserstein type distance between Gaussian
mixture distributions. This distance is a natural "Gromovization" of Problem
(6). Indeed, as it has already been observed in the literature Chen et al.,
(2018); Lambert et al., (2022), any Gaussian mixture in dimension $d$ can be
identified with a probability distribution on
$\mathbb{R}^{d}\times\mathbb{S}^{d}_{+}$, i.e. the product space of means and
covariance matrices. Equivalently, a finite Gaussian mixture can be seen as a
discrete probability distribution on the space of Gaussian distributions
$\mathcal{N}(\mathbb{R}^{d})$333$\mathcal{N}(\mathbb{R}^{d})$ includes the
degenerate Gaussian distributions, as for instance the Dirac distributions.,
which has been proven to be a complete metric space when endowed with $W_{2}$
Takatsu, (2010) and is furthermore separable since it is a subspace of
$\mathcal{W}_{2}(\mathbb{R}^{d})$ which is itself a separable metric space
when endowed with $W_{2}$ Bolley, (2008). Since the theory of optimal
transport still applies on measures over non-Euclidean spaces Villani, (2008),
it follows that Problem (6) can formally be thought as a simple OT problem
between two discrete measures in $\mathcal{P}(\mathcal{N}(\mathbb{R}^{d}))$.
Thus, one can define directly its Gromov version.

######  Definition .

Let $\mu=\sum_{k}a_{k}\mu_{k}$ and $\nu=\sum_{l}b_{l}\nu_{l}$ be two Gaussian
mixtures respectively on $\mathbb{R}^{d}$ and $\mathbb{R}^{d^{\prime}}$, we
define

$MGW_{2}^{2}(\mu,\nu)=\inf_{\omega\in\Pi(a,b)}\sum_{i,j,k,l}|W_{2}^{2}(\mu_{i},\mu_{k})-W_{2}^{2}(\nu_{j},\nu_{l})|^{2}\omega_{i,j}\omega_{k,l}\;.$
($MGW_{2}$)

Unlike $MW_{2}$, there is no straightforward equivalent formulation of this
latter problem. In particular, it is not clear whether Problem ($MGW_{2}$) is
equivalent or not to the continuous GW problem between $\mu$ and $\nu$ \- seen
as continuous measures on $\mathbb{R}^{d}$ and $\mathbb{R}^{d^{\prime}}$ \-
where the set of admissible couplings is restricted to Gaussian mixture
distributions. In the rest of the chapter, we distinguish the distribution
$\tilde{\mu}\in\mathcal{P}(\mathcal{N}(\mathbb{R}^{d}))$ from its associated
GMM $\mu\in GMM_{\infty}(\mathbb{R}^{d})$. Thanks to the identifiability
property of the set of finite Gaussian mixture, we have that each $\mu\in
GMM_{\infty}(\mathbb{R}^{d})$ is associated with a unique discrete
distribution $\tilde{\mu}\in\mathcal{P}(\mathcal{N}(\mathbb{R}^{d}))$ and
$MGW_{2}$ between $\mu$ and $\nu$ coincides with $GW_{2}$ with squared $W_{2}$
as cost functions between the associated measures $\tilde{\mu}$ and
$\tilde{\nu}$. More generally, one can define for any metric measure space of
the form $(\mathcal{N}(\mathbb{R}^{d}),W_{2}^{2},\tilde{\mu})$ and
$(\mathcal{N}(\mathbb{R}^{d^{\prime}}),W_{2}^{2},\tilde{\nu})$, the following
continuous GW problem,

$\inf_{\pi\in\Pi(\tilde{\mu},\tilde{\nu})}\int_{\mathcal{N}(\mathbb{R}^{d})\times\mathcal{N}(\mathbb{R}^{d^{\prime}})}\int_{\mathcal{N}(\mathbb{R}^{d})\times\mathcal{N}(\mathbb{R}^{d^{\prime}})}|W_{2}^{2}(\gamma,\gamma^{\prime})-W_{2}^{2}(\zeta,\zeta^{\prime})|^{2}\mathrm{d}\pi(\gamma,\zeta)\mathrm{d}\pi(\gamma^{\prime},\zeta^{\prime})\;,$
(10)

where $\tilde{\mu}$ and $\tilde{\nu}$ can be possibly thought as infinite
mixture of Gaussians. However, there is in general no identifiability property
for infinite Gaussian mixture and so for a given GMM $\mu$ on
$\mathbb{R}^{d}$, they might be more than one associated measure $\tilde{\mu}$
on $\mathcal{N}(\mathbb{R}^{d})$. For instance, the standard Normal
distribution $\mathrm{N}(0,1)$ can naturally be identified in
$\mathcal{P}(\mathcal{N}(\mathbb{R}))$ with the Dirac distribution at
$\mathrm{N}(0,1)$, but also with the Normal distribution $\mathrm{N}(0,1/2)$
over the parametrized line
$\\{\mathrm{N}(\theta,1/2)\in\mathcal{N}(\mathbb{R})\,:\;\theta\in\mathbb{R}\\}$,
or with $\mathrm{N}(0,1)$ over the parametrized line
$\\{\delta_{\theta}\in\mathcal{N}(\mathbb{R})\,:\;\theta\in\mathbb{R}\\}$.

###  3.1 Metric properties

Here we study the metric property of $MGW_{2}$ that mainly arises from the
Gromov-Wasserstein structure of Problem ($MGW_{2}$). Indeed, the following
result is a direct consequence of the theory developped by Sturm, (2012).

######  Theorem 1.

In the following, $\mu=\sum_{k}a_{k}\mu_{k}$ and $\nu=\sum_{l}b_{l}\nu_{l}$
are two GMMs respectively in $GMM_{K}(\mathbb{R}^{d})$ and
$GMM_{L}(\mathbb{R}^{d^{\prime}})$.

  * (i)

$MGW_{2}$ is non-negative and symmetric.

  * (ii)

$MGW_{2}$ satisfies the triangle inequality, i.e. for any $\xi\in
GMM_{S}(\mathbb{R}^{d^{\prime\prime}})$,

$MGW_{2}(\mu,\nu)\leq MGW_{2}(\mu,\xi)+MGW_{2}(\xi,\nu)\;.$ (11)

  * (iii)

$MGW_{2}(\mu,\nu)=0$ if and only if there exists a bijection
$\phi:\\{\mu_{k}\\}_{k}\rightarrow\\{\nu_{l}\\}_{l}$ such that
$\nu=\sum_{k}a_{k}\phi(\mu_{k})$ and $\phi$ is an isometry for $W_{2}$, i.e.
for all $k$ and $i$ smaller than $K$,
$W_{2}(\phi(\mu_{k}),\phi(\mu_{i}))=W_{2}(\mu_{k},\mu_{i})$.

###### Sketch of proof..

The proof of this results mainly consists in applying the theory of Sturm,
(2012), using the facts that $\mathcal{N}(\mathbb{R}^{d})$ is complete
Takatsu, (2010), separable Bolley, (2008), and metrizable with the Wasserstein
distance. ∎

$MGW_{2}$ defines thus a pseudometric on the set of all finite Gaussian
mixtures of arbitrary dimensions, i.e. the set,

$\textstyle{\mathcal{GMM}_{\infty}=\underset{d\geq
1}{\bigsqcup}GMM_{\infty}(\mathbb{R}^{d})\;,}$ (12)

that is invariant to the mappings $\phi$ that transform a finite Gaussian
mixture $\sum_{k=1}a_{k}\mu_{k}$ into another finite Gaussian mixture of the
form $\sum_{k=1}^{K}a_{k}\nu_{k}$ such that for all $k$ and $i$ smaller than
$K$, $W_{2}(\nu_{k},\nu_{i})=W_{2}(\mu_{k},\mu_{i})$. A question that arises
is: are all these mappings $\phi$ associated with mappings $T$ that are
isometries for the Euclidean norm and such that $T_{\\#}$ coincides with
$\phi$? We can already state the following converse result.

######  Proposition .

Let $d\geq{d^{\prime}}$, and let
$T:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ be a mapping that is an
isometry for the Euclidean norm. Then the mapping
$\phi_{T}:GMM_{\infty}(\mathbb{R}^{d^{\prime}})\rightarrow\mathcal{P}(\mathbb{R}^{d})$
defined as $\phi_{T}(\mu)=T_{\\#}\mu$ for all $\mu\in
GMM_{\infty}(\mathbb{R}^{d^{\prime}})$, is such that for any $\mu$ of the form
$\Sigma_{k=1}^{K}a_{k}\mu_{k}$, $\phi_{T}(\mu)$ is in
$GMM_{\infty}(\mathbb{R}^{d^{\prime}})$ and is of the form
$\Sigma_{k=1}^{K}a_{k}\nu_{k}$, with $\\{\nu_{k}\\}_{k=1}^{K}$ being such
that, for all $k$ and $i$ smaller than $K$,
$W_{2}(\nu_{k},\nu_{i})=W_{2}(\mu_{k},\mu_{i})$ and so
$MGW_{2}(\mu,T_{\\#}\mu)=0$.

###### Sketch of proof..

The proof of this result mainly consists in showing that for any $\mu\in
GMM_{\infty}(\mathbb{R}^{d^{\prime}})$, $\phi_{T}(\mu)$ is in
$GMM_{\infty}(\mathbb{R}^{d})$ because $T$ is necessarily affine, as a direct
consequence of the Mazur-Ulam theorem Mazur and Ulam, (1932) which implies
that any isometry from $\mathbb{R}^{d^{\prime}}$ to $\mathbb{R}^{d}$ (endowed
with the Euclidean norm) is necessarily affine. ∎

Hence, if $T:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ is an isometry
for the Euclidean norm, then $MGW_{2}$ is invariant to the mapping
$\phi_{T}:GMM_{\infty}(\mathbb{R}^{d^{\prime}})\rightarrow
GMM_{\infty}(\mathbb{R}^{d})$ such that for all $\mu\in
GMM_{\infty}(\mathbb{R}^{d^{\prime}})$, $\phi(\mu)=T_{\\#}\mu$. Yet, in
general, there exist mappings
$\phi:\mathcal{W}_{2}(\mathbb{R}^{d^{\prime}})\rightarrow\mathcal{W}_{2}(\mathbb{R}^{d})$
that are isometries for $W_{2}$ and that are not induced by any mapping
$T:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ that is an isometry for
the Euclidean norm. This has been proven by Kloeckner, (2010) in the general
case when considering isometries defined all over
$\mathcal{W}_{2}(\mathbb{R}^{d})$ but remains true when restricting to
isometries defined over subspaces of $\mathcal{N}(\mathbb{R}^{d})$ as the
following example suggests: let $\mathcal{N}_{++}(\mathbb{R})$ be the set of
one-dimensional Gaussian distributions with strictly positive mean. Let
$\phi:\mathcal{N}_{++}(\mathbb{R})\rightarrow\mathcal{N}_{++}(\mathbb{R})$ be
the mapping that swaps the mean and the standard deviation, i.e. such that for
any $\gamma=\mathrm{N}(m_{\gamma},\sigma_{\gamma}^{2})$ with $m_{\gamma}>0$
and $\sigma_{\gamma}>0$,
$\phi(\gamma)=\mathrm{N}(\sigma_{\gamma},m_{\gamma}^{2})$. Then $\phi$ is an
isometry for $W_{2}$. Observe indeed that for $\gamma$ and $\zeta$ in
$\mathcal{N}_{++}(\mathbb{R})$, we have

$W_{2}(\phi(\gamma),\phi(\zeta))=(\sigma_{\gamma}-\sigma_{\zeta})^{2}+(m_{\gamma}-m_{\zeta})^{2}=W_{2}(\gamma,\zeta)\;.$
(13)

Thus $\phi$ is an isometry for $W_{2}$, yet $\phi$ is not induced by any
isometry of $\mathbb{R}$. Hence there exist mappings from
$GMM_{\infty}(\mathbb{R}^{d^{\prime}})$ to $GMM_{\infty}(\mathbb{R}^{d})$ that
satisfy the conditions above but which are not induced by isometries for the
Euclidean norm from $\mathbb{R}^{d^{\prime}}$ to $\mathbb{R}^{d}$.

###  3.2 $MGW_{2}$ in practice

##### Using $MGW_{2}$ 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. In those applications, we aim to evaluate an OT distance between two
distributions of the form $\hat{\mu}=(1/M)\sum_{i}\delta_{x_{i}}$ and
$\hat{\nu}=(1/N)\sum_{j}\delta_{y_{j}}$ where $\\{x_{i}\\}_{i}$ and
$\\{y_{j}\\}_{j}$ are families of respectively $M$ and $N$ vectors of
$\mathbb{R}^{d}$ and $\mathbb{R}^{d^{\prime}}$. Though $\hat{\mu}$ and
$\hat{\nu}$ can be thought as mixtures of degenerate Gaussian distributions,
evaluating directly $MGW_{2}(\hat{\mu},\hat{\nu})$ is not particularly
interesting since we have in that case
$MGW_{2}(\hat{\mu},\hat{\nu})=GW_{2}(\|.\|^{2},\|.\|^{2},\hat{\mu},\hat{\nu})$.
However, we can design a pseudometric $MGW_{K,2}$ between $\hat{\mu}$ and
$\hat{\nu}$ by fitting two GMMs $\mu$ and $\nu$ with $K$ components on
$\hat{\mu}$ and $\hat{\nu}$ and then setting
$MGW_{K,2}(\hat{\mu},\hat{\nu})=MGW_{2}(\mu,\nu)$. The approximation of
$\hat{\mu}$ and $\hat{\nu}$ by $\mu$ and $\nu$ 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 $\hat{\mu}$ and
$\hat{\nu}$ will be of bad quality and we are likely to observe undesirable
behaviors, as for instance having $MGW_{K,2}(\hat{\mu},\hat{\nu})=0$ despite
$\hat{\mu}$ and $\hat{\nu}$ 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 illustrate the pratical use of
$MGW_{2}$ on a simple toy example, we draw $150$ samples from the spiral
dataset provided in the scikit-learn toolbox444The package is accessible here:
https://scikit-learn.org/stable/. Pedregosa et al., (2011) and we apply
rotations with various angles on this dataset. We then fit independently GMMs
with $20$ components on the initial and the target rotated datasets and we
compute $MGW_{2}$ between the two obtained GMMs. We also compute $GW_{2}$ with
inner-product as cost functions, $MW_{2}$ using also $20$ Gaussian components
and $W_{2}$. The results can be found in Figure 2. As expected, $MGW_{2}$ is
rotation-invariant as $GW_{2}$ which is not the case of $MW_{2}$ and $W_{2}$.

Spiral datasets |  Evolution of OT distances
---|---

|
---|---
|

Figure 2: Left first column: spiral datasets (in blue and red) composed of
$150$ points of $\mathbb{R}^{2}$. The red dataset corresponds to points
sampled from the distribution of the blue dataset rotated from $\uppi$. Left
second column: The two corresponding learned GMMs with $20$ components via EM
algorithm (each color corresponds to a Gaussian component of the GMMs). Right:
evolution of $MGW_{2}$, $GW_{2}$, $MW_{2}$, and $W_{2}$ 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 .

##### Difficulty of designing a transportation plan.

The $MGW_{2}$ problem can be used on discrete data to provide an optimal
coupling between the Gaussian components of the two Gaussian mixtures $\mu$
and $\nu$ that approximate the discrete data distributions $\hat{\mu}$ and
$\hat{\nu}$. However, some applications require an optimal coupling between
the points that compose $\hat{\mu}$ and $\hat{\nu}$. It is not straightforward
to design such a transportation plan associated with the plan that minimizes
the $MGW_{2}$ problem. More precisely, for two GMMs $\mu=\sum_{k}a_{k}\mu_{k}$
and $\nu=\sum_{l}b_{l}\nu_{l}$, the discrete $MW_{2}$ problem between the
associated distributions
$\tilde{\mu}\in\mathcal{P}(\mathcal{N}(\mathbb{R}^{d}))$ and
$\tilde{\nu}\in\mathcal{P}(\mathcal{N}(\mathbb{R}^{d}))$ is equivalent to
restricting the set of coupling to be GMMs in the continuous $W_{2}$ problem
between $\mu$ and $\nu$. Thus, there exists a direct relationship between the
optimal couplings $\omega^{*}$ and $\pi^{*}$ associated with these two latter
problems. Indeed, when the $mu_{k}$ are all non-degenerate distributions, we
have for any $x,y\in\mathbb{R}^{d}$,

$\pi^{*}(x,y)=\sum_{k,l}\omega^{*}_{k,l}p_{\mu_{k}}(x)\delta_{y=T^{k,l}_{W_{2}}(x)},$
(14)

where for $\mu_{k}=\mathrm{N}(m_{k},\Sigma_{k})$,
$p_{\mu_{k}}(x)=(2\uppi)^{-\frac{m}{2}}|\Sigma_{k}|^{-\frac{1}{2}}\exp\left[-\frac{1}{2}(x-m_{k})^{T}\Sigma^{-1}_{k}(x-m_{k})\right]$
is the density of $\mu_{k}$ in $x$, and $T^{k,l}_{W_{2}}$ is the optimal
affine transportation map between $\mu_{k}$ and
$\nu_{l}=\mathrm{N}(m_{l},\Sigma_{l})$ associated with $W_{2}$, i.e. for all
$x\in\mathbb{R}^{d}$,

$T^{k,l}_{W_{2}}(x)=m_{l}+\Sigma_{k}^{-\frac{1}{2}}(\Sigma_{k}^{\frac{1}{2}}\Sigma_{l}\Sigma_{k}^{\frac{1}{2}})^{\frac{1}{2}}\Sigma_{k}^{-\frac{1}{2}}(x-m_{k})\;.$
(15)

However, in the case of $MGW_{2}$, there doesn’t exist to the best of our
knowledge, any equivalent continuous formulation of the problem and so there
is formally no such plan $\pi^{*}$ associated with the discrete optimal
coupling $\omega^{*}$ that minimizes the $MGW_{2}$ problem. Yet, supposing
without any loss of generality that $d\geq{d^{\prime}}$, and that $\nu\in
GMM_{L}(\mathbb{R}^{d^{\prime}})$, one could still define by analogy a plan
$\pi$ relatively to $\omega^{*}$ using (14) and replacing $T_{W_{2}}$ by
suited transportation maps between Gaussian components. A naive approach would
be to simply replace $T_{W_{2}}$ by the affine transportation map
$T_{GGW_{2}}$ associated with the Gaussian plan that minimizes Problem
($GW_{p}$) with squared Euclidean distance costs restricted to Gaussian
couplings, that we have exhibited in Salmona et al., (2021):

$T^{k,l}_{GGW_{2}}(x)=m_{l}+P_{l}\left(\widetilde{\operatorname{Id}}_{d^{\prime}}D_{l}^{\frac{1}{2}}{D_{k}^{({d^{\prime}})}}^{-\frac{1}{2}}\right)^{[{d^{\prime}},d]}P_{k}^{T}(x-m_{k})\;,$
(16)

where $(P_{k},D_{k})$ and $(P_{l},D_{l})$ are the respective diagonalizations
of $\Sigma_{k}\ (=P_{k}D_{k}P_{k}^{T})$ and $\Sigma_{l}\
(=P_{l}D_{l}P_{l}^{T})$ that sort the eigenvalues in non-increasing order,
where $\widetilde{\operatorname{Id}}_{d^{\prime}}$ is any matrix of the form
$\mathrm{diag}((\pm 1)_{i\leq{d^{\prime}}})$, where $D_{k}^{({d^{\prime}})}$
denotes the submatrix of size ${d^{\prime}}\times{d^{\prime}}$ containing the
${d^{\prime}}$ first rows and the ${d^{\prime}}$ first columns of $D_{k}$, and
where for any matrix $A$ of size $r\times s$ with $r\leq{d^{\prime}}$ and
$s\leq d$, we denote $A^{[{d^{\prime}},d]}$ the matrix of size
${d^{\prime}}\times d$ of the form

$A^{[d,{d^{\prime}}]}=\begin{pmatrix}A&0\\\ 0&0\end{pmatrix}\;.$ (17)

Yet this approach implies that each Gaussian component is transported
independently of the others and so offers too many degrees of freedom. Observe
indeed that $T_{GGW_{2}}$ is defined up to
$\widetilde{\operatorname{Id}}_{d^{\prime}}$ that can be any matrix of the
form $\mathrm{diag}((\pm 1)_{i\leq{d^{\prime}}})$, implying that we have
$2^{d^{\prime}}$ possibilities for each Gaussian component. For each
component, since we want that points that are close to each other but
associated with different Gaussian components remain close when tranported, we
need to determine, relatively to all the other components, which of the
$2^{d^{\prime}}$ possibilities is the correct one. To illustrate this problem,
we show in Figure 3 a $2$-dimensional example where we derive two different
transport maps using $T_{GGW_{2}}$ but each time with a different
$\widetilde{\operatorname{Id}}_{2}$. If the transport on the middle of Figure
3 preserves the global structure of the distribution since two points that are
close to each other but associated with different Gaussian components remain
close when tranported, this is not the case of the transport on the right.

|  |
---|---|---
Figure 3: Left: two discrete distributions $\hat{\mu}$ (in gradient of colors)
and $\hat{\nu}$ (in blue) that have been drawn from two GMMs. The colors have
been added to $\hat{\mu}$ in order to visualize the couplings between
$\hat{\mu}$ and $\hat{\nu}$. Middle: transport of $\hat{\mu}$ obtained by
plugging the discrete plan that minimizes $MGW_{2}$ in (14), then using
$T_{GGW_{2}}$ with $\widetilde{\operatorname{Id}}_{2}=\operatorname{Id}_{2}$
for all components to transport the points. Right: transport of $\hat{\mu}$
obtained the same way as previously, but with another
$\widetilde{\operatorname{Id}}_{2}$.

Therefore, in order to design a transport plan $\pi_{\mathrm{MGW_{2}}}$
associated with the $MGW_{2}$ problem, it is in general necessary to determine
for each pair of indices $k,l$, which
$\widetilde{\operatorname{Id}}_{d^{\prime}}$ preserves, relatively to the
others, the global structure of the GMM, which becomes a difficult
combinatorial problem in itself as soon as the dimension ${d^{\prime}}$ is
large. Alternatively, a less tedious 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 $MGW_{2}$ problem. This is
the idea behind the embedded Wasserstein distance that we introduce in the
following section.

##  4 Embedded Wasserstein distance

In this section, we define an alternative distance to Gromov-Wasserstein also
invariant to isometries which explicits the isometric transformation applied
to one of the measure when computing the distance.

######  Definition .

Let $\mu\in\mathcal{P}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{P}(\mathbb{R}^{d^{\prime}})$. For $r\geq 1$ and $s\geq 1$, let
us denote $\mathrm{Isom}_{s}(\mathbb{R}^{r})$ the set of all isometries - for
the Euclidean norm - from $\mathbb{R}^{s}$ to $\mathbb{R}^{r}$. We define

$EW_{2}(\mu,\nu)=\inf\left\\{\inf_{\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})}W_{2}(\mu,\phi_{\\#}\nu),\inf_{\psi\in\mathrm{Isom}_{d}(\mathbb{R}^{d^{\prime}})}W_{2}(\psi_{\\#}\mu,\nu)\right\\}\;,$
($EW_{2}$)

with the convention that the infinimum over an empty set is equal to
$+\infty$.

Observe that if $d>{d^{\prime}}$, the set
$\mathrm{Isom}_{d}(\mathbb{R}^{d^{\prime}})$ is empty and so
$EW_{2}(\mu,\nu)=\inf_{\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})}W_{2}(\mu,\phi_{\\#}\nu)$.
In contrast, if $d<{d^{\prime}}$, $\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})$
is empty and so
$EW_{2}(\mu,\nu)=\inf_{\psi\in\mathrm{Isom}_{d}(\mathbb{R}^{d^{\prime}})}W_{2}(\psi_{\\#}\mu,\nu)$.
When $d=d^{\prime}$, the two infinimums are equivalent. In all what follows,
we will suppose without any loss of generality that $d\geq{d^{\prime}}$. More
generally, one can define, given two - not necessarily Euclidean - Polish
spaces $\mathcal{X}$ and $\mathcal{Y}$ each endowed with respective distances
$d_{\mathcal{X}}$ and $d_{\mathcal{Y}}$, and given two measures
$\mu\in\mathcal{P}(\mathcal{X})$ and $\nu\in\mathcal{P}(\mathcal{Y})$, the
following OT distance

$EW_{p}(\mu,\nu)=\inf\left\\{\inf_{\phi\in\mathrm{Isom}_{\mathcal{Y}}(\mathcal{X})}W_{p}(\mu,\phi_{\\#}\nu),\inf_{\psi\in\mathrm{Isom}_{\mathcal{X}}(\mathcal{Y})}W_{p}(\psi_{\\#}\mu,\nu)\right\\}\;,$
($EW_{p}$)

where $p\geq 1$ and $\mathrm{Isom}_{\mathcal{Y}}(\mathcal{X})$ (respectively
$\mathrm{Isom}_{\mathcal{X}}(\mathcal{Y}$)) is the set of all isometries from
$\mathcal{Y}$ to $\mathcal{X}$ (respectively from $\mathcal{X}$ to
$\mathcal{Y}$), i.e. such that
$d_{\mathcal{X}}(\phi(y),\phi(y^{\prime}))=d_{\mathcal{Y}}(y,y^{\prime})$ for
$y,y^{\prime}\in\mathcal{Y}$, and $W_{p}$ is the Wasserstein distance of order
$p$ on $\mathcal{P}(\mathcal{X})$. However there might be cases where both set
$\mathrm{Isom}_{\mathcal{Y}}(\mathcal{X})$ and
$\mathrm{Isom}_{\mathcal{X}}(\mathcal{Y})$ are empty. In that case, by
convention $EW_{p}(\mu,\nu)=+\infty$.

###  4.1 Equivalent formulations

First we show that when $d\geq{d^{\prime}}$, $EW_{2}$ is equivalent to the OT
distance introduced by Alvarez-Melis et al., (2019) described in Section 2.3
for a particular choice of transformation space $\mathcal{H}$. In all what
follows, we denote $\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ the _Stiefel
manifold_ James, (1976), i.e. the set of rectangular othogonal matrices of
size $d\times{d^{\prime}}$ such that $P^{T}P=\operatorname{Id}_{d^{\prime}}$.
We show the following result.

######  Proposition .

Let $\mu\in\mathcal{W}_{2}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{W}_{2}(\mathbb{R}^{d^{\prime}})$ and let suppose without any
loss of generality $d\geq{d^{\prime}}$. Then,

$EW^{2}_{2}(\mu,\nu)=\inf_{\pi\in\Pi(\mu,\nu)}\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d}),\
b\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-Py-b\|^{2}\mathrm{d}\pi(x,y)\;.$
(18)

###### Sketch of proof..

This is principally a consequence of the Mazur-Ulam theorem Mazur and Ulam,
(1932) which implies that any isometry from $\mathbb{R}^{d^{\prime}}$ to
$\mathbb{R}^{d}$ (endowed with the Euclidean norm) is necessarily affine. ∎

Note that it is easy to show from Section 4.1 that the infinimum in ($EW_{2}$)
is always achieved, see Appendix C for details. Now we can derive an
equivalent formulation of the $EW_{2}$ problem. For any measure
$\mu\in\mathcal{W}_{2}(\mathbb{R}^{d})$, let us denote $\bar{\mu}$ its
associated centered measure, i.e. the measure such that if $X\sim\mu$, then
$X-\mathbb{E}_{X\sim\mu}[X]\sim\bar{\mu}$.

######  Proposition .

Let $\mu\in\mathcal{W}_{2}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{W}_{2}(\mathbb{R}^{d^{\prime}})$ and let suppose
$d\geq{d^{\prime}}$. Problem ($EW_{2}$) is equivalent to

$\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\sup_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\mathrm{tr}(P^{T}K_{\pi})\;,$
(19)

where $\bar{\mu}$ and $\bar{\nu}$ denotes the centered measures associated
with $\mu$ and $\nu$ and
$K_{\pi}=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}\mathrm{d}\pi(x,y)$.
Furthermore, solving this latter problem in $P$ for a fixed $\pi$, (19)
reduces to

$\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\left\|\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}\mathrm{d}\pi(x,y)\right\|_{*}\;,$
(20)

and this is achieved at

$P_{\pi}^{*}=U_{\pi}\operatorname{Id}^{[d,{d^{\prime}}]}_{d^{\prime}}V_{\pi}^{T}\;,$
(21)

where $U_{\pi}\in\mathbb{O}(\mathbb{R}^{d})$ and
$V_{\pi}\in\mathbb{O}(\mathbb{R}^{d^{\prime}})$ are the left and right
orthogonal matrices associated with the Singular Value Decomposition (SVD) of
$K_{\pi}$.

###### Sketch of proof..

The equivalence between Problems (18) and (19) can be established with
straightforward computation. The equivalence between Problems (19) and (20)
can then be proved directly using (Alvarez-Melis et al.,, 2019, lemma 4.2).
For the sake of completeness, we offer an alternative proof of this latter
equivalence based on Lagragian analysis. ∎

Observe that $P_{\pi}^{*}$ is the projection of $K_{\pi}$ on the Stiefel
manifold $\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ since
$\mathrm{tr}(P^{T}K_{\pi})$ is the Frobenius inner-product between $P$ and
$K_{\pi}$ and so maximizing the inner-product is equivalent to minimizing the
distance between $K_{\pi}$ and $P$ since $\|P\|_{\mathcal{F}}$ is necessarily
equal to ${d^{\prime}}$. Using Propositions 4.1 and 4.1, one can show that
$EW_{2}$ also coincides with the OT distance between metric measure spaces
introduced by Sturm, (2006), that reads as,

$D((\mathcal{X},d_{\mathcal{X}},\mu),(\mathcal{Y},d_{\mathcal{Y}},\nu))=\inf_{\mathcal{Z},\psi,\phi}W_{2}(\psi_{\\#}\mu,\phi_{\\#}\nu)\;,$
(22)

where $(\mathcal{X},d_{\mathcal{X}},\mu)$ and
$(\mathcal{Y},d_{\mathcal{Y}},\nu)$ are two metric measure spaces as defined
in Section 2.2, $\mathcal{Z}$ is a third Polish space, and
$\psi:\mathcal{X}\rightarrow\mathcal{Z}$ and
$\phi:\mathcal{Y}\rightarrow\mathcal{Z}$ are two isometric mappings. Indeed,
we show the following result.

######  Proposition .

Let $\mu\in\mathcal{W}_{2}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{W}_{2}(\mathbb{R}^{d^{\prime}})$ with $d$ not necessarily
greater than ${d^{\prime}}$. Let $r\geq\max\\{d,{d^{\prime}}\\}$ and let
$\psi\in\mathrm{Isom}_{d}(\mathbb{R}^{r})$. Then,
$EW_{2}(\mu,\nu)=EW_{2}(\psi_{\\#}\mu,\nu)$.

###### Sketch of proof..

The proof of this result consists mainly in observing that transforming the
measures with an isometric mapping does not affect the singular values of the
cross-covariance matrix $\int xy^{T}\mathrm{d}\pi(x,y)$. ∎

Note that it is straightforward using Section 4.1 to show that $EW_{2}$
defines a pseudometric on $\bigsqcup_{k\geq 1}\mathcal{W}_{2}(\mathbb{R}^{k})$
that is invariant to isometries for the Euclidean norm, see Appendix C for
details. Finally, to complete this section, we emphasize that $EW_{2}$ is
different of the distance proposed in Cai and Lim, (2022), that we call here
$PW_{2}$ for Projection Wasserstein discrepancy. Between two measures
$\mu\in\mathcal{P}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{P}(\mathbb{R}^{d^{\prime}})$, and supposing that
$d\geq{d^{\prime}}$, they define

$PW_{2}(\mu,\nu)=\inf_{\phi\in\Gamma_{d}(\mathbb{R}^{d^{\prime}})}W_{2}(\phi_{\\#}\mu,\nu)\;,$
($PW_{2}$)

where $\Gamma_{d}(\mathbb{R}^{d^{\prime}})$ is the set of affine mappings
$\phi$ from $\mathbb{R}^{d}$ to $\mathbb{R}^{d^{\prime}}$ such that for all
$x\in\mathbb{R}^{d}$, $\phi(y)=P^{T}(x-b)$ where $b\in\mathbb{R}^{d}$ and
$P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$. Observe that in that case
$\phi$ is not an isometry since it is not injective. As a result, the term
only depending of the marginal $\mu$ in the developpement of the square of the
distance will now depends on $P$ and thus one can show - see Appendix D for
details - that Problem ($PW_{2}$) is equivalent to the following problem

$\inf_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\left(\mathrm{tr}(P^{T}\Sigma_{x}P)-2\mathrm{tr}(P^{T}K_{\pi})\right)\;,$
(23)

where
$\Sigma_{x}=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}xx^{T}d\bar{\mu}(x)$,
$K_{\pi}=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}d\pi(x,y)$.
Observe that Problem (23) can be interpreted as a regularization in $P$ of
Problem (19). It can also be interpreted as a $W_{2}$ problem between $\nu$
and a measure $\mu^{\prime}$ which has a different second-order moment than
$\mu$. More details about the difference between $PW_{2}$ and $EW_{2}$ ,
especially between Gaussian distributions, can be found in Appendix D.

###  4.2 Link with Gromov-Wasserstein

Before studying $EW_{2}$ between Gaussian mixture models, we exhibits here two
particuliar cases where the Problem ($EW_{2}$) shares close connections with
the Gromov-Wasserstein problem ($GW_{p}$) of order $2$ with inner-product
costs between the centered measures $\bar{\mu}$ and $\bar{\nu}$, i.e.

$GW_{2}^{2}(\langle.\rangle_{d},\langle.\rangle_{d^{\prime}},\bar{\mu},\bar{\nu})=\inf_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}(\langle
x,x^{\prime}\rangle_{d}-\langle
y,y^{\prime}\rangle_{d^{\prime}})^{2}\mathrm{d}\pi(x,y)\pi(x^{\prime},y^{\prime})\;.$
($GW_{2}$-IP)

These cases are when one of the two distributions is 1D and when both
distributions are Gaussian.

##### Univariate case.

When $\nu$ is a one-dimensional distribution on $\mathbb{R}$, the $EW_{2}$
problem is equivalent to the $GW_{2}$ problem with inner-products as cost
functions. Indeed, since
$K_{\pi}=\int_{\mathbb{R}^{d}\times\mathbb{R}}xy^{T}\mathrm{d}\pi(x,y)$ is of
size $d\times 1$, it has a unique singular value $\lambda_{\pi}>0$, and so one
can observe that Problems ($\mathcal{F}$-COV) and ($*$-COV), i.e.

$\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\left\|K_{\pi}\right\|_{\mathcal{F}}\qquad\text{and}\qquad\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\left\|K_{\pi}\right\|_{*}\;,$
(24)

that are respectively equivalent to Problems ($GW_{2}$-IP) and ($EW_{2}$), can
both be rewritten

$\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\lambda_{\pi}\;,$ (25)

implying thus that the two problems are equivalent. Observe that this allows
to recover almost directly the result of (Vayer,, 2020, Theorem 4.2.4) that
states that the GW problem with inner-product costs between two centered one-
dimensional distributions $\mu\in\mathcal{W}_{2}(\mathbb{R})$ and
$\nu\in\mathcal{W}_{2}(\mathbb{R})$ such that $\mu$ is absolutely continuous
with respect of the Lesbegue measure that are respectively of the form
$(\operatorname{Id}_{\mathbb{R}},T_{\mathrm{GW}}^{\shortuparrow})_{\\#}\mu$
and
$(\operatorname{Id}_{\mathbb{R}},T_{\mathrm{GW}}^{\shortdownarrow})_{\\#}\mu$
where $T_{\mathrm{GW}}^{\shortuparrow}=F^{-1}_{\nu}\circ
F^{\shortuparrow}_{\mu}$ and
$T_{\mathrm{GW}}^{\shortdownarrow}=F^{-1}_{\nu}\circ
F^{\shortdownarrow}_{\mu}$ with $F^{-1}_{\nu}$ being the generalized quantile
function associated with $\nu$ and $F^{\shortuparrow}_{\mu}$ and
$F^{\shortdownarrow}_{\mu}$ are respectively the _cumulative_ and _anti-
cumulative_ distribution function associated with $\mu$, i.e for all
$x\in\mathbb{R}$,

$F^{\shortuparrow}_{\mu}(x)=\mu(\left(-\infty,x\right])\qquad\text{and}\qquad
F^{\shortdownarrow}_{\mu}(x)=\mu(\left[-x,+\infty\right))\;.$ (26)

Interessingly, this result doesn’t hold anymore in the case the costs
functions in the GW problem are the squared Euclidean distances instead of the
inner-product. Indeed, in the case where are discrete distributions, Beinert
et al., (2022) have exhibited a counter-example. A procedure to exhibit other
counter-examples has been proposed in Dumont et al., (2022).

##### Gaussian distributions.

When $\mu=\mathrm{N}(m_{0},\Sigma_{0})$ with $m_{0}\in\mathbb{R}^{d}$ and
$\Sigma_{0}\in\mathbb{S}^{d}_{++}$ and $\nu=\mathrm{N}(m_{1},\Sigma_{1})$ with
$m_{1}\in\mathbb{R}^{d^{\prime}}$ and
$\Sigma_{1}\in\mathbb{S}^{d^{\prime}}_{+}$, one can show - see Section C.2 for
details - that Problems ($GW_{2}$-IP) and ($EW_{2}$) share some common
solutions which are of the ones exhibited in Salmona et al., (2021), i.e. with
$\pi$ of the form $(\operatorname{Id}_{d},T)$, with $T$ affine such that for
all $x\in\mathbb{R}^{d}$,

$T(x)=m_{1}+P_{1}\left(\widetilde{\operatorname{Id}}_{d^{\prime}}D_{1}^{\frac{1}{2}}{D_{0}^{({d^{\prime}})}}^{-\frac{1}{2}}\right)^{[{d^{\prime}},d]}P^{T}_{0}(x-m_{0})\;,$
(27)

where $(P_{0},D_{0})$ and $(P_{1},D_{1})$ are the respective diagonalizations
of $\Sigma_{0}\ (=P_{0}D_{0}P_{0}^{T})$ and $\Sigma_{1}\
(=P_{1}D_{1}P_{1}^{T})$ that sort the eigenvalues in non-increasing order and
$\widetilde{\operatorname{Id}}_{d^{\prime}}$ is any matrix of the form
$\mathrm{diag}((\pm 1)_{1\leq i\leq{d^{\prime}}})$. Note that it is not clear
if the two problems are strictly equivalent or only share some solutions,
since we are unable to prove that the solutions of the GW problem with inner-
product costs exhibited in Salmona et al., (2021) are the only solutions of
the problem. Thus, there might exist solutions of Problem ($GW_{2}$-IP) which
are not of the form above and which are not solution of Problem ($EW_{2}$).

###  4.3 Embedded Wasserstein distance between GMMs

Similarly to Delon and Desolneux, (2020), one can define an OT distance
derived from $EW_{2}$ when $\mu$ and $\nu$ are GMMs by restricting the set of
admissible couplings to be themselves GMMs.

######  Definition .

Let $\mu\in GMM_{K}(\mathbb{R}^{d})$ and $\nu\in
GMM_{L}(\mathbb{R}^{d^{\prime}})$ and suppose that $d\geq{d^{\prime}}$. We
define

$MEW_{2}(\mu,\nu)=\inf\left\\{\inf_{\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})}MW_{2}(\mu,\phi_{\\#}\nu),\inf_{\psi\in\mathrm{Isom}_{d}(\mathbb{R}^{d^{\prime}})}MW_{2}(\psi_{\\#}\mu,\nu)\right\\}\;.$
(28)

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\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ and $b\in\mathbb{R}^{d}$.
Similarly to $EW_{2}$, one can show that the infinimum in $\phi$ is always
achieved and that $MEW_{2}$ satisfies all the properties of a pseudometric on
$\mathcal{GMM}_{\infty}$. Supposing without any loss of generality that
$d\geq{d^{\prime}}$ and using the equivalent discrete formulation (6) of the
$MW_{2}$ problem, we get that for $\mu=\sum_{k}a_{k}\mu_{k}$ and
$\nu=\sum_{l}b_{l}\nu_{l}$, the problem is equivalent to

$\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\inf_{\omega\in\Pi(a,b)}\sum_{k,l}\omega_{k,l}W^{2}_{2}(\mu^{\prime}_{k},P_{\\#}\nu^{\prime}_{l})\;,$
($MEW_{2}$)

where for any $k\leq K$ and $l\leq L$, $\mu^{\prime}_{k}$ and
$\nu^{\prime}_{l}$ are the Gaussian components respectively associated to the
centered GMMs $\bar{\mu}$ and $\bar{\nu}$. Note that $\mu^{\prime}_{k}$ and
$\nu^{\prime}_{l}$ are not necessarily themselves centered.

####  4.3.1 Numerical solver

This time, it is not possible to derive analytically the closed form of the
optimal $P^{*}$ for Problem ($MEW_{2}$). However, one can still solve the
problem numerically using an alternate minimization scheme. Indeed, Problem
($MEW_{2}$) is not convex in $P$ and $\omega$, but is convex in $\omega$ if
$P$ is fixed and is furthermore a 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_{2}$) is not convex in $P$ for a
fixed $\omega$ because the feasible set, i.e. the Stiefel manifold
$\mathbb{V}_{{d^{\prime}}}(\mathbb{R}^{d})$, is not convex. For a fixed
$\omega$, the minimization in $P$ can be done by projected gradient descent
Calamai and Moré, (1987), i.e. for a given iterate $P^{\\{i\\}}$ and a given
$\omega$, the next iterate $P^{\\{i+1\\}}$ is given by

$P^{\\{i+1\\}}=\kappa_{\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\left(P^{\\{i\\}}-\eta\frac{\partial
J_{\omega}(P^{\\{i\\}})}{\partial P}\right)\;,$ (29)

where $\kappa_{\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}$ is the projection
mapping on the Stiefel manifold, where $\eta>0$ and where for all matrices $P$
of size ${d^{\prime}}\times d$,
$J_{\omega}(P)=\sum_{k,l}\omega_{k,l}W_{2}(\mu^{\prime}_{k},P_{\\#}\nu^{\prime}_{l})$.
As we have seen above in Section 4.1, for all $P$ of size ${d^{\prime}}\times
d$, the projection $\kappa_{\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}$ can be
written

$\kappa_{\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}(P)=U_{P}\operatorname{Id}_{d^{\prime}}^{[d,{d^{\prime}}]}V_{P}^{T}\;,$
(30)

where $U_{P}\in\mathbb{O}(\mathbb{R}^{d})$ and
$V_{P}\in\mathbb{O}(\mathbb{R}^{d^{\prime}})$ are respectively the left and
right orthogonal matrices associated with the SVD of $P$. In a nutshell, this
yields to Algorithm 1.

Algorithm 1 Mixture embedded Wasserstein solver

1:$\mu=\sum_{k}^{K}a_{k}\mu_{k}$, $\nu=\sum_{l}^{L}b_{l}\nu_{l}$,
$P^{\\{0\\}}\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$, $\eta>0$.

2:while not converged do

3: $[C]_{k,l}\leftarrow W^{2}_{2}(\mu_{k},\nu_{l})$ for $k=1,\dots,K$;
$l=1,\dots,L$

4: $\omega^{\\{i\\}}\leftarrow\textsc{Solve-OT}(a,b,C)$$\triangleright$ Solve
a classic OT problem.

5: while not converged do$\triangleright$ Do projected gradient descent on
$P$.

6: $A\leftarrow P^{\\{i-1\\}}-\eta\partial
J_{\omega^{\\{i\\}}}(P^{\\{i-1\\}})/\partial P$

7: $U,\Sigma,V^{T}\leftarrow\textsc{SVD}(A)$

8: $P^{\\{i\\}}\leftarrow
U\operatorname{Id}_{d^{\prime}}^{[d,{d^{\prime}}]}V^{T}$

9: end while

10:end while

11:return $\omega$, $P$

When $\mu$ and $\nu$ are only composed of non-degenerate Gaussian components,
one can compute $\partial J_{\omega}(P)/\partial P$ either by using automatic
differentiation Baydin et al., (2018) or by using the following technical
result, whose proof is postponed to Appendix B.

######  Lemma .

Let for any $k\leq K$, $\mu_{k}=\mathrm{N}(m_{0k},\Sigma_{0k})$ with
$m_{0k}\in\mathbb{R}^{d}$ and $\Sigma_{0k}\in\mathbb{S}^{d}_{++}$ and for any
$l\leq L$, $\nu_{l}=\mathrm{N}(m_{1l},\Sigma_{1l})$ with
$m_{1l}\in\mathbb{R}^{d^{\prime}}$ and
$\Sigma_{1l}\in\mathbb{S}^{d^{\prime}}_{++}$. For any $\omega$ in the $K\times
L$ simplex, let
$J_{\omega}:\mathbb{R}^{d\times{d^{\prime}}}\rightarrow\mathbb{R}$ be the
functional defined, for all matrix $P$ of size $d\times{d^{\prime}}$, by

$J_{\omega}(P)=\sum_{k,l}\omega_{k,l}W^{2}_{2}(\mu_{k},P_{\\#}\nu_{l})\;.$
(31)

Then for any full-rank matrix $P$ of size $d\times{d^{\prime}}$, we have

$\frac{\partial J_{\omega}(P)}{\partial
P}=2\sum_{k,l}\omega_{k,l}\left[Pm_{1l}m^{T}_{1l}-m_{0k}m^{T}_{1l}-\Sigma_{0k}P\Sigma_{1l}^{\frac{1}{2}}(\Sigma_{1l}^{\frac{1}{2}}P^{T}\Sigma_{0k}P\Sigma_{1l}^{\frac{1}{2}})^{-\frac{1}{2}}\Sigma_{1l}^{\frac{1}{2}}\right]\;.$
(32)

##### Initialization procedure.

Since the problem is non-convex, the solution to which Algorithm 1 converges
strongly depends on the initialization of $P$. It is therefore crucial to
design a good initialization procedure. To do so, we 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-regularized $W_{2}$ problem between the
two discrete measures $\mu^{\circ}=\sum_{k}a_{k}\delta_{m_{0k}}$ and
$\nu^{\circ}=\sum_{k}b_{l}\delta_{m_{1l}}$ with a large value of
regularization $\varepsilon_{0}$ in order to obtain a coupling
$\omega^{\\{1\\}}$. Then we set

$P^{\\{1\\}}=\textstyle{\kappa_{\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\left(\sum_{k,l}\omega^{\\{1\\}}_{k,l}m_{0k}m_{1l}^{T}\right)}\;.$
(33)

We then solve another entropic-regularized $W_{2}$ problem, this time between
$\mu^{\circ}$ and $P^{\\{1\\}}_{\\#}\nu^{\circ}$, using a smaller value of
regularization $\varepsilon_{1}=\alpha\times\varepsilon_{0}$ with
$\alpha\in(0,1)$. We obtain thus a new coupling $\omega^{\\{2\\}}$ and we can
then derive $P^{\\{2\\}}$ as previously. We repeat this procedure $N_{it}$
times until the regularization term $\varepsilon_{N_{it}}$ becomes small
enough. This boils down to Algorithm 2.

Algorithm 2 Annealed initialization procedure for mixture embedded Wasserstein

1:$a$, $b$, $\\{m_{0k}\\}_{k}^{K}$, $\\{m_{1l}\\}_{l}^{L}$,
$\varepsilon_{0}>0$, $\alpha\in(0,1)$,
$P^{\\{0\\}}=\operatorname{Id}_{d^{\prime}}^{[d,{d^{\prime}}]}$

2:for $i=1,\dots,N_{it}$ do

3: $[C]_{k,l}\leftarrow\|m_{0k}-P^{\\{i-1\\}}m_{1l}\|^{2}$

4:
$\omega^{\\{i\\}}\leftarrow\textsc{$\varepsilon$-OT}(a,b,C,\varepsilon_{i-1})$$\triangleright$
Solve a regularized OT problem.

5: $A\leftarrow\sum_{k,l}\omega^{\\{i\\}}_{k,l}m_{0k}m_{1l}^{T}$

6: $U,\Sigma,V^{T}\leftarrow\textsc{SVD}(A)$

7: $P^{\\{i\\}}=U\operatorname{Id}_{{d^{\prime}}}^{[d,{d^{\prime}}]}V^{T}$

8: $\varepsilon_{i}\leftarrow\alpha\varepsilon_{i-1}$$\triangleright$
Annealing scheme.

9:end for

10:return $P$

In practice, we set in all our experiments $\alpha=0.95$ and
$\varepsilon_{0}=1$ as in Alvarez-Melis et al., (2019). Furthermore we
observed that in most cases, setting $N_{it}=10$ was sufficient to obtain a
good initialization of $P$ for Algorithm 1.

####  4.3.2 Transportation plans and transportation maps

Since ($MEW_{2}$) has a continous equivalent formulation (28), one can derive
from any optimal solution $(\omega^{*},P^{*})$ of the former, an optimal
solution $(\pi^{*},\phi^{*})$ of the latter. More precisely, we have on the
one hand for all $y\in\mathbb{R}^{d^{\prime}}$, $\phi^{*}(y)=P^{*}y+b^{*}$,
where $b^{*}=\mathbb{E}_{X\sim\mu}[X]-P^{*}\mathbb{E}_{Y\sim\nu}[Y]$, and on
the other hand for all $(x,y)\in\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}$,

$\pi^{*}(x,y)=\sum_{k,l}\omega^{*}_{k,l}p_{\mu_{k}}(x)\delta_{y=\psi^{*}\circ
T^{k,l}_{W_{2}}(x)}\;,$ (34)

where $T^{k,l}_{W_{2}}(x)$ is the optimal $W_{2}$ transport map between
$\mu^{\prime}_{k}$ and $P^{*}_{\\#}\nu^{\prime}_{l}$ and
$\psi^{*}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d^{\prime}}$ is defined for all
$x\in\mathbb{R}^{d}$ as $\psi^{*}(x)=P^{*T}(x-b^{*})$. As in Delon and
Desolneux, (2020), it is possible to define a unique assignment of each $x$ by
setting for all $x\in\mathbb{R}^{d}$,

$T_{\mathrm{mean}}(x)=\mathbb{E}_{(X,Y)\sim\pi^{*}}[Y|X=x]=\textstyle{\frac{\sum_{k,l}\omega^{*}_{k,l}p_{\mu_{k}}(x)\psi^{*}\circ
T^{k,l}_{W_{2}}(x)}{\sum_{k}a_{k}p_{\mu_{k}}(x)p_{\mu_{k}}(x)}}\;.$ (35)

Note that $T_{\mathrm{mean}}$ is not a Monge map since $\pi^{*}$ is not of the
form $(\operatorname{Id}_{d},T)_{\\#}\mu$. In particular,
$T_{\mathrm{mean}\\#}\mu$ is not equal to $\nu$ and $T_{\mathrm{mean}\\#}\mu$
is not necessarily the gradient of a convex function. Another possible way to
define an assignment proposed by Delon and Desolneux, (2020) is to define it
as a random assignment for a fixed $x$, i.e.

$T_{\mathrm{rand}}(x)=\psi^{*}\circ T^{k,l}_{W_{2}}(x)\qquad\text{with
probability}\textstyle{\qquad
p_{k,l}(x)=\frac{\omega_{k,l}^{*}p_{\mu_{k}}(x)}{\sum_{i}a_{i}p_{\mu_{i}}(x)}}\;.$
(36)

Hence, when using $MEW_{2}$ to obtain an assignment between two sets
$\\{x_{i}\\}^{M}_{i}$ and $\\{y_{j}\\}^{N}_{j}$ of respectively $M$ and $N$
vectors of $\mathbb{R}^{d}$ and $\mathbb{R}^{d^{\prime}}$, one can compute
either $T_{\mathrm{mean}}(x)$ or $T_{\mathrm{rand}}(x)$ for each $x_{i}$, and
then determine for each $x_{i}$ which $y_{j}$ is the closest of
$T_{\mathrm{mean}}(x_{i})$ \- or $T_{\mathrm{rand}}(x_{i})$ \- using a
nearest-neighbor algorithm Fix and Hodges, (1951).

####  4.3.3 Improving the $MGW_{2}$ method

Inspired by the $MEW_{2}$ method presented above, we propose in this section
to improve the $MGW_{2}$ 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_{2}$ similarly to
(34).

##### Annealing scheme.

Since Problem ($MGW_{2}$) is non-convex, we are only guaranteed to converge
towards a local minimum when solving it with a classic non-regularized GW
solver Peyré et al., (2016). Furthermore, the convergence towards a particular
minimum depends strongly on the initialization of the coupling $\omega$. Since
the discrete GW problem in $MGW_{2}$ is of very small scale and so not costly
in itself, we propose, by anology with $MEW_{2}$, 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.

Algorithm 3 Annealed mixture Gromov-Wasserstein solver

1:$\mu=\sum_{k}^{K}a_{k}\mu_{k}$, $\nu=\sum_{l}^{L}b_{l}\nu_{l}$,
$\alpha\in(0,1)$, $\varepsilon_{0}$, $\omega^{\\{0\\}}=ab^{T}$

2:$[C^{x}]_{k,i}\leftarrow W_{2}^{2}(\mu_{k},\mu_{i})$ for $k=1,\dots,K$,
$i=1,\dots,K$

3:$[C^{y}]_{l,j}\leftarrow W_{2}^{2}(\nu_{l},\nu_{j})$ for $l=1,\dots,L$,
$j=1,\dots,L$

4:for $n=1,\dots,N_{it}$ do

5:
$\omega^{\\{n\\}}\leftarrow\textsc{$\varepsilon$-GW}(a,b,C^{x},C^{y},\varepsilon_{n-1},\omega^{\\{n-1\\}})$$\triangleright$
Solve a regularized GW problem.

6: $\varepsilon_{n}\leftarrow\alpha\varepsilon_{n-1}$$\triangleright$
Annealing scheme.

7:end for

8:return $\textsc{Solve-
GW}(a,b,C^{x},C^{y},\omega^{\\{N_{it}\\}})$$\triangleright$ Solve the non-
regularized GW problem.

As previously, we set in our experiments $\alpha=0.95$ and $\varepsilon_{0}=1$
as in Alvarez-Melis et al., (2019) and we observed that, in toy cases where we
know what the global minimum is, that $N_{it}=10$ seemed to be a sufficient
number of iterations to prevent the algorithm from converging towards a sub-
optimal minimum.

##### Designing a transportation plan.

Still by analogy with $MEW_{2}$, one can design a transportation plan for
$MGW_{2}$ by defining a matrix
$P_{MGW_{2}}\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ and a vector
$b_{MGW_{2}}\in\mathbb{R}^{d}$, and then replacing $T_{W_{2}}\circ\psi^{*}$ in
(34) by $T_{W_{2}}\circ\psi_{MGW_{2}}$, where for all $x\in\mathbb{R}^{d}$,
$\psi_{MGW_{2}}(x)=P_{MGW_{2}}^{T}(x-b_{MGW_{2}})$. Given two GMMs
$\mu=\sum_{k}a_{k}\mu_{k}$ and $\nu=\sum_{l}b_{l}\nu_{l}$ respectively in
$GMM_{K}(\mathbb{R}^{d})$ and $GMM_{L}(\mathbb{R}^{d^{\prime}})$ and given the
optimal discrete plan $\omega^{*}$ solution of Problem ($MGW_{2}$), one can
define the matrix $P_{MGW_{2}}$ as the solution of the following problem

$\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\sum_{k,l}\omega^{*}_{k,l}W_{2}^{2}(\mu^{\prime}_{k},P_{\\#}\nu^{\prime}_{l})\;,$
(37)

where $\mu^{\prime}_{k}$ and $\nu^{\prime}_{l}$ are the Gaussian component of
the centered GMMs $\bar{\mu}$ and $\bar{\nu}$, then we can set
$b_{MGW_{2}}=\mathbb{E}_{X\sim\mu}[X]-P_{MGW_{2}}\mathbb{E}_{Y\sim\mu}[Y]$. As
above, this problem can be solved numerically by performing a projected
gradient descent on $P$, using either automatic differentiation or Section
4.3.1. This is also a non-convex optimization problem since
$\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ 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^{\\{0\\}}_{MGW_{2}}=\textstyle{\kappa_{\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\left(\sum_{k,l}\omega^{*}_{k,l}m_{0k}m^{T}_{1l}\right)}\;.$
(38)

Finally, using $P_{MGW_{2}}$ one can define a continous plan $\pi_{MGW_{2}}$
associated with the discrete optimal plan $\omega^{*}$ solution of the
$MGW_{2}$ problem similarly to (34). We can therefore use $MGW_{2}$ to
transport distributions, using as previously either $T_{\mathrm{mean}}$ or
$T_{\mathrm{rand}}$. We can also, as for $MEW_{2}$, use $MGW_{2}$ to obtain an
assignment between two sets of points.

##  5 Experiments

In what follows, we use $MGW_{2}$ and $MEW_{2}$ to solve Gromov-Wasserstein
related tasks on various datasets. More precisely, we apply first the two
methods on simple toy low-dimensional GMMs. Then, we show that both methods
can be used to solve relatively efficiently GW related tasks on real datasets
in large scale settings involving sometimes several tens of thousands of
points. We apply thus our methods to two shape matching problems, then to
color transfer on hyperspectral images. In all our experiments, we use the
numerical solvers provided by the Python Optimal Transport (POT) package555The
package is accessible here: https://pythonot.github.io/. Flamary et al.,
(2021) that implements solvers for the non-regularized and regularized classic
OT and GW problems.

###  5.1 Low dimensional GMMs

In Figure 4, we use again the example of Figure 3 and we derive an optimal
transport plan for the $MGW_{2}$ problem as described in Section 4.3.3. We
also show the plan obtained by solving the $EW_{2}$ problem. One can see 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.

data |  $\bm{MGW_{2}}$ |  $\bm{MEW_{2}}$
---|---|---

|  |
---|---|---

Figure 4: Left: two discrete distributions $\hat{\mu}$ (in gradient of colors)
and $\hat{\nu}$ (in blue) that have been drawn from two GMMs. The colors have
been added to $\hat{\mu}$ in order to visualize the couplings between
$\hat{\mu}$ and $\hat{\nu}$. Middle: transport of $\hat{\mu}$ obtained by
solving the $MGW_{2}$ problem, then deriving
$P_{MGW_{2}}\in\mathbb{V}_{2}(\mathbb{R}^{2})$ by solving Problem (37). Right:
transport of $\hat{\mu}$ obtained by solving the $MEW_{2}$ problem.

###  5.2 Application to shape matching

We apply now our methods on two shape matching problems. The first one
consists in reproducing an experiment originally conducted in Rustamov et al.,
(2013) and presented in Solomon et al., (2016) with the use of entropic-
regularized GW, that aims to revover the cyclical nature of a horse’s gallop.
The second problem consists in drawing correspondences betwen human shaped
meshes from the SHREC’19 dataset666The SHREC’19 dataset is accessible here:
http://profs.scienze.univr.it/ marin/shrec19/ Melzi et al., (2019) in the
sense that we aim to assign a hand with a hand, a foot with a foot, etc. Note
that the goal of these experiments is not to obtain state-of-the-art results
in shape-matching, but rather to demonstrate the usability of our methods in
moderate-to-large scale settings.

##### 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 et
al., (2016) with the use of entropic-regularized GW. The goal of this
experiment is to compute a matrix of pairwise distances between the 45 meshes
representing a galloping horse, and then to conduct a Multi-Dimensional
Scaling (MDS) Borg and Groenen, (2005) \- which roughly can be thought as a
generalization of PCA - of the pairwise distances in order to plot each mesh
as a $2$-dimensional point. The results can be found in Figure 5. As in
Solomon et al., (2016), the cyclical nature of the motion is recovered in both
cases when $MGW_{2}$ or $MEW_{2}$ is used to compare the meshes. Each mesh is
composed of approximately $9000$ vertices and the average time to compute one
distance when using the POT implementation of the entropic-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
$MGW_{2}$, and around one hour using $MEW_{2}$, these times including the
fitting of all the GMMs.

$\bm{MGW_{2}}$ |  $\bm{MEW_{2}}$ |  Meshes
---|---|---

$\begin{array}[]{cc}\includegraphics[width=212.47617pt]{Figures/MDS_gromov_20.pdf}&\includegraphics[width=212.47617pt]{Figures/MDS_ew2_20.pdf}\end{array}$

$\begin{array}[]{cc}\includegraphics[width=208.13574pt]{Figures/cheval_0.png}&\includegraphics[width=208.13574pt]{Figures/cheval_4.png}\\\
\includegraphics[width=208.13574pt]{Figures/cheval_7.png}&\includegraphics[width=208.13574pt]{Figures/cheval_10.png}\end{array}$

Figure 5: MDS on the galloping horse animation using the $MGW_{2}$ distance
(left), and the $MEW_{2}$ distance (middle). Each point corresponds to a given
mesh and the meshes are colored in function of their number in the sequence.
The computation of both distances have been done by first fitting GMMs with
$20$ components on each mesh independently.

##### Local minima.

To highlight the importance of using an annealing scheme when deriving
$MGW_{2}$ or $MEW_{2}$, 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 $MGW_{2}^{2}$ and
$MEW_{2}^{2}$ between one given fixed mesh and all the others. In both cases,
the annealing scheme seems to be useful to prevent the solver to converge
towards sub-optimal mininima. However, if the $MGW_{2}$ solver seems to often
converge to the same optimum regardless the use of the annealing scheme, this
is not the case of $MEW_{2}$ which, without the annealing initialization
procedure (Algorithm 2), converges most of the time to a sub-optimal minimum,
so much that the periodical aspect doesn’t even appear in that case. Beside to
highlight the importance of using a good initialization, this experiment also
emphasizes the fact that when solving a 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-optimal.

|
---|---

Figure 6: Left: Evolution of $MGW^{2}_{2}$ 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 $MEW^{2}_{2}$ 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 have been done by first fitting GMMs with $20$ components on each
mesh independently.

##### Matching human shaped meshes.

To demonstrate the usability of our methods in larger scale settings, we use
the SHREC’19 dataset 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 represents 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 vertices. The
results are displayed on Figure 7. Observe that in most cases, both methods
seem to be able to match correctly the colored parts. Yet in the last row,
$MEW_{2}$ matches a leg at the left in red to an arm at the right. This
probably implies 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
$MGW_{2}$ or $MEW_{2}$ fail to find correct correspondences and exhibit
behaviors similar to $MEW_{2}$ in the last row, which suggests that the
methods converge sometimes to sub-optimal minima despite the annealing
schemes.

Source | Target ($\bm{MGW_{2}}$) | Target ($\bm{MEW_{2}}$)
---|---|---
|  |
|  |
|  |

Figure 7: Shape matching between human-shaped meshes using $MGW_{2}$ (middle)
and $MEW_{2}$ (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 red correspond
to Gaussian components that are matched together when solving $MGW_{2}$ or
$MEW_{2}$. From left to right and top to bottom, the meshes are composed
respectively of $84912$, $30300$, $75000$, $273624$, $360678$, and $360357$
vertices.

###  5.3 Application to hyperspectral image color transfer

The goal here is to reproduce the experiment of color transfer conducted in
Delon and Desolneux, (2020), but this time using a hyperspectral image, i.e an
image with more than $3$ color channels. More precisely, we aim to create an
RGB image from an hyperspectral image $u$ using the colors of another RGB
image $v$. To do so, we consider images as empirical distributions in the
color spaces and we solve a Gromov-Wasserstein problem between the
distributions $\hat{\mu}=\frac{1}{M}\sum_{k}^{M}\delta_{u_{k}}$ and
$\hat{\nu}=\frac{1}{N}\sum_{l}^{N}\delta_{v_{l}}$, where $M$ and $N$ are the
number of pixels in respectively the hyperspectral image and the RGB image we
use as color palette, and $\\{u_{k}\\}_{k}^{M}$ and $\\{v_{l}\\}_{l}^{N}$ are
the values at each pixel, i.e for here all $l$, $v_{l}\in\mathbb{R}^{3}$ and
for all $k$, $u_{k}\in\mathbb{R}^{d}$ with $d>3$. We thus fit two GMMs $\mu$
and $\nu$ and respectively $\hat{\mu}$ and $\hat{\nu}$ and we use $MGW_{2}$ or
$MEW_{2}$ to derive a mapping
$T_{\mathrm{mean}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{3}$, as described in
Section 4.3.2. We apply this process to a hyperspectral image of $512\times
512$ pixels with $15$ channels that are displayed in Figure 8 top left. We use
as color palettes two paintings by Gauguin and Renoir, displayed in Figure 8
top right, that are respectively _Manhana no atua_ (top) and _Le déjeuner des
canotiers_ (bottom). These two images are composed of $1024\times 768$ pixels.
The resulting images $T_{\mathrm{mean}}(u)$ are displayed in Figure 8 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 $T_{\mathrm{mean}}$. This experiment
shows that $MGW_{2}$ and $MEW_{2}$ can be used in large scale settings:
observe indeed that the color distributions $\hat{\mu}$ and $\hat{\nu}$ are
composed respectively of approximatively $300000$ and $800000$ points, which
makes the problem intractable with entropic-GW solvers such as Peyré et al.,
(2016) or Solomon et al., (2016). Furthermore, note also that $d=15$ is
already a relatively high dimension in the context of Gromov-Wasserstein. In
term of computation time, the fitting of the GMM for the hyperspectral image
takes aproximatively $15$ minutes against one minute 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 $MEW_{2}$. Finally, for both methods, it takes around
$20$ minutes to compute the whole RGB image $T_{\mathrm{mean}}(u)$.

---

|  |  |
---|---|---|---

$\bm{MGW_{2}}$ |  $\bm{MEW_{2}}$
---|---

Figure 8: Color transfers between a hyperspectral image with $15$ channels
(top left) and two paintings by Gauguin and Renoir (top right, top to bottom).
Bottom line: the obtained RGB images using $MGW_{2}$ and $MEW_{2}$. 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)).

##  6 Discussion

In this chapter, we have introduced two new OT distances on the set of
Gaussian mixture models, $MGW_{2}$ and $MEW_{2}$, and we have shown that they
both can be used to solve relatively efficiently Gromov-Wasserstein related
problems on Euclidean spaces, especially in moderate-to-large scale settings
involving several tens of thousands of points. These OT distances are also by
design particularly suited to settings where there already exists a kind of
clustering structure in the data. This being said, if $MEW_{2}$ 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 $MGW_{2}$ for roughly the same quality
of results, and so we believe that $MGW_{2}$ is a better choice in practice.
This latter distance is part of the families of Gromov-Wasserstein type OT
distances that reduce the size of the GW problem by quantization Chowdhury et
al., (2021) or by clustering Xu et al., (2019); Blumberg et al., (2020). To
the best of our knowledge, however, no such methods specific to the Euclidean
case had already been proposed in the literature. $MGW_{2}$ could also 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 applied problems, an interesting extension on our work
could be to 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 Expectation-Maximization (EM) algorithm
Dempster et al., (1977) which can become relatively 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 algorithms that seems to outperform it. These approaches
are based on Riemannian stochastic optimization, leveraging the rich
Riemmanian structure of the set of positive definite matrices. Another
interesting alternative that has been shown to outperform the EM algorithm has
been proposed by Kolouri et al., (2018) and is based on the minimization of
the sliced-Wasserstein distance. Integrating this in our method could result
thus in an approach fully-based on optimal transport.

Another possible limitation of our work lies in the fact that the $MGW_{2}$
solver converges sometimes to sub-optimal local minima. If the annealed
procedure introduced in Section 4.3.3 seems to reduce this issue, we generally
have no guarantee that the solution we have converged to is optimal. This is
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 a sub-optimal local minimum because we
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 Euclidean distances. A future perspective
of work could be therefore to study if a similar idea could be applied for
solving the $MGW_{2}$ problem.

## References

  * Altschuler et al., (2019) 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 et al., (2018) Altschuler, J., Bach, F., Rudi, A., and Weed, J. (2018).  Approximating the quadratic transportation metric in near-linear time.  arXiv preprint arXiv:1810.10046.
  * Alvarez-Melis and Jaakkola, (2018) Alvarez-Melis, D. and Jaakkola, T. (2018).  Gromov–Wasserstein alignment of word embedding spaces.  In Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing, pages 1881–1890.
  * Alvarez-Melis et al., (2019) 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 Statistics, pages 1870–1879. PMLR.
  * Arjovsky et al., (2017) Arjovsky, M., Chintala, S., and Bottou, L. (2017).  Wasserstein generative adversarial networks.  In International Conference on Machine Learning, pages 214–223. PMLR.
  * Baker, (1971) Baker, J. (1971).  Isometries in normed spaces.  The American Mathematical Monthly, 78(6):655–658.
  * Baydin et al., (2018) Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M. (2018).  Automatic differentiation in machine learning: a survey.  Journal of Marchine Learning Research, 18:1–43.
  * Beinert et al., (2022) Beinert, R., Heiss, C., and Steidl, G. (2022).  On assignment problems related to Gromov–Wasserstein distances on the real line.  arXiv preprint arXiv:2205.09006.
  * Blumberg et al., (2020) Blumberg, A. J., Carriere, M., Mandell, M. A., Rabadan, R., and Villar, S. (2020).  MREC: a fast and versatile framework for aligning and matching point clouds with applications to single cell molecular data.  stat, 1050:20.
  * Bolley, (2008) Bolley, F. (2008).  Separability and completeness for the Wasserstein distance.  Lecture Notes in Mathematics-Springer-Verlag-, 1934:371.
  * Borg and Groenen, (2005) Borg, I. and Groenen, P. J. (2005).  Modern multidimensional scaling: theory and applications.  Springer Science & Business Media.
  * Brogat-Motte et al., (2022) Brogat-Motte, L., Flamary, R., Brouard, C., Rousu, J., and d’Alché Buc, F. (2022).  Learning to predict graphs with fused Gromov–Wasserstein barycenters.  In International Conference on Machine Learning, pages 2321–2335. PMLR.
  * Bunne et al., (2019) 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 and Lim, (2022) Cai, Y. and Lim, L.-H. (2022).  Distances between probability distributions of different dimensions.  IEEE Transactions on Information Theory, 68(6):4020–4031.
  * Calamai and Moré, (1987) Calamai, P. H. and Moré, J. J. (1987).  Projected gradient methods for linearly constrained problems.  Mathematical programming, 39(1):93–116.
  * Chen et al., (2018) Chen, Y., Georgiou, T. T., and Tannenbaum, A. (2018).  Optimal transport for Gaussian mixture models.  IEEE Access, 7:6269–6278.
  * Chowdhury and Mémoli, (2019) Chowdhury, S. and Mémoli, F. (2019).  The Gromov–Wasserstein distance between networks and stable network invariants.  Information and Inference: A Journal of the IMA, 8(4):757–787.
  * Chowdhury et al., (2021) Chowdhury, S., Miller, D., and Needham, T. (2021).  Quantized Gromov–Wasserstein.  In Machine Learning and Knowledge Discovery in Databases. Research Track: European Conference, ECML PKDD 2021, Bilbao, Spain, September 13–17, 2021, Proceedings, Part III 21, pages 811–827. Springer.
  * Cohen and Guibasm, (1999) 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, (1920) Courant, R. (1920).  Über die eigenwerte bei den differentialgleichungen der mathematischen physik.  Mathematische Zeitschrift, 7(1-4):1–57.
  * Courty et al., (2018) Courty, N., Flamary, R., and Ducoffe, M. (2018).  Learning Wasserstein embeddings.  In ICLR 2018-6th International Conference on Learning Representations, pages 1–13.
  * Courty et al., (2016) 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. IEEE.
  * Cuturi, (2013) Cuturi, M. (2013).  Sinkhorn distances: lightspeed computation of optimal transport.  Advances in neural information processing systems, 26.
  * Delon and Desolneux, (2020) 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–970.
  * Delon et al., (2022) Delon, J., Gozlan, N., and Saint-Dizier, A. (2022).  Generalized Wasserstein barycenters between probability measures living on different subspaces.  Annals of Applied Probability.
  * Dempster et al., (1977) Dempster, A. P., Laird, N. 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.
  * Dumont et al., (2022) Dumont, T., Lacombe, T., and Vialard, F.-X. (2022).  On the existence of Monge maps for the Gromov–Wasserstein problem.
  * Fischer, (1905) Fischer, E. (1905).  Über quadratische formen mit reellen koeffizienten.  Monatshefte für Mathematik und Physik, 16:234–249.
  * Fix and Hodges, (1951) Fix, E. and Hodges, J. (1951).  Discriminatory analysis: nonparametric discrimination: consistency properties. report. 4.  T. USAF School of Aviation Medicine.
  * Flamary et al., (2021) 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, 22(1):3571–3578.
  * Forrow et al., (2019) 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 2454–2465. PMLR.
  * Genevay et al., (2018) 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 et al., (1984) Givens, C. R., Shortt, R. M., et al. (1984).  A class of Wasserstein metrics for probability distributions.  In Michigan Mathematical Journal, volume 31, pages 231–240. the University of Michigan.
  * Hosseini and Sra, (2020) Hosseini, R. and Sra, S. (2020).  An alternative to EM for Gaussian mixture models: batch and stochastic Riemannian optimization.  Mathematical programming, 181(1):187–223.
  * James, (1976) James, I. M. (1976).  The topology of Stiefel manifolds, volume 24.  Cambridge University Press.
  * Kloeckner, (2010) Kloeckner, B. (2010).  A geometric study of Wasserstein spaces: Euclidean spaces.  Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 9(2):297–323.
  * Kolouri et al., (2019) Kolouri, S., Nadjahi, K., Simsekli, U., Badeau, R., and Rohde, G. (2019).  Generalized sliced Wasserstein distances.  Advances in neural information processing systems, 32.
  * Kolouri et al., (2018) Kolouri, S., Rohde, G. K., and Hoffmann, H. (2018).  Sliced Wasserstein distance for learning Gaussian mixture models.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3427–3436.
  * Lambert et al., (2022) 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 et al., (2022) Leclaire, A., Delon, J., and Desolneux, A. (2022).  Optimal transport between GMMs for texture synthesis.
  * Luzi et al., (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 on Applications of Computer Vision, pages 279–288.
  * Magnus and Neudecker, (2019) Magnus, J. R. and Neudecker, H. (2019).  Matrix differential calculus with applications in statistics and econometrics.  John Wiley & Sons.
  * Mazur and Ulam, (1932) Mazur, S. and Ulam, S. (1932).  Sur les transformations isométriques d’espaces vectoriels normés.  CR Acad. Sci. Paris, 194(946-948):116.
  * Melzi et al., (2019) 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 Eurographics Association.
  * Mémoli, (2009) Mémoli, F. (2009).  Spectral Gromov–Wasserstein distances for shape matching.  In 2009 IEEE 12th International Conference on Computer Vision Workshops, ICCV Workshops, pages 256–263. IEEE.
  * Mémoli, (2011) Mémoli, F. (2011).  Gromov–Wasserstein distances and the metric approach to object matching.  In Foundations of Computational Mathematics, volume 11, pages 417–487. Springer.
  * Pasande et al., (2022) Pasande, M., Hosseini, R., and Araabi, B. N. (2022).  Stochastic first-order learning for large-scale flexibly tied Gaussian mixture model.  arXiv preprint arXiv:2212.05402.
  * Pedregosa et al., (2011) 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, 12:2825–2830.
  * Pele and Taskar, (2013) 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 et al., (2008) Petersen, K. B., Pedersen, M. S., et al. (2008).  The matrix cookbook.  Technical University of Denmark, 7(15):510.
  * Peyré et al., (2016) Peyré, G., Cuturi, M., and Solomon, J. (2016).  Gromov–Wasserstein averaging of kernel and distance matrices.  In International conference on machine learning, pages 2664–2672. PMLR.
  * Rabin et al., (2014) 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 et al., (2012) 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, Ein-Gedi, Israel, May 29–June 2, 2011, Revised Selected Papers 3, pages 435–446. Springer.
  * Rustamov et al., (2013) 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 et al., (2023) Ryner, M., Kronqvist, J., and Karlsson, J. (2023).  Globally solving the Gromov–Wasserstein problem for point clouds in low dimensional euclidean spaces.  arXiv preprint arXiv:2307.09057.
  * Salmona et al., (2021) Salmona, A., Delon, J., and Desolneux, A. (2021).  Gromov–Wasserstein distances between Gaussian distributions.  Journal of Applied Probability, 59(4).
  * Scetbon and Cuturi, (2020) Scetbon, M. and Cuturi, M. (2020).  Linear time Sinkhorn divergences using positive features.  Advances in Neural Information Processing Systems, 33:13468–13480.
  * Scetbon et al., (2021) Scetbon, M., Cuturi, M., and Peyré, G. (2021).  Low-rank Sinkhorn factorization.  In International Conference on Machine Learning, pages 9344–9354. PMLR.
  * Scetbon et al., (2022) Scetbon, M., Peyré, G., and Cuturi, M. (2022).  Linear-time Gromov–Wasserstein distances using low rank couplings and costs.  In International Conference on Machine Learning, pages 19347–19365. PMLR.
  * Seguy et al., (2017) Seguy, V., Damodaran, B. B., Flamary, R., Courty, N., Rolet, A., and Blondel, M. (2017).  Large-scale optimal transport and mapping estimation.  arXiv preprint arXiv:1711.02283.
  * Sembach et al., (2022) 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 and Knopp, (1967) Sinkhorn, R. and Knopp, P. (1967).  Concerning nonnegative matrices and doubly stochastic matrices.  Pacific Journal of Mathematics, 21(2):343–348.
  * Solomon et al., (2015) 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 et al., (2016) Solomon, J., Peyré, G., Kim, V. G., and Sra, S. (2016).  Entropic metric alignment for correspondence problems.  ACM Transactions on Graphics (ToG), 35(4):1–13.
  * Sturm, (2006) Sturm, K.-T. (2006).  On the geometry of metric measure spaces. i.  Acta Math, 196:65–131.
  * Sturm, (2012) 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, (2010) Takatsu, A. (2010).  On Wasserstein geometry of Gaussian measures.  In Probabilistic approach to geometry, pages 463–472. Mathematical Society of Japan.
  * Tolstikhin et al., (2018) 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, (2020) Vayer, T. (2020).  A contribution to Optimal Transport on incomparable spaces.  PhD thesis, Lorient.
  * (70) 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.
  * (71) Vayer, T., Flamary, R., Courty, N., Tavenard, R., and Chapel, L. (2019b).  Sliced Gromov–Wasserstein.  Advances in Neural Information Processing Systems, 32.
  * Villani, (2008) Villani, C. (2008).  Optimal transport: old and new, volume 338.  Springer Science & Business Media.
  * Xu et al., (2019) Xu, H., Luo, D., and Carin, L. (2019).  Scalable Gromov–Wasserstein learning for graph partitioning and matching.  Advances in neural information processing systems, 32.
  * Xu et al., (2018) 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 and Spragins, (1968) Yakowitz, S. J. and Spragins, J. D. (1968).  On the identifiability of finite mixtures.  The Annals of Mathematical Statistics, 39(1):209–214.
  * Zhang and Chen, (2020) Zhang, Q. and Chen, J. (2020).  A unified framework for Gaussian mixture reduction with composite transportation distance.  arXiv preprint arXiv:2002.08410.

## Organization of the supplementary

The supplementary is organized as follows. First, in Appendix A, we show five
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. In
Appendix C, we give more details on the metric properties of $EW_{2}$ and on
the behavior of this latter distance on Gaussian distributions. Finally, in
Appendix D, we give more details on the OT distance introduced by Cai and Lim,
(2022) that we call here projection Wasserstein discrepancy.

##  Appendix A Technical lemmas

Before turning to the proofs of the theoretical results, we state here five
technical lemmas that will be used throughout the proofs of the results of the
paper.

###  A.1 A property of couplings between measures living in different
dimensions

First we start by recalling the following result by Delon et al., (2022).

######  Lemma A0 Delon et al., (2022).

Let $\mu\in\mathcal{W}_{2}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{W}_{2}(\mathbb{R}^{d^{\prime}})$ with $d$ not necessarily
greater than ${d^{\prime}}$, and let
$T:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ be a measurable map. Then
$\pi^{\prime}\in\Pi(\mu,T_{\\#}\nu)$ if and only if there is some
$\pi\in\Pi(\mu,\nu)$ such that
$\pi^{\prime}=(\operatorname{Id}_{d},T)_{\\#}\pi$. In particular, if there
exist $a,b\geq 0$ such that $\|T(y)\|\leq a+b\|y\|$ for all
$y\in\mathbb{R}^{d^{\prime}}$, then

$\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-T(y)\|^{2}\mathrm{d}\pi(x,y)=\inf_{\pi\in\Pi(\mu,T_{\\#}\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x-z\|^{2}\mathrm{d}\pi(x,z)\;.$
(39)

###  A.2 Isometries in Euclidean spaces

We show the following result, that states that any isometry
$T:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ for the Euclidean norms
is affine and of the form, for all $y\in\mathbb{R}^{d^{\prime}}$, $T(y)=Py+b$,
where $b\in\mathbb{R}^{d}$ and $P$ is in the Stiefel manifold
$\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ as defined in Section 4.1.

######  Lemma A0.

Suppose $d\geq{d^{\prime}}$. Then
$\phi:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ is an isometry for the
Euclidean norm if and only if there exists
$P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ and $b\in\mathbb{R}^{d}$ such
that for all $y\in\mathbb{R}^{d^{\prime}}$, $\phi$ is of the form

$\phi(y)=Py+b\;.$ (40)

###### Proof.

First observe that for $P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ and
$b\in\mathbb{R}^{d}$, $y\mapsto Py+b$ is an isometry since we have, for any
$y$ and $y^{\prime}$ in $\mathbb{R}^{d^{\prime}}$

$\|Py+b-Py^{\prime}-b\|^{2}=\|P(y-y^{\prime})\|^{2}=(y-y^{\prime})^{T}P^{T}P(y-y^{\prime})=(y-y^{\prime})^{T}(y-y^{\prime})=\|y-y^{\prime}\|^{2}\;.$
(41)

The converse is a consequence of the Mazur–Ulam theorem Mazur and Ulam, (1932)
that states \- in the version of Baker, (1971) \- that an isometry from a real
normed space to a _strictly convex_ normed space, i.e. 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 $\\{x\in\mathbb{R}^{d}\,:\;\|x\|\leq 1\\}$ is a
strictly convex set, we get that for all $x\in\mathbb{R}^{d^{\prime}}$, $\phi$
is of the form $y\mapsto Py+b$ with $P$ being a matrix of size
$d\times{d^{\prime}}$, and $b\in\mathbb{R}^{d}$. Moreover we have for all
$y,y^{\prime}\in\mathbb{R}^{d^{\prime}}$

$\|\phi(y)-\phi(y^{\prime})\|^{2}=\|Py-
Py^{\prime}\|^{2}=\|P(y-y^{\prime})\|^{2}=(y-y^{\prime})^{T}P^{T}P(y-y^{\prime})\;.$
(42)

Since $\phi$ is an isometry, it follows that
$\|y-y^{\prime}\|^{2}=(y-y^{\prime})^{T}P^{T}P(y-y^{\prime})$ and so
$P^{T}P=\operatorname{Id}_{d^{\prime}}$, which concludes the proof. ∎

###  A.3 Centering of measures

Finally, we show the following result.

######  Lemma A0.

Let $\mu\in\mathcal{W}_{2}(\mathbb{R}^{d})$ and
$\nu\in\mathcal{W}_{2}(\mathbb{R}^{d^{\prime}})$ with $d$ not necessarily
greater than ${d^{\prime}}$. Let $\bar{\mu}$ and $\bar{\nu}$ denote the
centered measures associated to $\mu$ and $\nu$ and let $\mathfrak{P}$ be any
subset of matrices of size $d\times{d^{\prime}}$. Then,

$\inf_{\pi\in\Pi(\mu,\nu)}\inf_{P\in\mathfrak{P},\
b\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-Py-b\|^{2}\mathrm{d}\pi(x,y)=\inf_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\inf_{P\in\mathfrak{P}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-Py\|^{2}\mathrm{d}\pi(x,y)\;.$
(43)

###### Proof.

Denoting $m_{0}=\mathbb{E}_{X\sim\mu}[X]$, $m_{1}=\mathbb{E}_{Y\sim\nu}[Y]$,
$\tilde{x}=x-m_{0}$, and $\tilde{y}=y-m_{1}$, we have for any
$\pi\in\Pi(\mu,\nu)$,

$\displaystyle\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-Py-b\|^{2}\mathrm{d}\pi(x,y)$
$\displaystyle=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|\tilde{x}-P\tilde{y}-b+m_{0}-Pm_{1}\|^{2}\mathrm{d}\pi(x,y)$
(44)
$\displaystyle=\|m_{0}-b-Pm_{1}\|^{2}+\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|\tilde{x}-P\tilde{y}\|^{2}\mathrm{d}\pi(x,y)\;,$
(45)

since
$\int\langle\tilde{x}-P\tilde{y},m_{0}-b-Pm_{1}\rangle\mathrm{d}\pi(x,y)=0$.
Thus it follows,

$\displaystyle\inf_{\pi\in\Pi(\mu,\nu)}\inf_{P\in\mathfrak{P},\
b\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}$
$\displaystyle\|x-Py-b\|^{2}\mathrm{d}\pi(x,y)$ (47)
$\displaystyle=\inf_{P\in\mathfrak{P}}\left(\inf_{b\in\mathbb{R}^{d}}\|m_{0}-Pm_{1}-b\|^{2}+\inf_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-Py\|^{2}\mathrm{d}\pi(x,y)\right)\;.$
(48)

Observe now that for any $P\in\mathfrak{P}$, $\|m_{0}-Pm_{1}-b\|^{2}=0$ if
$b=m_{0}-Pm_{1}$, which concludes the proof. ∎

###  A.4 Some properties of symmetric matrices

Here we state two technical results on symmetric matrices that will be useful
during the proofs of the results on Gaussian distributions.

######  Lemma A0.

Let $A\in\mathbb{S}^{d}$. We denote $\lambda_{1}$ and $\lambda_{d}$ its
largest and smallest eigenvalues. For all $x\in\mathbb{R}^{d}$ such that
$\|x\|=1$, we have

  * (i)

$x$ is an eigenvector of $A$ associated to $\lambda_{1}$ if and only if
$x^{T}Ax=\lambda_{1}$.

  * (ii)

$x$ is an eigenvector of $A$ associated to $\lambda_{d}$ if and only if
$x^{T}Ax=\lambda_{d}$.

###### Proof.

Let $x\in\mathbb{R}^{d}$ such $\|x\|=1$. Since $A$ is symmetric, there exists
$O\in\mathbb{O}(\mathbb{R}^{d})$ and
$\Lambda=\mathrm{diag}((\lambda_{k})_{1\leq k\leq d})$ such that
$x^{T}Ax=x^{T}OAO^{T}x$. Denoting $z$ the vector $O^{T}x$, we get thus

$x^{T}Ax=z^{T}\Lambda z=\sum_{k=1}^{d}\lambda_{k}z_{k}^{2}\;.$ (49)

Hence it follows that

$\lambda_{d}\|z\|^{2}\leq x^{T}Ax\leq\lambda_{1}\|z\|^{2}\;,$ (50)

with equality if and only if $z$ is an eigenvector associated with
$\lambda_{1}$ or $\lambda_{d}$. ∎

######  Lemma A0.

Suppose that $d\geq{d^{\prime}}$. Let $\Sigma$ be a positive semi-definite
(PSD) matrix of size $d+{d^{\prime}}$ of the form

$\Sigma=\begin{pmatrix}\Sigma_{0}&K\\\ K^{T}&\Sigma_{1}\end{pmatrix}\;,$ (51)

with $\Sigma_{0}\in\mathbb{S}_{++}^{d}$,
$\Sigma_{1}\in\mathbb{S}_{+}^{d^{\prime}}$ and $K$ being a rectangular matrix
of size $d\times{d^{\prime}}$. Let $S=\Sigma_{1}-K^{T}\Sigma_{0}^{-1}K$ be the
Schur complement of $\Sigma$. Then there exists $r\leq{d^{\prime}}$ and
$B_{r}\in\mathbb{V}_{r}(\mathbb{R}^{d})$ such that

$K=\Sigma_{0}^{\frac{1}{2}}B_{r}\Lambda_{r}U_{r}^{T}\;,$ (52)

where $U_{r}\in\mathbb{V}_{r}(\mathbb{R}^{d^{\prime}})$ and $\Lambda_{r}$ is a
diagonal positive matrix of size $r$ such that

$\Sigma_{1}-S=U_{r}\Lambda^{2}_{r}U_{r}^{T}\;.$ (53)

###### Proof.

For a given Schur complement $S=\Sigma_{1}-K^{T}\Sigma_{0}^{-1}K$, we have
$K^{T}\Sigma_{0}^{-1}K=\Sigma_{1}-S$. Since
$\Sigma_{0}\in\mathbb{S}_{++}^{d}$, we can deduce that
$K^{T}\Sigma_{0}^{-1}K\in\mathbb{S}_{+}^{d^{\prime}}$ and so that
$\Sigma_{1}-S\in\mathbb{S}_{+}^{d^{\prime}}$. We note $r$ the rank of
$K^{T}\Sigma_{0}^{-1}K$. One can observe that

$r\leq{d^{\prime}}\leq d\;,$ (54)

where the left-hand side inequality follows from the fact that
$\text{rk}(AB)\leq\min\\{\text{rk}(A),\text{rk}(B)\\}$. Then, $\Sigma_{1}-S$
can be diagonalized

$\Sigma_{1}-S=K^{T}\Sigma_{0}^{-1}K=U\Lambda^{2}U^{T}=U_{r}\Lambda^{2}_{r}U_{r}^{T}\;,$
(55)

with
$\Lambda^{2}=\text{diag}(\lambda_{1}^{2},...,\lambda_{r}^{2})^{[{d^{\prime}}]}$,
$\Lambda_{r}^{2}=\text{diag}(\lambda_{1}^{2},...,\lambda_{r}^{2})$, and
$U_{r}\in\mathbb{V}_{r}(\mathbb{R}^{d^{\prime}})$ such that
$U=\begin{pmatrix}U_{r}&U_{{d^{\prime}}-r}\end{pmatrix}$. From (55), we can
deduce that

$(\Sigma_{0}^{-\frac{1}{2}}KU_{r}\Lambda_{r}^{-1})^{T}\Sigma_{0}^{-\frac{1}{2}}KU_{r}\Lambda_{r}^{-1}=\operatorname{Id}_{r}\;,$
(56)

where $\Lambda_{r}$ is the unique PSD square-root of $\Lambda^{2}_{r}$. Let us
set $B_{r}=\Sigma_{0}^{-\frac{1}{2}}KU_{r}\Lambda_{r}^{-1}$ such that
$B_{r}\in\mathbb{V}_{r}(\mathbb{R}^{d})$. It follows that

$KU_{r}=\Sigma_{0}^{\frac{1}{2}}B_{r}\Lambda_{r}\;.$ (57)

Moreover, since $U_{d-r}^{T}K^{T}\Sigma_{0}^{-1}KU_{d-r}=0$ and $\Sigma_{0}\in
S_{d}^{++}(\mathbb{R})$, it follows that $KU_{{d^{\prime}}-r}=0$ and so

$K=KUU^{T}=KU_{r}U_{r}^{T}=\Sigma_{0}^{\frac{1}{2}}B_{r}\Lambda_{r}U_{r}^{T}\;,$
(58)

which concludes the proof. ∎

##  Appendix B Proofs of the theoretical results

###  B.1 Proof of Theorem 1

###### Proof of Theorem 1.

Takatsu, (2010) has shown that the space of Gaussian distributions
$\mathcal{N}(\mathbb{R}^{d})$ is a complete metric space when endowed with
$W_{2}$. Moreover, $\mathcal{N}(\mathbb{R}^{d})$ is separable since it is a
subspace of $\mathcal{W}_{2}(\mathbb{R}^{d})$ which is itself a separable
metric space when endowed with $W_{2}$ Bolley, (2008). Thus,
$\mathcal{N}(\mathbb{R}^{d})$ is Polish and we can directly apply the Gromov-
Wasserstein theory developped in Sturm, (2012). Let
$(\mathcal{N}(\mathbb{R}^{d}),W_{2},\tilde{\mu})$ and
$(\mathcal{N}(\mathbb{R}^{d^{\prime}}),W_{2},\tilde{\nu})$ be two metric
measure spaces respectively in $\mathbb{M}_{4}$. Let us define

$D(\tilde{\mu},\tilde{\nu})=\inf_{\pi\in\Pi(\tilde{\mu},\tilde{\nu})}\int_{\mathcal{N}(\mathbb{R}^{d})\times\mathcal{N}(\mathbb{R}^{d^{\prime}})}\int_{\mathcal{N}(\mathbb{R}^{d})\times\mathcal{N}(\mathbb{R}^{d^{\prime}})}|W_{2}^{2}(\gamma,\gamma^{\prime})-W_{2}^{2}(\zeta,\zeta^{\prime})|^{2}\mathrm{d}\pi(\gamma,\zeta)\mathrm{d}\pi(\gamma^{\prime},\zeta^{\prime})\;.$
(59)

Applying Sturm, (2012, Corollary 9.3), we get that $D$ defines a metric over
the space of metric measure spaces of the form
$(\mathcal{N}(\mathbb{R}^{d}),W_{2},\tilde{\mu})$ quotiented by the strong
isomorphisms, and thus we get directly that $D$ is symmetric, non-negative,
satisfies the triangle inequality and $D(\tilde{\mu},\tilde{\nu})=0$ if and
only if there exists a bijection
$\phi:\mathrm{supp}(\tilde{\mu})\rightarrow\mathrm{supp}(\tilde{\nu})$ such
that $\tilde{\nu}=\phi_{\\#}\tilde{\mu}$, where for any $\gamma$ and
$\gamma^{\prime}$ in $\mathrm{supp}(\tilde{\mu})$,
$W_{2}(\phi(\gamma),\phi(\gamma^{\prime}))=W_{2}(\gamma,\gamma^{\prime})$. Now
observe that if $\mu=\sum_{k}a_{k}\mu_{k}$ and $\nu=\sum_{l}b_{l}\nu_{l}$ are
respectively in $GMM_{K}(\mathbb{R}^{d})$ and
$GMM_{L}(\mathbb{R}^{d^{\prime}})$ and
$\tilde{\mu}=\sum_{k}a_{k}\delta_{\mu_{k}}$ and
$\tilde{\nu}=\sum_{l}b_{l}\delta_{\nu_{l}}$ are respectively in
$\mathcal{P}(\mathcal{N}(\mathbb{R}^{d}))$ and
$\mathcal{P}(\mathcal{N}(\mathbb{R}^{d^{\prime}}))$, we have

$\int_{\mathcal{N}(\mathbb{R}^{d})\times\mathcal{N}(\mathbb{R}^{d})}W_{2}^{4}(\gamma,\gamma^{\prime})\mathrm{d}\tilde{\mu}(\gamma)\mathrm{d}\tilde{\mu}(\gamma^{\prime})=\sum_{k,i}a_{k}a_{i}W_{2}^{4}(\mu_{k},\mu_{i})<+\infty\;,$
(60)

and

$\int_{\mathcal{N}(\mathbb{R}^{d^{\prime}})\times\mathcal{N}(\mathbb{R}^{d^{\prime}})}W_{2}^{4}(\zeta,\zeta^{\prime})\mathrm{d}\tilde{\nu}(\zeta)\mathrm{d}\tilde{\nu}(\zeta^{\prime})=\sum_{l,j}b_{l}b_{j}W_{2}^{4}(\nu_{l},\nu_{j})<+\infty\;,$
(61)

so $(\mathcal{N}(\mathbb{R}^{d}),W_{2},\tilde{\mu})$ and
$(\mathcal{N}(\mathbb{R}^{d^{\prime}}),W_{2},\tilde{\nu})$ are both in
$\mathbb{M}_{4}$. Furthermore, we have
$MGW_{2}(\mu,\nu)=D(\tilde{\mu},\tilde{\nu})$. Hence $MGW_{2}$ inherits the
metric properties of $D$, which concludes the proof. ∎

###  B.2 Proof of Section 3.1

###### Proof of Section 3.1..

First recall that the push-foward measure $T_{\\#}\mu$ with $\mu$ on
$\mathbb{R}^{d^{\prime}}$ and
$T:\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d}$ is defined as the measure
on $\mathbb{R}^{d}$ such that for every Borel set $\mathsf{A}$ of
$\mathbb{R}^{d}$, $T_{\\#}\mu(\mathsf{A})=\mu(T^{-1}(\mathsf{A}))$.
Equivalently, for any measurable map $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$,
we have

$\int_{\mathbb{R}^{d}}h(x)\mathrm{d}(T_{\\#}\mu)(x)=\int_{\mathbb{R}^{d^{\prime}}}(h\circ
T)(y)\mathrm{d}\mu(y)\;.$ (62)

Now observe that for any finite GMM $\mu$ on $\mathbb{R}^{d^{\prime}}$ of the
form $\mu=\sum_{k}^{K}a_{k}\mu_{k}$, we have

$\displaystyle\textstyle{\int_{\mathbb{R}^{d^{\prime}}}(h\circ
T)(y)\mathrm{d}\mu(y)}$
$\displaystyle=\textstyle{\int_{\mathbb{R}^{d^{\prime}}}(h\circ
T)(y)\mathrm{d}\left(\sum_{k}^{K}a_{k}\mu_{k}(y)\right)}$ (63)
$\displaystyle=\textstyle{\sum_{k}^{K}a_{k}\int_{\mathbb{R}^{d^{\prime}}}(h\circ
T)(y)\mathrm{d}\mu_{k}(y)}$ (64)
$\displaystyle=\textstyle{\sum_{k}^{K}a_{k}\int_{\mathbb{R}^{d}}h(x)\mathrm{d}(T_{\\#}\mu_{k})(x)}$
(65)
$\displaystyle=\textstyle{\int_{\mathbb{R}^{d}}h(x)\mathrm{d}\left(\sum_{k}^{K}a_{k}(T_{\\#}\mu_{k})(x)\right)}\;,$
(66)

and so $T_{\\#}\mu$ is of the form $\sum_{k}^{K}a_{k}(T_{\\#}\mu_{k})$ with
$T_{\\#}\mu_{k}$ Gaussian since $T$ is necessarily affine as a consequence of
Section A.2. Thus, $T_{\\#}\mu$ is in $GMM_{\infty}(\mathbb{R}^{d})$. This
proves that $\phi_{T}$ takes its values only in $GMM_{\infty}(\mathbb{R}^{d})$
and that $\phi_{T}(\sum_{k=1}^{K}a_{k}\mu_{k})$ is of the form
$\sum_{k=1}a_{k}\nu_{k}$. Now observe that, for every $k$ and $i$ smaller than
$K$,

$W^{2}_{2}(\phi_{T}(\mu_{k}),\phi_{T}(\mu_{i}))=\inf_{\pi\in\Pi(T_{\\#}\mu_{k},T_{\\#}\mu_{i})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x-y\|^{2}\mathrm{d}\pi(x,y)\;.$
(67)

Using two times succesively Section A.1 using the fact that $T$ is an isometry
an so for any $y\in\mathbb{R}^{d^{\prime}}$, $\|T(y)\|=\|y\|$, it follows

$\inf_{\pi\in\Pi(T_{\\#}\mu_{k},T_{\\#}\mu_{i})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x-x^{\prime}\|^{2}\mathrm{d}\pi(x,x^{\prime})=\inf_{\pi\in\Pi(\mu_{k},\mu_{i})}\int_{\mathbb{R}^{d^{\prime}}\times\mathbb{R}^{d^{\prime}}}\|y-y^{\prime}\|^{2}\mathrm{d}\pi(y,y^{\prime})=W_{2}(\mu_{k},\mu_{i})\;.$
(68)

Thus, $MGW_{2}(\mu,T_{\\#}\mu)=0$ as a direct consequence of Theorem 1, which
concludes the proof. ∎

###  B.3 Proof of Section 4.1

###### Proof of Section 4.1.

Since we suppose $d\geq{d^{\prime}}$, we have

$EW_{2}^{2}(\mu,\nu)=\inf_{\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})}W_{2}^{2}(\mu,\phi_{\\#}\nu)\;.$
(69)

Let $\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})$ for the Euclidean
norm. Using Section A.2, we get that there exists
$P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$ and $b\in\mathbb{R}^{d}$ such
that for all $y\in\mathbb{R}^{d^{\prime}}$, $\phi(y)=Py+b$. Moreover, we have,
using Section A.1,

$\displaystyle EW_{2}^{2}(\mu,\nu)$
$\displaystyle=\inf_{\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})}\inf_{\pi\in\Pi(\mu,\phi_{\\#}\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x-y\|^{2}\mathrm{d}\pi(x,y)$
(70)
$\displaystyle=\inf_{\phi\in\mathrm{Isom}_{d^{\prime}}(\mathbb{R}^{d})}\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{d^{\prime}}\times\mathbb{R}^{d}}\|x-\phi(y)\|^{2}\mathrm{d}\pi(x,y)$
(71)
$\displaystyle=\inf_{\pi\in\Pi(\mu,\nu)}\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d}),\
b\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d^{\prime}}\times\mathbb{R}^{d}}\|x-Py-b\|^{2}\mathrm{d}\pi(x,y)\;,$
(72)

which concludes the proof. ∎

###  B.4 Proof of Section 4.1

Before turning to the proof of Section 4.1, we will prove the following lemma.

######  Lemma B0.

Let $K$ be a matrix of size $d\times{d^{\prime}}$ with Singular Value
Decomposition (SVD) $K=U_{K}\Sigma_{K}V_{K}^{T}$ and let $\mathfrak{P}$ be any
compact set of matrices of size $d\times{d^{\prime}}$. Then,

$\sup_{P\in\mathfrak{P}}\mathrm{tr}(P^{T}K)=\max_{P\in\mathfrak{P}}\mathrm{tr}(\Sigma^{T}_{P}\Sigma_{K})\;,$
(73)

where $\Sigma_{P}=\mathrm{diag}^{[d,{d^{\prime}}]}(\bm{\sigma}(P))$ with
$\bm{\sigma}(P)\in\mathbb{R}_{+}^{d^{\prime}}$ denoting the vector of singular
values of P. Furthermore it is achieved at $P$ of the form,

$P=U_{K}\Sigma_{P}V_{K}^{T}\;.$ (74)

###### Proof of Section B.4.

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 $\mathfrak{P}$ is compact and the mapping
$P\mapsto\mathrm{tr}(P^{T}K)$ is continuous. For a given $P\in\mathfrak{P}$,
let $U_{P}\Sigma_{P}V_{P}^{T}$ be the SVD of $P$. The problem can be rewritten
as

$\max_{P\in\mathfrak{P}}\mathrm{tr}(V_{P}\Sigma^{T}_{P}U_{P}^{T}U_{K}\Sigma_{K}V_{K}^{T})\;.$
(75)

Now, let us denote $U=U_{P}^{T}U_{K}$ and $V=V_{P}^{T}V_{K}$. Observe that $U$
is in $\mathbb{O}(\mathbb{R}^{d})$ and $V$ is in
$\mathbb{O}(\mathbb{R}^{d^{\prime}})$. Using the cyclical permutation of the
trace operator, the problem becomes

$\max_{P\in\mathfrak{P}}\mathrm{tr}(\Sigma^{T}_{P}U\Sigma_{K}V^{T})\;.$ (76)

Now, for a given fixed $\Sigma_{P}$, we determine which $U$ and $V$ maximize
$\mathrm{tr}(\Sigma^{T}_{P}U\Sigma_{K}V^{T})$. This problem reads as

$\max_{U\in\mathbb{O}(\mathbb{R}^{d}),\
V\in\mathbb{O}(\mathbb{R}^{d^{\prime}})}\mathrm{tr}(\Sigma^{T}_{P}U\Sigma_{K}V^{T})\;.$
(77)

The Lagrangian of this problem reads as

$\mathcal{L}(U,V,C_{0},C_{1})=-\mathrm{tr}(\Sigma^{T}_{P}U\Sigma_{K}V^{T})+\mathrm{tr}(C_{0}(U^{T}U-\operatorname{Id}_{d}))+\mathrm{tr}(C_{1}(V^{T}V-\operatorname{Id}_{d^{\prime}}))\;,$
(78)

where $C_{0}\in\mathbb{S}^{d}$ and $C_{1}\in\mathbb{S}^{d^{\prime}}$ are the
Lagrange multipliers respectively associated with the constraints
$U\in\mathbb{O}(\mathbb{R}^{d})$ and
$V\in\mathbb{O}(\mathbb{R}^{d^{\prime}})$. The first order condition gives

$\left\\{\begin{array}[]{ll}\Sigma_{P}V\Sigma_{K}^{T}=2UC_{0}\\\
\Sigma_{P}^{T}U\Sigma_{K}=2VC_{1}\;,\end{array}\right.$ (79)

or equivalently

$\left\\{\begin{array}[]{ll}U^{T}\Sigma_{P}V\Sigma_{K}^{T}=2C_{0}\\\
\Sigma_{P}^{T}U\Sigma_{K}V^{T}=2VC_{1}V^{T}\;.\end{array}\right.$ (80)

Since $C_{0}$ and $C_{1}$ are symmetric matrices (because they are associated
with symmetric constraints), we get that both left-hand terms are symmetric.
This gives the following conditions

$\left\\{\begin{array}[]{ll}U^{T}\Sigma_{P}V\Sigma_{K}^{T}=\Sigma_{K}V^{T}\Sigma_{P}^{T}U\\\
\Sigma_{P}^{T}U\Sigma_{K}V^{T}=V\Sigma_{K}^{T}U^{T}\Sigma_{P}\;.\end{array}\right.$
(81)

Now, observe that when multiplying the first condition at right by
$U^{T}\Sigma_{P}$ and multiplying the second condition at left by
$\Sigma_{K}V^{T}$, we get by combining the two conditions

$\left\\{\begin{array}[]{ll}U\Sigma_{K}V^{T}\Sigma_{P}^{T}\Sigma_{P}=\Sigma_{P}\Sigma_{P}^{T}U\Sigma_{K}V^{T}\\\
U^{T}\Sigma_{P}V\Sigma_{K}^{T}\Sigma_{K}=\Sigma_{K}\Sigma_{K}^{T}U^{T}\Sigma_{P}V^{T}\;,\end{array}\right.$
(82)

or equivalently,

$\left\\{\begin{array}[]{ll}U\Sigma_{K}V^{T}D_{P}=D_{P}^{[d]}U\Sigma_{K}V^{T}\\\
U^{T}\Sigma_{P}VD_{K}=D_{K}^{[d]}U^{T}\Sigma_{P}V^{T}\;,\end{array}\right.$
(83)

where $D_{P}=\mathrm{diag}(\bm{\sigma}(P))$ and
$D_{K}=\mathrm{diag}(\bm{\sigma}(K))$. Multiplying the first condition at left
by $V\Sigma_{K}^{T}U^{T}$ and the second condition at right by
$V\Sigma_{P}^{T}U$, this yields to

$\left\\{\begin{array}[]{ll}VD_{K}V^{T}D_{P}=V\Sigma
U^{T}D_{P}^{[d]}U\Sigma_{K}V^{T}\\\
D_{K}^{[d]}U^{T}D_{P}^{[d]}U=U^{T}\Sigma_{P}VD_{K}V\Sigma_{P}^{T}U\;.\end{array}\right.$
(84)

It follows that $VD_{K}V^{T}D_{P}$ and $D_{K}^{[d]}U^{T}D_{P}^{[d]}U$ are
symmetric matrices and so $VD_{K}V^{T}$ commutes with $D_{P}$ and
$U^{T}D_{P}^{[d]}U$ commutes with $D_{K}^{[d]}$. Thus we can deduce that $U$
and $V$ are permutation matrices. Since the singular values are ordered in
non-increasing order, we deduce that the problem is maximized when
$U=\operatorname{Id}_{d}$ and $V=\operatorname{Id}_{{d^{\prime}}}$. This
implies that $U_{P}=U_{K}$ and $V_{P}=V_{K}$, which concludes the proof. ∎

Note that Section B.4 is especially useful when the constraint of belonging to
the set $\mathfrak{P}$ can be expressed as a constraint on the singular
values. Observe that this is the case of
$\mathbb{V}_{{d^{\prime}}}(\mathbb{R}^{d})$ since for all
$P\in\mathbb{V}_{{d^{\prime}}}(\mathbb{R}^{d})$, we have
$P^{T}P=\operatorname{Id}_{{d^{\prime}}}$ and so an equivalent condition of
belonging in $\mathbb{V}_{{d^{\prime}}}(\mathbb{R}^{d})$ is that
$\bm{\sigma}(P)=\mathbbm{1}_{d^{\prime}}$. Now we are ready to turn to the
proof of Section 4.1.

###### Proof of Section 4.1.

Using Lemmas A.3, Problem (18) can be rewritten

$\displaystyle EW_{2}^{2}(\mu,\nu)$
$\displaystyle=\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\inf_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\|x-Py\|^{2}\mathrm{d}\pi(x,y)$
(85)
$\displaystyle=\inf_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\inf_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\left(\|x\|^{2}+\|Py\|^{2}-2\langle
x,Py\rangle\right)\mathrm{d}\pi(x,y)\;.$ (86)

Since for all $P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$, $\|Py\|$ doesn’t
depend on $P$, we get that the problem is equivalent to

$\sup_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\langle
x,Py\rangle\mathrm{d}\pi(x,y)\;.$ (87)

Now observe that for all $\pi\in\Pi(\bar{\mu},\bar{\nu})$,

$\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\langle
x,Py\rangle\mathrm{d}\pi(x,y)=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\mathrm{tr}(xy^{T}P^{T})\mathrm{d}\pi(x,y)=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}\mathrm{tr}(P^{T}xy^{T})\mathrm{d}\pi(x,y)\;,$
(88)

where we used the cyclical permutation property of the trace operator. Finally
using the linearity of the trace, we get that the problem is equivalent to

$\sup_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\mathrm{tr}\left(P^{T}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}\mathrm{d}\pi(x,y)\right)\;,$
(89)

or equivalently,

$\sup_{P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})}\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\left\langle
P,\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}\mathrm{d}\pi(x,y)\right\rangle\;.$
(90)

Now, using Section B.4 and using the fact that if
$P\in\mathbb{V}_{d^{\prime}}(\mathbb{R}^{d})$,
$\bm{\sigma}(P)=\mathbbm{1}_{d^{\prime}}$, we get that the problem reduces to

$\sup_{\pi\in\Pi(\bar{\mu},\bar{\nu})}\left\|\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}\mathrm{d}\pi(x,y)\right\|_{*}\;,$
(91)

and this is achieved for
$P^{*}=U_{\pi}\operatorname{Id}_{{d^{\prime}}}^{[d,{d^{\prime}}]}V_{\pi}^{T}$,
where $U_{\pi}\in\mathbb{O}(\mathbb{R}^{d})$ and
$V_{\pi}\in\mathbb{O}(\mathbb{R}^{d^{\prime}})$ are respectively the left and
right orthogonal matrices of the SVD of
$\int_{\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}}xy^{T}\mathrm{d}\pi(x,y)$,
which concludes the proof. ∎

###  B.5 Proof of Section 4.1

###### Proof of Section 4.1.

First, using Section A.2, we get that there exists
$P_{1}\in\mathbb{V}_{d}(\mathbb{R}^{r})$ and $b_{1}\in\mathbb{R}^{r}$ such