A Wasserstein-type distance in the space of Gaussian Mixture Models

TL;DR

This paper introduces a new Wasserstein-type distance specifically designed for Gaussian mixture models, enabling efficient computation and practical applications in high-dimensional image processing tasks.

Contribution

It proposes a novel Wasserstein-type distance restricted to Gaussian mixtures, with a simple discrete formulation suitable for high-dimensional problems.

Findings

Simple discrete formulation for the distance

Properties of the Wasserstein-type distance analyzed

Practical applications demonstrated in image processing

Abstract

In this paper we introduce a Wasserstein-type distance on the set of Gaussian mixture models. This distance is defined by restricting the set of possible coupling measures in the optimal transport problem to Gaussian mixture models. We derive a very simple discrete formulation for this distance, which makes it suitable for high dimensional problems. We also study the corresponding multi-marginal and barycenter formulations. We show some properties of this Wasserstein-type distance, and we illustrate its practical use with some examples in image processing.

Equations

$$\underset{Y_0\sim {\mu }_0;Y_1\sim {\mu }_1}{\inf }\mathbb{E}\left(c(Y_0,Y_1)\right),$$

$${\int }_{{\mathbb{R}}^d}\Vert x{\Vert }^{p}d\mu (x)<+\infty ,$$

$$\forall x,y\in {\mathbb{R}}^d,{\mathrm{P}}_t(x,y)=(1-t)x+ty\in {\mathbb{R}}^d.$$

$$W_p^p({\mu }_0,{\mu }_1):=\underset{Y_0\sim {\mu }_0;Y_1\sim {\mu }_1}{\inf }\mathbb{E}\left(\Vert Y_0-Y_1{\Vert }^p\right)=\underset{\gamma \in \Pi ({\mu }_0,{\mu }_1)}{\inf }{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\Vert y_0-y_1{\Vert }^{p}d\gamma (y_0,y_1).$$

$$\gamma =(\mathrm{Id},T)\#{\mu }_0,$$

$$\forall t\in [0,1],{\mu }_t:={\mathrm{P}}_t\#\gamma$$

$${\mu }_t\in {\mathrm{argmin}}_{\rho }(1-t)W_2({\mu }_0,\rho {)}^2+t{W}_2({\mu }_1,\rho {)}^2.$$

$$B(x)=\sum _{i=0}^{J-1}{\lambda }_i{x}_i={\mathrm{argmin}}_{y\in {\mathbb{R}}^d}\sum _{i=0}^{J-1}{\lambda }_i\Vert x_i-y{\Vert }^2$$

$$\underset{\nu \in {\mathcal{P}}_2({\mathbb{R}}^d)}{\inf }\sum _{j=0}^{J-1}{\lambda }_j{W}_2^2({\mu }_j,\nu ).$$

$$mm{W}_2^2({\mu }_0,\dots ,{\mu }_{J-1})$$

$$W_2^2({\mu }_0,{\mu }_1)=\Vert m_0-m_1{\Vert }^2+\mathrm{tr}\left({\Sigma }_0+{\Sigma }_1-2{\left({\Sigma }_0^{\frac{1}{2}}{\Sigma }_1{\Sigma }_0^{\frac{1}{2}}\right)}^{\frac{1}{2}}\right),$$

$$\forall x\in {\mathbb{R}}^d,T(x)=m_1+{\Sigma }_0^{-\frac{1}{2}}{\left({\Sigma }_0^{\frac{1}{2}}{\Sigma }_1{\Sigma }_0^{\frac{1}{2}}\right)}^{\frac{1}{2}}{\Sigma }_0^{-\frac{1}{2}}(x-m_0)=m_1+{\Sigma }_0^{-1}({\Sigma }_0{\Sigma }_1{)}^{\frac{1}{2}}(x-m_0),$$

$${\Sigma }_t=((1-t){\mathrm{I}}_d+tC){\Sigma }_0((1-t){\mathrm{I}}_d+tC),$$

$${\gamma }^{\ast }(x_0,\dots ,x_{J-1})=g_{m_0,{\Sigma }_0}(x_0){\delta }_{(x_1,\dots ,x_{J-1})=(S_1{S}_0^{-1}{x}_0,\dots ,S_{J-1}{S}_0^{-1}{x}_0)}$$

$$\sum _{j=0}^{J-1}{\lambda }_j{\left({\Sigma }_{\ast }^{1/2}{\Sigma }_j{\Sigma }_{\ast }^{1/2}\right)}^{1/2}={\Sigma }_{\ast }.$$

$${\Gamma }_N=\{\pi \in {\mathbb{R}}_{+}^N;\pi 1_N=\sum _{k=1}^N{\pi }_k=1\}.$$

$$\mu =\sum _{k=1}^K{\pi }_k{\mu }_k\text{ where }{\mu }_k=\mathcal{N}(m_k,{\Sigma }_k)\text{ and }\pi \in {\Gamma }_K.$$

$$GM{M}_d(\infty )={\cup }_{K\ge 0}GM{M}_d(K).$$

$$\left\{\sum _{k=1}^N{\pi }_k{\delta }_{y_k};N\in \mathbb{N},(y_k{)}_k\in ({\mathbb{R}}^d{)}^N,({\pi }_k{)}_k\in {\Gamma }_N\right\}$$

$${\cup }_k\left(m_k+\mathrm{Span}({\Sigma }_k)\right)=\mathrm{graph}(T).$$

$$f_t(x)=\frac{1}{1-t}\left(g\left(\frac{x+t}{1-t}\right)1_{x<-t}+g\left(\frac{x-t}{1-t}\right)1_{x>t}\right),$$

$$M{W}_2^2({\mu }_0,{\mu }_1):=\underset{\gamma \in \Pi ({\mu }_0,{\mu }_1)\cap GM{M}_{2d}(\infty )}{\inf }{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\Vert y_0-y_1{\Vert }^{2}d\gamma (y_0,y_1).$$

$$M{W}_2({\mu }_0,{\mu }_1)\ge W_2({\mu }_0,{\mu }_1).$$

$$\Pi ({\pi }_0,{\pi }_1)$$

$$=\{w\in {\mathcal{M}}_{K_0,K_1}({\mathbb{R}}^{+});\forall k,\sum _j{w}_{kj}={\pi }_0^k\text{ and }\forall j,\sum _k{w}_{kj}={\pi }_1^j\}.$$

$$M{W}_2^2({\mu }_0,{\mu }_1)=\underset{w\in \Pi ({\pi }_0,{\pi }_1)}{\min }\sum _{k,l}{w}_{kl}{W}_2^2({\mu }_0^k,{\mu }_1^l).$$

$${\gamma }^{\ast }(x,y)=\sum _{k,l}{w}_{k,l}^{\ast }{g}_{m_0^k,{\Sigma }_0^k}(x){\delta }_{y=T_{k,l}(x)}$$

$$\underset{w\in \Pi ({\pi }_0,{\pi }_1)}{\inf }\sum _{k,l}{w}_{kl}{W}_2^2({\mu }_0^k,{\mu }_1^l).$$

$${\gamma }_{kl}={\mathrm{argmin}}_{\gamma \in \Pi ({\mu }_0^k,{\mu }_1^l)}{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\Vert y_0-y_1{\Vert }^{2}d\gamma (y_0,y_1)$$

$${\gamma }^{\ast }=\sum _{k,l}{w}_{kl}^{\ast }{\gamma }_{kl}.$$

$$\sum _{k,l}{w}_{kl}^{\ast }{W}_2^2({\mu }_0^k,{\mu }_1^l)$$

Full text

{AMS}

65K10, 65K05, 90C05, 62-07, 68Q25, 68U10, 68U05, 68R10

1 Introduction

Nowadays, Gaussian Mixture Models (GMM) have become ubiquitous in statistics and machine learning. These models are especially useful in applied fields to represent probability distributions of real datasets. Indeed, as linear combinations of Gaussian distributions, they are perfect to model complex multimodal densities and can approximate any continuous density when the numbers of components is chosen large enough. Their parameters are also easy to infer with algorithms such as the Expectation-Maximization (EM) algorithm [13]. For instance, in image processing, a large body of works use GMM to represent patch distributions in images111Patches are small image pieces, they can be seen as vectors in a high dimensional space., and use these distributions for various applications, such as image restoration [36, 28, 35, 32, 19, 12] or texture synthesis [16].

The optimal transport theory provides mathematical tools to compare or interpolate between probability distributions. For two probability distributions $\mu_{0}$ and $\mu_{1}$ on $\mathbb{R}^{d}$ and a positive cost function $c$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$, the goal is to solve the optimization problem

[TABLE]

where the notation $Y\sim\mu$ means that $Y$ is a random variable with probability distribution $\mu$. When $c(x,y)=\|x-y\|^{p}$ for $p\geq 1$, Equation (1) (to a power $1/p$) defines a distance between probability distributions that have a moment of order $p$, called the Wasserstein distance $W_{p}$.

While this subject has gathered a lot of theoretical work (see [30, 31, 27] for three reference monographies on the topic), its success in applied fields was slowed down for many years by the computational complexity of numerical algorithms which were not always compatible with large amount of data. In recent years, the development of efficient numerical approaches has been a game changer, widening the use of optimal transport to various applications notably in image processing, computer graphics and machine learning [23]. However, computing Wasserstein distances or optimal transport plans remains intractable when the dimension of the problem is too high.

Optimal transport can be used to compute distances or geodesics between Gaussian mixture models, but optimal transport plans between GMM, seen as probability distributions on a higher dimensional space, are usually not Gaussian mixture models themselves, and the corresponding Wasserstein geodesics between GMM do not preserve the property of being a GMM. In order to keep the good properties of these models, we define in this paper a variant of the Wasserstein distance by restricting the set of possible coupling measures to Gaussian mixture models. The idea of restricting the set of possible coupling measures has already been explored for instance in [3], where the distance is defined on the set of the probability distributions of strong solutions to stochastic differential equations. The goal of the authors is to define a distance which keeps the good properties of $W_{2}$ while being numerically tractable.

In this paper, we show that restricting the set of possible coupling measures to Gaussian mixture models transforms the original infinitely dimensional optimization problem into a finite dimensional problem with a simple discrete formulation, depending only on the parameters of the different Gaussian distributions in the mixture. When the ground cost is simply $c(x,y)=\|x-y\|^{2}$, this yields a geodesic distance, that we call $MW_{2}$ (for Mixture Wasserstein), which is obviously larger than $W_{2}$, and is always upper bounded by $W_{2}$ plus a term depending only on the trace of the covariance matrices of the Gaussian components in the mixture. The complexity of the corresponding discrete optimization problem does not depend on the space dimension, but only on the number of components in the different mixtures, which makes it particularly suitable in practice for high dimensional problems. Observe that this equivalent discrete formulation has been proposed twice recently in the machine learning literature, but with a very different point of view, by two independent teams [8, 9] and [6, 7].

Our original contributions in this paper are the following:

We derive an explicit formula for the optimal transport between two GMM restricted to GMM couplings, and we show several properties of the resulting distance, in particular how it compares to the classical Wasserstein distance. 2. 2.

We study the multi-marginal and barycenter formulations of the problem, and show the link between these formulations. 3. 3.

We propose a generalized formulation to be used on distributions that are not GMM. 4. 4.

We provide two applications in image processing, respectively to color transfer and texture synthesis.

The paper is organized as follows. Section 2 is a reminder on Wasserstein distances and barycenters between probability measures on $\mathbb{R}^{d}$. We also recall the explicit formulation of $W_{2}$ between Gaussian distributions. In Section 3, we recall some properties of Gaussian mixture models, focusing on an identifiabiliy property that will be necessary for the rest of the paper. We also show that optimal transport plans for $W_{2}$ between GMM are generally not GMM themselves. Then, Section 4 introduces the $MW_{2}$ distance and derives the corresponding discrete formulation. Section 4.5 compares $MW_{2}$ with $W_{2}$. Section 5 focuses on the corresponding multi-marginal and barycenter formulations. In Section 6, we explain how to use $MW_{2}$ in practice on distributions that are not necessarily GMM. We conclude in Section 7 with two applications of the distance $MW_{2}$ to image processing. To help the reproducibility of the results we present in this paper, we have made our Python codes available on the Github website https://github.com/judelo/gmmot.

Notations

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

• •

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

• •

If $\mu$ is a positive measure on a space $\mathcal{X}$ and $T:\mathcal{X}\to\mathcal{Y}$ is an application, $T\#\mu$ stands for the push-forward measure of $\mu$ by $T$, i.e. the measure on $\mathcal{Y}$ such that $\forall A\subset\mathcal{Y}$, $(T\#\mu)(A)=\mu(T^{-1}(A))$.

• •

The notation $\mathrm{tr}(M)$ denotes the trace of the matrix $M$.

• •

The notation $\mathrm{Id}$ is the identity application.

• •

$\langle\xi,\xi^{\prime}\rangle$ denotes the Euclidean scalar product between $\xi$ and $\xi^{\prime}$ in $\mathbb{R}^{d}$

• •

$\mathcal{M}_{n,m}(\mathbb{R})$ is the set of real matrices with $n$ lines and $m$ columns, and we denote by $\mathcal{M}_{n_{0},n_{1},\dots,n_{J-1}}(\mathbb{R})$ the set of $J$ dimensional tensors of size $n_{k}$ in dimension $k$.

• •

$\mathbf{1}_{n}=(1,1,\dots,1)^{t}$ denotes a column vector of ones of length $n$.

• •

For a given vector $m$ in $\mathbb{R}^{d}$ and a $d\times d$ covariance matrix $\Sigma$, $g_{m,\Sigma}$ denotes the density of the Gaussian (multivariate normal) distribution $\mathcal{N}(m,\Sigma)$.

• •

When $a_{i}$ is a finite sequence of $K$ elements (real numbers, vectors or matrices), we denote its elements as $a_{i}^{0},\ldots,a_{i}^{K-1}$.

2 Reminders: Wasserstein distances and barycenters between probability measures on $\mathbb{R}^{d}$

Let $d\geq 1$ be an integer. We recall in this section the definition and some basic properties of the Wasserstein distances between probability measures on $\mathbb{R}^{d}$. We write $\mathcal{P}(\mathbb{R}^{d})$ the set probability measures on $\mathbb{R}^{d}$. For $p\geq 1$, the Wasserstein space $\mathcal{P}_{p}(\mathbb{R}^{d})$ is defined as the set of probability measures $\mu$ with a finite moment of order $p$, i.e. such that

[TABLE]

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

For $t\in[0,1]$, we define $\mathrm{P}_{t}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ by

[TABLE]

Observe that $\mathrm{P}_{0}$ and $\mathrm{P}_{1}$ are the projections from $\mathbb{R}^{d}\times\mathbb{R}^{d}$ onto $\mathbb{R}^{d}$ such that $\mathrm{P}_{0}(x,y)=x$ and $\mathrm{P}_{1}(x,y)=y$.

2.1 Wasserstein distances

Let $p\geq 1$, and let $\mu_{0},\mu_{1}$ be two probability measures in $\mathcal{P}_{p}(\mathbb{R}^{d})$. Define $\Pi(\mu_{0},\mu_{1})\subset\mathcal{P}_{p}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ as being the subset of probability distributions $\gamma$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginal distributions $\mu_{0}$ and $\mu_{1}$, i.e. such that $\mathrm{P}_{0}\#\gamma=\mu_{0}$ and $\mathrm{P}_{1}\#\gamma=\mu_{1}$. The $p$-Wasserstein distance $W_{p}$ between $\mu_{0}$ and $\mu_{1}$ is defined as

[TABLE]

This formulation is a special case of (1) when $c(x,y)=\|x-y\|^{p}$. It can be shown (see for instance [31]) that there is always a couple $(Y_{0},Y_{1})$ of random variables which attains the infimum (hence a minimum) in the previous energy. Such a couple is called an optimal coupling. The probability distribution $\gamma$ of this couple is called an optimal transport plan between $\mu_{0}$ and $\mu_{1}$. This plan distributes all the mass of the distribution $\mu_{0}$ onto the distribution $\mu_{1}$ with a minimal cost, and the quantity $W_{p}^{p}(\mu_{0},\mu_{1})$ is the corresponding total cost.

As suggested by its name ($p$-Wasserstein distance), $W_{p}$ defines a metric on $\mathcal{P}_{p}(\mathbb{R}^{d})$. It also metrizes the weak convergence222A sequence $(\mu_{k})_{k}$ converges weakly to $\mu$ in $\mathcal{P}_{p}(\mathbb{R}^{d})$ if it converges to $\mu$ in the sense of distributions and if $\int\|y\|^{p}d\mu_{k}(y)$ converges to $\int\|y\|^{p}d\mu(y)$. in $\mathcal{P}_{p}(\mathbb{R}^{d})$ (see [31], chapter 6). It follows that $W_{p}$ is continuous on $\mathcal{P}_{p}(\mathbb{R}^{d})$ for the topology of weak convergence.

From now on, we will mainly focus on the case $p=2$, since $W_{2}$ has an explicit formulation if $\mu_{0}$ and $\mu_{1}$ are Gaussian measures.

2.2 Transport map, transport plan and displacement interpolation

Assume that $p=2$. When $\mu_{0}$ and $\mu_{1}$ are two probability distributions on $\mathbb{R}^{d}$ and assuming that $\mu_{0}$ is absolutely continuous, then it can be shown that the optimal transport plan $\gamma$ for the problem (2) is unique and has the form

[TABLE]

where $T:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ is an application called optimal transport map and satisfying $T\#\mu_{0}=\mu_{1}$ (see [31]).

If $\gamma$ is an optimal transport plan for $W_{2}$ between two probability distributions $\mu_{0}$ and $\mu_{1}$, the path $(\mu_{t})_{t\in[0,1]}$ given by

[TABLE]

defines a constant speed geodesic in $\mathcal{P}_{2}(\mathbb{R}^{d})$ (see for instance [27] Ch.5, Section 5.4). The path $(\mu_{t})_{t\in[0,1]}$ is called the displacement interpolation between $\mu_{0}$ and $\mu_{1}$ and it satisifes

[TABLE]

This interpolation, often called Wasserstein barycenter in the literature, can be easily extended to more than two probability distributions, as recalled in the next paragraphs.

2.3 Multi-marginal formulation and barycenters

For $J\geq 2$, for a set of weights $\lambda=(\lambda_{0},\dots,\lambda_{J-1})\in(\mathbb{R}_{+})^{J}$ such that $\lambda\mathbf{1}_{J}=\lambda_{0}+\ldots+\lambda_{J-1}=1$ and for $x=(x_{0},\dots,x_{J-1})\in(\mathbb{R}^{d})^{J}$, we write

[TABLE]

the barycenter of the $x_{i}$ with weights $\lambda_{i}$.

For $J$ probability distributions $\mu_{0},\mu_{1}\dots,\mu_{J-1}$ on $\mathbb{R}^{d}$, we say that $\nu^{*}$ is the barycenter of the $\mu_{j}$ with weights $\lambda_{j}$ if $\nu^{*}$ is solution of

[TABLE]

Existence and unicity of barycenters for $W_{2}$ has been studied in depth by Agueh and Carlier in [1]. They show in particular that if one of the $\mu_{j}$ has a density, this barycenter is unique. They also show that the solutions of the barycenter problem are related to the solutions of the multi-marginal transport problem (studied by Gangbo and Świéch in [17])

[TABLE]

where $\Pi(\mu_{0},\mu_{1},\dots,\mu_{J-1})$ is the set of probability measures on $(\mathbb{R}^{d})^{J}$ having $\mu_{0},\mu_{1},\dots,\mu_{J-1}$ as marginals. More precisely, they show that if (7) has a solution $\gamma^{*}$, then $\nu^{*}=B\#\gamma^{*}$ is a solution of (6), and the infimum of (7) and (6) are equal.

2.4 Optimal transport between Gaussian distributions

Computing optimal transport plans between probability distributions is usually difficult. In some specific cases, an explicit solution is known. For instance, in the one dimensional ($d=1$) case, when the cost $c$ is a convex function of the Euclidean distance on the line, the optimal plan consists in a monotone rearrangement of the distribution $\mu_{0}$ into the distribution $\mu_{1}$ (the mass is transported monotonically from left to right, see for instance Ch.2, Section 2.2 of [30] for all the details). Another case where the solution is known for a quadratic cost is the Gaussian case in any dimension $d\geq 1$.

2.4.1 Distance $W_{2}$ between Gaussian distributions

If $\mu_{i}=\mathcal{N}(m_{i},\Sigma_{i})$, $i\in\{0,1\}$ are two Gaussian distributions on $\mathbb{R}^{d}$, the $2$-Wasserstein distance $W_{2}$ between $\mu_{0}$ and $\mu_{1}$ has a closed-form expression, which can be written

[TABLE]

where, for every symmetric semi-definite positive matrix $M$, the matrix $M^{\frac{1}{2}}$ is its unique semi-definite positive square root.

If $\Sigma_{0}$ is non-singular, then the optimal map $T$ between $\mu_{0}$ and $\mu_{1}$ turns out to be affine and is given by

[TABLE]

and the optimal plan $\gamma$ is then a Gaussian distribution on $\mathbb{R}^{d}\times\mathbb{R}^{d}=\mathbb{R}^{2d}$ that is degenerate since it is supported by the affine line $y=T(x)$. These results have been known since [14].

Moreover, if $\Sigma_{0}$ and $\Sigma_{1}$ are non-degenerate, the geodesic path $(\mu_{t})$, $t\in(0,1)$, between $\mu_{0}$ and $\mu_{1}$ is given by $\mu_{t}=\mathcal{N}(m_{t},\Sigma_{t})$ with $m_{t}=(1-t)m_{0}+tm_{1}$ and

[TABLE]

with $\mathrm{I}_{d}$ the $d\times d$ identity matrix and $C=\Sigma_{1}^{\frac{1}{2}}\left(\Sigma_{1}^{\frac{1}{2}}\Sigma_{0}\Sigma_{1}^{\frac{1}{2}}\right)^{{-\frac{1}{2}}}\Sigma_{1}^{\frac{1}{2}}$.

This property still holds if the covariance matrices are not invertible, by replacing the inverse by the Moore-Penrose pseudo-inverse matrix, see Proposition 6.1 in [33]. The optimal map $T$ is not generalized in this case since the optimal plan is usually not supported by the graph of a function.

2.4.2 $W_{2}$-Barycenters in the Gaussian case

For $J\geq 2$, let $\lambda=(\lambda_{0},\dots,\lambda_{J-1})\in(\mathbb{R}_{+})^{J}$ be a set of positive weights summing to $1$ and let $\mu_{0},\mu_{1}\dots,\mu_{J-1}$ be $J$ Gaussian probability distributions on $\mathbb{R}^{d}$. For $j=$ $0\dots J-1$, we denote by $m_{j}$ and $\Sigma_{j}$ the expectation and the covariance matrix of $\mu_{j}$. Theorem 2.2 in [26] tells us that if the covariances $\Sigma_{j}$ are all positive definite, then the solution of the multi-marginal problem (7) for the Gaussian distributions $\mu_{0},\mu_{1}\dots,\mu_{J-1}$ can be written

[TABLE]

where $S_{j}=\Sigma_{j}^{1/2}\left(\Sigma_{j}^{1/2}\Sigma_{*}\Sigma_{j}^{1/2}\right)^{-1/2}\Sigma_{j}^{1/2}$ with $\Sigma_{*}$ a solution of the fixed-point problem

[TABLE]

The barycenter $\nu^{*}$ of all the $\mu_{j}$ with weights $\lambda_{j}$ is the distribution $\mathcal{N}(m_{*},\Sigma_{*})$, with $m_{*}=\sum_{j=0}^{J-1}\lambda_{j}m_{j}$. Equation (11) provides a natural iterative algorithm (see [2]) to compute the fixed point $\Sigma_{*}$ from the set of covariances $\Sigma_{j}$, $j\in\{0,\dots,J-1\}$.

3 Some properties of Gaussian Mixtures Models

The goal of this paper is to investigate how the optimisation problem (2) is transformed when the probability distributions $\mu_{0}$, $\mu_{1}$ are finite Gaussian mixture models and the transport plan $\gamma$ is forced to be a Gaussian mixture model. This will be the aim of Section 4. Before, we first need to recall a few basic properties on these mixture models, and especially a density property and an identifiability property.

In the following, for $N\geq 1$ integer, we define the simplex

[TABLE]

Definition 1**.**

Let $K\geq 1$ be an integer. A (finite) Gaussian mixture model of size $K$ on $\mathbb{R}^{d}$ is a probability distribution $\mu$ on $\mathbb{R}^{d}$ that can be written

[TABLE]

We write $GMM_{d}(K)$ the subset of $\mathcal{P}(\mathbb{R}^{d})$ made of probability measures on $\mathbb{R}^{d}$ which can be written as Gaussian mixtures with less than $K$ components (such mixtures are obviously also in $\mathcal{P}_{p}(\mathbb{R}^{d})$ for any $p\geq 1$). For $K<K^{\prime}$, $GMM_{d}(K)\subset GMM_{d}(K^{\prime}).$ The set of all finite Gaussian mixture distributions is written

[TABLE]

3.1 Two properties of GMM

The following lemma states that any measure in $\mathcal{P}_{p}(\mathbb{R}^{d})$ can be approximated with any precision for the distance $W_{p}$ by a finite convex combination of Dirac masses. This classical result will be useful in the rest of the paper.

Lemma 3.1**.**

The set

[TABLE]

*is dense in $\mathcal{P}_{p}(\mathbb{R}^{d})$ for the metric $W_{p}$, for any $p\geq 1$. *

For the sake of completeness, we provide a proof in Appendix, adapted from the proof of Theorem 6.18 in [31]. Since Dirac masses can be seen as degenerate Gaussian distributions, a direct consequence of Lemma 3.1 is the following proposition.

Proposition 1**.**

$GMM_{d}(\infty)$* is dense in $\mathcal{P}_{p}(\mathbb{R}^{d})$ for the metric $W_{p}$. *

Another important property will be necessary, related to the identifiability of Gaussian mixture models. It is clear such models are not stricto sensu identifiable, since reordering the indexes of a mixture changes its parametrization without changing the underlying probability distribution, or also because a component with mass $1$ can be divided in two identical components with masses $\frac{1}{2}$, for example. However, if we write mixtures in a “compact” way (forbidding two components of the same mixture to be identical), identifiability holds, up to a reordering of the indexes. This property is reminded below.

Proposition 2**.**

*The set of finite Gaussian mixtures is identifiable, in the sense that two mixtures $\mu_{0}=\sum_{k=1}^{K_{0}}\pi_{0}^{k}\mu_{0}^{k}$ and $\mu_{1}=\sum_{k=1}^{K_{1}}\pi_{1}^{k}\mu_{1}^{k}$, written such that all $\{\mu_{0}^{k}\}_{k}$ (resp. all $\{\mu_{1}^{j}\}_{j}$) are pairwise distinct, are equal if and only if $K_{0}=K_{1}$ and we can reorder the indexes such that for all $k$, $\pi_{0}^{k}=\pi_{1}^{k}$, $m_{0}^{k}=m_{1}^{k}$ and $\Sigma_{0}^{k}=\Sigma_{1}^{k}$.

This result is also classical and the proof is provided in Appendix.

3.2 Optimal transport and Wasserstein barycenters between Gaussian Mixture Models

We are now in a position to investigate optimal transport between Gaussian mixture models (GMM). A first important remark is that given two Gaussian mixtures $\mu_{0}$ and $\mu_{1}$ on $\mathbb{R}^{d}$, optimal transport plans $\gamma$ between $\mu_{0}$ and $\mu_{1}$ are usually not GMM.

Proposition 3**.**

*Let $\mu_{0}\in GMM_{d}(K_{0})$ and $\mu_{1}\in GMM_{d}(K_{1})$ be two Gaussian mixtures such that $\mu_{1}$ cannot be written $T\#\mu_{0}$ with $T$ affine. Assume also that $\mu_{0}$ is absolutely continuous with respect to the Lebesgue measure. Let $\gamma\in\Pi(\mu_{0},\mu_{1})$ be an optimal transport plan between $\mu_{0}$ and $\mu_{1}$. Then $\gamma$ does not belongs to $GMM_{2d}(\infty)$. *

Proof 3.2**.**

Since $\mu_{0}$ is absolutely continuous with respect to the Lebesgue measure, we know that the optimal transport plan is unique and is of the form $\gamma=(\mathrm{Id},T)\#\mu_{0}$ for a measurable map $T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ that satisfies $T\#\mu_{0}=\mu_{1}$. Thus, if $\gamma$ belongs to $GMM_{2d}(\infty)$, all of its components must be degenerate Gaussian distributions $\mathcal{N}(m_{k},\Sigma_{k})$ such that

[TABLE]

*It follows that $T$ must be affine on $\mathbb{R}^{d}$, which contradicts the hypotheses of the proposition. *

When $\mu_{0}$ is not absolutely continuous with respect to the Lebesgue measure (which means that one of its components is degenerate), we cannot write $\gamma$ under the form (3), but we conjecture that the previous result usually still holds. A notable exception is the case where all Gaussian components of $\mu_{0}$ and $\mu_{1}$ are Dirac masses on $\mathbb{R}^{d}$, in which case $\gamma$ is also a GMM composed of Dirac masses on $\mathbb{R}^{2d}$.

We conjecture that since optimal plans $\gamma$ between two GMM are usually not GMM, the barycenters $(\mathrm{P}_{t})\#\gamma$ between $\mu_{0}$ and $\mu_{1}$ are also usually not GMM either (with the exception of $t=0,1$). Take the one dimensional example of $\mu_{0}=\mathcal{N}(0,1)$ and $\mu_{1}=\frac{1}{2}(\delta_{-1}+\delta_{1})$. Clearly, an optimal transport map between $\mu_{0}$ and $\mu_{1}$ is defined as $T(x)=\mathrm{sign}(x)$. For $t\in(0,1)$, if we denote by $\mu_{t}$ the barycenter between $\mu_{0}$ with weight $1-t$ and $\mu_{1}$ with weight $t$, then it is easy to show that $\mu_{t}$ has a density

[TABLE]

where $g$ is the density of $\mathcal{N}(0,1)$. The density $f_{t}$ is equal to [math] on the interval $(-t,t)$ and therefore cannot be the density of a GMM.

4 $MW_{2}$: a distance between Gaussian Mixture Models

In this section, we define a Wasserstein-type distance between Gaussian mixtures ensuring that barycenters between Gaussian mixtures remain Gaussian mixtures. To this aim, we restrict the set of admissible transport plans to Gaussian mixtures and show that the problem is well defined. Thanks to the identifiability results proved in the previous section, we will show that the corresponding optimization problem boils down to a very simple discrete formulation.

4.1 Definition of $MW_{2}$

Definition 2**.**

Let $\mu_{0}$ and $\mu_{1}$ be two Gaussian mixtures. We define

[TABLE]

First, observe that the problem is well defined since $\Pi(\mu_{0},\mu_{1})\cap GMM_{2d}(\infty)$ contains at least the product measure $\mu_{0}\otimes\mu_{1}$. Notice also that from the definition we directly have that

[TABLE]

4.2 An equivalent discrete formulation

Now, we can show that this optimisation problem has a very simple discrete formulation. For $\pi_{0}\in\Gamma_{K_{0}}$ and $\pi_{1}\in\Gamma_{K_{1}}$, we denote by $\Pi(\pi_{0},\pi_{1})$ the subset of the simplex $\Gamma_{K_{0}\times K_{1}}$ with marginals $\pi_{0}$ and $\pi_{1}$, i.e.

[TABLE]

Proposition 4**.**

Let $\mu_{0}=\sum_{k=1}^{K_{0}}\pi_{0}^{k}\mu_{0}^{k}$ and $\mu_{1}=\sum_{k=1}^{K_{1}}\pi_{1}^{k}\mu_{1}^{k}$ be two Gaussian mixtures, then

[TABLE]

Moreover, if $w^{*}$ is a minimizer of (16), and if $T_{k,l}$ is the $W_{2}$-optimal map between $\mu_{0}^{k}$ and $\mu_{1}^{l}$, then $\gamma^{*}$ defined as

[TABLE]

*is a minimizer of (13). *

Proof 4.1**.**

First, let $w^{*}$ be a solution of the discrete linear program

[TABLE]

For each pair $(k,l)$, let

[TABLE]

and

[TABLE]

Clearly, $\gamma^{*}\in\Pi(\mu_{0},\mu_{1})\cap GMM_{2d}(K_{0}K_{1})$. It follows that

[TABLE]

because $GMM_{2d}(K_{0}K_{1})\subset GMM_{2d}(\infty)$.

Now, let $\gamma$ be any element of $\Pi(\mu_{0},\mu_{1})\cap GMM_{2d}(\infty)$. Since $\gamma$ belongs to $GMM_{2d}(\infty)$, there exists an integer $K$ such that $\gamma=\sum_{j=1}^{K}w_{j}\gamma_{j}$. Since $\mathrm{P}_{0}{\#}\gamma=\mu_{0}$, it follows that

[TABLE]

Thanks to the identifiability property shown in the previous section, we know that these two Gaussian mixtures must have the same components, so for each $j$ in $\{1,\dots K\}$, there is $1\leq k\leq K_{0}$ such that $\mathrm{P}_{0}{\#}\gamma_{j}=\mu_{0}^{k}$. In the same way, there is $1\leq l\leq K_{1}$ such that $\mathrm{P}_{1}{\#}\gamma_{j}=\mu_{1}^{l}$. It follows that $\gamma_{j}$ belongs to $\Pi(\mu_{0}^{k},\mu_{1}^{l})$. We conclude that the mixture $\gamma$ can be written as a mixture of Gaussian components $\gamma_{kl}\in\Pi(\mu_{0}^{k},\mu_{1}^{l})$, i.e $\gamma=\sum_{k=1}^{K_{0}}\sum_{l=1}^{K_{1}}w_{kl}\gamma_{kl}$. Since $\mathrm{P}_{0}{\#}\gamma=\mu_{0}$ and $\mathrm{P}_{1}{\#}\gamma=\mu_{1}$, we know that $w\in\Pi(\pi_{0},\pi_{1})$. As a consequence,

[TABLE]

*This inequality holds for any $\gamma$ in $\Pi(\mu_{0},\mu_{1})\cap GMM_{2d}(\infty)$, which concludes the proof. *

It happens that the discrete form (16), which can be seen as an aggregation of simple Wasserstein distances between Gaussians, has been recently proposed as an ingenious alternative to $W_{2}$ in the machine learning literature, both in [8, 9] and [6, 7]. Observe that the point of view followed here in our paper is quite different from these works, since $MW_{2}$ is defined in a completely continuous setting as an optimal transport between GMMs with a restriction on couplings, following the same kind of approach as in [3]. The fact that this restriction leads to an explicit discrete formula, the same as the one proposed independently in [8, 9] and [6, 7], is quite striking. Observe also that thanks to the “identifiability property” of GMMs, this continuous formulation (13) is obviously non ambiguous, in the sense that the value of the minimium is the same whatever the parametrization of the Gaussian mixtures $\mu_{0}$ and $\mu_{1}$. This was not obvious from the discrete versions. We will see in the following sections how this continuous formulation can be extended to multi-marginal and barycenter formulations, and how it can be generalized or used in the case of more general distributions.

Notice that we do not use in the definition and in the proof the fact that the ground cost is quadratic. Definition 2 can thus be generalized to other cost functions $c:\mathbb{R}^{2d}\mapsto\mathbb{R}$. The reason why we focus on the quadratic cost is that optimal transport plans between Gaussian measures for $W_{2}$ can be computed explicitely. It follows from the equivalence between the continuous and discrete forms of $MW_{2}$ that the solution of (13) is very easy to compute in practice. Another consequence of this equivalence is that there exists at least one optimal plan $\gamma^{*}$ for (13) containing less than $K_{0}+K_{1}-1$ Gaussian components.

Corollary 1**.**

*Let $\mu_{0}=\sum_{k=1}^{K_{0}}\pi_{0}^{k}\mu_{0}^{k}$ and $\mu_{1}=\sum_{k=1}^{K_{1}}\pi_{1}^{k}\mu_{1}^{k}$ be two Gaussian mixtures on $\mathbb{R}^{d}$, then the infimum in (13) is attained for a given $\gamma^{*}\in\Pi(\mu_{0},\mu_{1})\cap GMM_{2d}(K_{0}+K_{1}-1)$. *

Proof 4.2**.**

*This follows directly from the proof that there exists at least one optimal $w^{*}$ for (13) containing less than $K_{0}+K_{1}-1$ Gaussian components (see [23]). *

4.3 An example in one dimension

In order to illustrate the behavior of the optimal maps for $MW_{2}$, we focus here on a very simple example in one dimension, where $\mu_{0}$ and $\mu_{1}$ are the following mixtures of two Gaussian components

[TABLE]

[TABLE]

Figure 1 shows the optimal transport plans between $\mu_{0}$ (in blue) and $\mu_{1}$ (in red), both for the Wasserstein distance $W_{2}$ and for $MW_{2}$. As we can observe, the optimal transport plan for $MW_{2}$ (a probability measure on $\mathbb{R}\times\mathbb{R}$) is a mixture of three degenerate Gaussians measures supported by 1D lines.

4.4 Metric properties of $MW_{2}$ and displacement interpolation

4.4.1 Metric properties of $MW_{2}$

Proposition 5**.**

$MW_{2}$* defines a metric on $GMM_{d}(\infty)$ and the space $GMM_{d}(\infty)$ equipped with the distance $MW_{2}$ is a geodesic space. *

This proposition can be proved very easily by making use of the discrete formulation (16) of the distance (see for instance [7]). For the sake of completeness, we provide in the following a proof of the proposition using only the continuous formulation of $MW_{2}$.

Proof 4.3**.**

First, observe that $MW_{2}$ is obviously symmetric and positive. It is also clear that for any Gaussian mixture $\mu$, $MW_{2}(\mu,\mu)=0$. Conversely, assume that $MW_{2}(\mu_{0},\mu_{1})=0$, it implies that $W_{2}(\mu_{0},\mu_{1})=0$ and thus $\mu_{0}=\mu_{1}$ since $W_{2}$ is a distance.

It remains to show that $MW_{2}$ satisfies the triangle inequality. This is a classical consequence of the gluing lemma, but we must be careful to check that the constructed measure remains a Gaussian mixture. Let $\mu_{0},\mu_{1},\mu_{2}$ be three Gaussian mixtures on $\mathbb{R}^{d}$. Let $\gamma_{01}$ and $\gamma_{12}$ be optimal plans respectively for $(\mu_{0},\mu_{1})$ and $(\mu_{1},\mu_{2})$ for the problem $MW_{2}$ (which means that $\gamma_{01}$ and $\gamma_{12}$ are both GMM on $\mathbb{R}^{2d}$). The classical gluing lemma consists in disintegrating $\gamma_{01}$ and $\gamma_{12}$ into

[TABLE]

and to define

[TABLE]

which boils down to assume independence conditionnally to the value of $y_{1}$. Since $\gamma_{01}$ and $\gamma_{12}$ are Gaussian mixtures on $\mathbb{R}^{2d}$, the conditional distributions $d\gamma_{01}(y_{0}|y_{1})$ and $d\gamma_{12}(y_{2}|y_{1})$ are also Gaussian mixtures for all $y_{1}$ in the support of $\mu_{1}$ (recalling that $\mu_{1}$ is the marginal on $y_{1}$ of both $\gamma_{01}$ and $\gamma_{12}$). If we define a distribution $\gamma_{02}$ by integrating $\gamma_{012}$ over the variable $y_{1}$, i.e.

[TABLE]

then $\gamma_{02}$ is obviously also a Gaussian mixture on $\mathbb{R}^{2d}$ with marginals $\mu_{0}$ and $\mu_{2}$. The rest of the proof is classical. Indeed, we can write

[TABLE]

Writing $\|y_{0}-y_{2}\|^{2}=\|y_{0}-y_{1}\|^{2}+\|y_{1}-y_{2}\|^{2}+2\langle y_{0}-y_{1},y_{1}-y_{2}\rangle$ (with $\langle\;,\;\rangle$ the Euclidean scalar product on $\mathbb{R}^{d}$), and using the Cauchy-Schwarz inequality, it follows that

[TABLE]

The triangle inequality follows by taking for $\gamma_{01}$ (resp. $\gamma_{12}$) the optimal plan for $MW_{2}$ between $\mu_{0}$ and $\mu_{1}$ (resp. $\mu_{1}$ and $\mu_{2}$).

Now, let us show that $GMM_{d}(\infty)$ equipped with the distance $MW_{2}$ is a geodesic space. For a path $\rho=(\rho_{t})_{t\in[0,1]}$ in $GMM_{d}(\infty)$ (meaning that each $\rho_{t}$ is a GMM on $\mathbb{R}^{d}$), we can define its length for $MW_{2}$ by

[TABLE]

Let $\mu_{0}=\sum_{k}\pi_{0}^{k}\mu_{0}^{k}$ and $\mu_{1}=\sum_{l}\pi_{1}^{l}\mu_{1}^{l}$ be two GMM. Since $MW_{2}$ satifies the triangle inequality, we always have that $\mathrm{Len}_{MW_{2}}(\rho)\geq MW_{2}(\mu_{0},\mu_{1})$ for all paths $\rho$ such that $\rho_{0}=\mu_{0}$ and $\rho_{1}=\mu_{1}$. To prove that $(GMM_{d}(\infty),MW_{2})$ is a geodesic space we just have to exhibit a path $\rho$ connecting $\mu_{0}$ to $\mu_{1}$ and such that its length is equal to $MW_{2}(\mu_{0},\mu_{1})$.

We write $\gamma^{\ast}$ the optimal transport plan between $\mu_{0}$ and $\mu_{1}$. For $t\in(0,1)$ we can define

[TABLE]

Let $t<s\in[0,1]$ and define $\gamma_{t,s}^{\ast}=(\mathrm{P}_{t},\mathrm{P}_{s})\#\gamma^{\ast}$. Then $\gamma_{t,s}^{\ast}\in\Pi(\mu_{t},\mu_{s})\cap GMM_{2d}(\infty)$ and therefore

[TABLE]

Thus we have that $MW_{2}(\mu_{t},\mu_{s})\leq(s-t)MW_{2}(\mu_{0},\mu_{1})$ Now, by the triangle inequality,

[TABLE]

*Therefore all inequalities are equalities, and $MW_{2}(\mu_{t},\mu_{s})=(s-t)MW_{2}(\mu_{0},\mu_{1})$ for all $0\leq t\leq s\leq 1$. This implies that the $MW_{2}$ length of the path $(\mu_{t})_{t}$ is equal to $MW_{2}(\mu_{0},\mu_{1})$. It allows us to conclude that $(GMM_{d}(\infty),MW_{2})$ is a geodesic space, and we have also given the explicit expression of the geodesic. *

The following Corollary is a direct consequence of the previous results.

Corollary 2**.**

The barycenters between $\mu_{0}=\sum_{k}\pi_{0}^{k}\mu_{0}^{k}$ and $\mu_{1}=\sum_{l}\pi_{1}^{l}\mu_{1}^{l}$ all belong to $GMM_{d}(\infty)$ and can be written explicitely as

[TABLE]

where $w^{\ast}$ is an optimal solution of (16), and $\mu_{t}^{k,l}$ is the displacement interpolation between $\mu_{0}^{k}$ and $\mu_{1}^{l}$. When $\Sigma_{0}^{k}$ is non-singular, it is given by

[TABLE]

*with $T_{k,l}$ the affine transport map between $\mu_{0}^{k}$ and $\mu_{1}^{l}$ given by Equation (9). These barycenters have less than $K_{0}+K_{1}-1$ components. *

4.4.2 1D and 2D barycenter examples

One dimensional case

Figure 2 shows barycenters $\mu_{t}$ for $t=0.2,0.4,0.6,0.8$ between the $\mu_{0}$ and $\mu_{1}$ defined in Section 4.3, for both the metric $W_{2}$ and $MW_{2}$. Observe that the barycenters computed for $MW_{2}$ are a bit more regular (we know that they are mixtures of at most 3 Gaussian components) than those obtained for $W_{2}$.

Two dimensional case

Figure 3 shows barycenters $\mu_{t}$ between the following two dimensional mixtures

[TABLE]

[TABLE]

where $I_{2}$ is the $2\times 2$ identity matrix. Notice that the $MW_{2}$ geodesic looks much more regular, each barycenter is a mixture of less than three Gaussians.

4.5 Comparison between $MW_{2}$ and $W_{2}$

Proposition 6**.**

Let $\mu_{0}\in GMM_{d}(K_{0})$ and $\mu_{1}\in GMM_{d}(K_{1})$ be two Gaussian mixtures, written as in (12). Then,

[TABLE]

The left-hand side inequality is attained when for instance

• •

$\mu_{0}$* and $\mu_{1}$ are both composed of only one Gaussian component,*

• •

$\mu_{0}$* and $\mu_{1}$ are finite linear combinations of Dirac masses,*

• •

$\mu_{1}$* is obtained from $\mu_{0}$ by an affine transformation.*

As we already noticed it, the first inequality is obvious and follows from the definition of $MW_{2}$. It might not be completely intuitive that $MW_{2}$ can indeed be strictly larger than $W_{2}$ because of the density property of $GMM_{d}(\infty)$ in $\mathcal{P}_{2}(\mathbb{R}^{d})$. This follows from the fact that our optimization problem has constraints $\gamma\in\Pi(\mu_{0},\mu_{1})$. Even if any measure $\gamma$ in $\Pi(\mu_{0},\mu_{1})$ can be approximated by a sequence of Gaussian mixtures, this sequence of Gaussian mixtures will generally not belong to $\Pi(\mu_{0},\mu_{1})$, hence explaining the difference between $MW_{2}$ and $W_{2}$.

In order to show that $MW_{2}$ is always smaller than the sum of $W_{2}$ plus a term depending on the trace of the covariance matrices of the two Gaussian mixtures, we start with a lemma which makes more explicit the distance $MW_{2}$ between a Gaussian mixture and a mixture of Dirac distributions.

Lemma 4.4**.**

Let $\mu_{0}=\sum_{k=1}^{K_{0}}\pi_{0}^{k}\mu_{0}^{k}$ with $\mu_{0}^{k}=\mathcal{N}(m_{0}^{k},\Sigma_{0}^{k})$ and $\mu_{1}=\sum_{k=1}^{K_{1}}\pi_{1}^{k}\delta_{m_{1}^{k}}$. Let $\tilde{\mu}_{0}=\sum_{k=1}^{K_{0}}\pi_{0}^{k}\delta_{m_{0}^{k}}$ ($\tilde{\mu}_{0}$ only retains the means of $\mu_{0}$). Then,

[TABLE]

Proof 4.5**.**

[TABLE]

In other words, the squared distance $MW_{2}^{2}$ between $\mu_{0}$ and $\mu_{1}$ is the sum of the squared Wasserstein distance between $\tilde{\mu}_{0}$ and $\mu_{1}$ and a linear combination of the traces of the covariance matrices of the components of $\mu_{0}$. We are now in a position to show the other inequality between $MW_{2}$ and $W_{2}$.

Proof 4.6** (Proof of Proposition 6).**

Let $(\mu_{0}^{n})_{n}$ and $(\mu_{1}^{n})_{n}$ be two sequences of mixtures of Dirac masses respectively converging to $\mu_{0}$ and $\mu_{1}$ in $\mathcal{P}_{2}(\mathbb{R}^{d})$. Since $MW_{2}$ is a distance,

[TABLE]

We study in the following the limits of these three terms when $n\rightarrow+\infty$.

First, observe that $MW_{2}(\mu_{0}^{n},\mu_{1}^{n})=W_{2}(\mu_{0}^{n},\mu_{1}^{n})\longrightarrow_{n\rightarrow\infty}W_{2}(\mu_{0},\mu_{1})$ since $W_{2}$ is continuous on $\mathcal{P}_{2}(\mathbb{R}^{d})$.

Second, using Lemma 4.4, for $i=0,1$,

[TABLE]

Define the measure $d\gamma(x,y)=\sum_{k=1}^{K_{i}}\pi_{i}^{k}\delta_{m_{i}^{k}}(y)g_{m_{i}^{k},\Sigma_{i}^{k}}(x)dx$, with $g_{m_{i}^{k},\Sigma_{i}^{k}}$ the probability density function of the Gaussian distribution $\mathcal{N}(m_{i}^{k},\Sigma_{i}^{k})$. The probability measure $\gamma$ belongs to $\Pi(\mu_{i},\tilde{\mu}_{i})$, so

[TABLE]

We conclude that

[TABLE]

This ends the proof of the proposition.

Observe that if $\mu$ is a Gaussian distribution $\mathcal{N}(m,\Sigma)$ and $\mu^{n}$ a distribution supported by a finite number of points which converges to $\mu$ in $\mathcal{P}_{2}(\mathbb{R}^{d})$, then

[TABLE]

and

[TABLE]

Let us also remark that if $\mu_{0}$ and $\mu_{1}$ are Gaussian mixtures such that $\max_{k,i}\mathrm{trace}(\Sigma_{i}^{k})\leq\varepsilon$, then

[TABLE]

4.6 Generalization to other mixture models

A natural question is to know if the methodology we have developped here, and that restricts the set of possible coupling measures to Gaussian mixtures, can be extended to other families of mixtures. Indeed, in the image processing litterature, as well as in many other fields, mixture models beyond Gaussian ones are widely used, such as Generalized Gaussian Mixture Models [11] or mixtures of T-distributions [29], for instance. Now, to extend our methodology to other mixtures, we need two main properties: (a) the identifiability property (that will ensure that there is a canonical way to write a distribution as a mixture); and (b) a marginal consistency property (we need all the marginal of an element of the family to remain in the same family). These two properties permit in particular to generalize the proof of Proposition 4. In order to make the discrete formulation convenient for numerical computations, we also need that the $W_{2}$ distance between any two elements of the family must be easy to compute.

Starting from this last requirement, we can consider a family of elliptical distributions, where the elements are of the form

[TABLE]

where $m\in\mathbb{R}^{d}$, $\Sigma$ is a positive definite symmetric matrix and $h$ is a given function from $[0,+\infty)$ to $[0,+\infty)$. Gaussian distributions are an example, with $h(t)=\exp(-t/2)$. Generalized Gaussian distributions are obtained with $h(t)=\exp(-t^{\beta})$, with $\beta$ not necessarily equal to $1$. T-distributions are also in this family, with $h(t)=(1+t/\nu)^{-(\nu+d)/2}$, etc. Thanks to their elliptical contoured property, the $W_{2}$ distance between two elements in such a family (i.e. $h$ fixed) can be explicitely computed (see Gelbrich [18]), and yields a formula that is the same as the one in the Gaussian case (Equation (8)). In such a family, the identifiability property can be checked, using the asymptotic behavior in all directions of $\mathbb{R}^{d}$. Now, if we want the marginal consistency property to be also satisfied (which is necessary if we want the coupling restriction problem to be well-defined), the choice of $h$ is very limited. Indeed, Kano in [21], proved that the only elliptical distributions with the marginal consistency property are the ones which are a scale mixture of normal distributions with a mixing variable that is unrelated to the dimension $d$. So, generalized Gaussian distributions don’t satisfy this marginal consistency property, but T-distributions do.

5 Multi-marginal formulation and barycenters

5.1 Multi-marginal formulation for $MW_{2}$

Let $\mu_{0},\mu_{1}\dots,\mu_{J-1}$ be $J$ Gaussian mixtures on $\mathbb{R}^{d}$, and let $\lambda_{0},\dots\lambda_{J-1}$ be $J$ positive weights summing to $1$. The multi-marginal version of our optimal transport problem restricted to Gaussian mixture models can be written

[TABLE]

where

[TABLE]

and where $\Pi(\mu_{0},\mu_{1},\dots,\mu_{J-1})$ is the set of probability measures on $(\mathbb{R}^{d})^{J}$ having $\mu_{0}$, $\mu_{1}$, $\dots$, $\mu_{J-1}$ as marginals.

Writing for every $j$, $\mu_{j}=\sum_{k=1}^{K_{j}}\pi_{j}^{k}\mu_{j}^{k}$, and using exactly the same arguments as in Proposition 4, we can easily show the following result.

Proposition 7**.**

The optimisation problem (18) can be rewritten under the discrete form

[TABLE]

where $\Pi(\pi_{0},\pi_{1},\dots,\pi_{J-1})$ is the subset of tensors $w$ in $\mathcal{M}_{K_{0},K_{1},\dots,K_{J-1}}(\mathbb{R}^{+})$ having $\pi_{0}$, $\pi_{1}$, $\dots$, $\pi_{J-1}$ as discrete marginals, i.e. such that

[TABLE]

Moreover, the solution $\gamma^{*}$ of (18) can be written

[TABLE]

*where $w^{*}$ is solution of (20) and $\gamma^{*}_{k_{0}k_{1}\dots k_{J-1}}$ is the optimal multi-marginal plan between the Gaussian measures $\mu_{0}^{k_{0}},\dots,\mu_{J-1}^{k_{J-1}}$ (see Section 2.4.2). *

From Section 2.4.2, we know how to construct the optimal multi-marginal plans $\gamma^{*}_{k_{0}k_{1}\dots k_{J-1}}$, which means that computing a solution for (18) boils down to solve the linear program (20) in order to find $w^{*}$.

5.2 Link with the $MW_{2}$-barycenters

We will now show the link between the previous multi-marginal problem and the barycenters for $MW_{2}$.

Proposition 8**.**

The barycenter problem

[TABLE]

*has a solution given by $\nu^{*}=B\#\gamma^{*}$, where $\gamma^{*}$ is an optimal plan for the multi-marginal problem (18). *

Proof 5.1**.**

For any $\gamma\in\Pi(\mu_{0},\dots,\mu_{J-1})\cap GMM_{Jd}(\infty)$, we define $\gamma_{j}=(P_{j},B)\#\gamma$, with $B$ the barycenter application defined in (5) and $P_{j}:(\mathbb{R}^{d})^{J}\mapsto\mathbb{R}^{d}$ such that $P(x_{0},\dots,x_{J-1})=x_{j}$. Observe that $\gamma_{j}$ belongs to $\Pi(\mu_{j},\nu)$ with $\nu=B\#\gamma$. The probability measure $\gamma_{j}$ also belongs to $GMM_{2d}(\infty)$ since $(P_{j},B)$ is a linear application. It follows that

[TABLE]

This inequality holds for any arbitrary $\gamma\in\Pi(\mu_{0},\dots,\mu_{J-1})\cap GMM_{Jd}(\infty)$, thus

[TABLE]

Conversely, for any $\nu$ in $GMM_{d}(\infty)$, we can write $\nu=\sum_{l=1}^{L}\pi_{\nu}^{l}\nu^{l}$, the $\nu^{l}$ being Gaussian probability measures. We also write $\mu_{j}=\sum_{k=1}^{K_{j}}\pi_{j}^{k}\mu_{j}^{k}$, and we call $w^{j}$ the optimal discrete plan for $MW_{2}$ between the mixtures $\mu_{j}$ and $\nu$ (see Equation (16)). Then,

[TABLE]

Now, if we define a $K_{0}\times\dots\times K_{J-1}\times L$ tensor $\alpha$ and a $K_{0}\times\dots\times K_{J-1}$ tensor $\overline{\alpha}$ by

[TABLE]

clearly $\alpha\in\Pi(\pi_{0},\dots,\pi_{J-1},\pi_{\nu})$ and $\overline{\alpha}\in\Pi(\pi_{0},\dots,\pi_{J-1})$. Moreover,

[TABLE]

*the last inequality being a consequence of Proposition 7. Since this holds for any arbitrary $\nu$ in $GMM_{d}(\infty)$, this ends the proof. *

The following corollary gives a more explicit formulation for the barycenters for $MW_{2}$, and shows that the number of Gaussian components in the mixture is much smaller than $\prod_{j=0}^{J-1}K_{j}$.

Corollary 3**.**

Let $\mu_{0},\dots,\mu_{J-1}$ be $J$ Gaussian mixtures such that all the involved covariance matrices are positive definite, then the solution of (23) can be written

[TABLE]

*where $\nu_{k_{0}\dots k_{J-1}}$ is the Gaussian barycenter for $W_{2}$ between the components $\mu_{0}^{k_{0}},\dots,\mu_{J-1}^{k_{J-1}}$, and $w^{*}$ is the optimal solution of (20). Moreover, this barycenter has less than $K_{0}+\dots+K_{J-1}-J+1$ non-zero coefficients. *

Proof 5.2**.**

*This follows directly from the proof of the previous propositions. The linear program (20) has $K_{0}+\dots+K_{J-1}-J+1$ affine constraints, and thus must have at least a solution with less than $K_{0}+\dots+K_{J-1}-J+1$ components. *

To conclude this section, it is important to emphasize that the problem of barycenters for the distance $MW_{2}$, as defined in (23), is completely different from

[TABLE]

Indeed, since $GMM_{d}(\infty)$ is dense in $\mathcal{P}_{2}(\mathbb{R}^{d})$ and the total cost on the right is continuous on $\mathcal{P}_{2}(\mathbb{R}^{d})$, the infimum in (25) is exactly the same as the infimum over $\mathcal{P}_{2}(\mathbb{R}^{d})$. Even if the barycenter for $W_{2}$ is not a mixture itself, it can be approximated by a sequence of Gaussian mixtures with any desired precision. Of course, these mixtures might have a very high number of components in practice.

5.3 Some examples

The previous propositions give us a very simple way to compute barycenters between Gaussian mixtures for the metric $MW_{2}$. For given mixtures $\mu_{0},\dots,\mu_{J-1}$, we first compute all the values $mmW_{2}(\mu_{0}^{k_{0}},\dots,\mu_{J-1}^{k_{J-1}})$ between their components (and these values can be computed iteratively, see Section 2.4.2) and the corresponding Gaussian barycenters $\nu_{k_{0}\dots k_{J-1}}$. Then we solve the linear program (20) to find $w^{*}$.

Figure 4 shows the barycenters between the following simple two dimensional mixtures

[TABLE]

where $I_{2}$ is the $2\times 2$ identity matrix. Each barycenter is a mixture of at most $K_{0}+K_{1}+K_{2}+K_{3}-4+1=11$ components. By thresholding the mixtures densities, this yields barycenters between 2-D shapes.

To go further, Figure 5 shows barycenters where more involved shapes have been approximated by mixtures of 12 Gaussian components each. Observe that, even if some of the original shapes (the star, the cross) have symmetries, these symmetries are not necessarily respected by the estimated GMM, and thus not preserved in the barycenters. This could be easily solved by imposing some symmetry in the GMM estimation for these shapes.

6 Using $MW_{2}$ in practice

6.1 Extension to probability distributions that are not GMM

Most applications of optimal transport involve data that do not follow a Gaussian mixture model and we can wonder how to make use of the distance $MW_{2}$ and the corresponding transport plans in this case. A simple solution is to approach these data by convenient Gaussian mixture models and to use the transport plan $\gamma$ (or one of the maps defined in the previous section) to displace the data.

Given two probability measures $\nu_{0}$ and $\nu_{1}$, we can define a pseudo-distance $MW_{K,2}(\nu_{0},\nu_{1})$ as the distance $MW_{2}({\mu}_{0},{\mu}_{1})$, where each ${\mu}_{i}$ ($i=0,1$) is the Gaussian mixture model with $K$ components which minimizes an appropriate “similarity measure” to $\nu_{i}$. For instance, if $\nu_{i}$ is a discrete measure $\nu_{i}=\frac{1}{J_{i}}\sum_{j=1}^{J_{i}}\delta_{x^{i}_{j}}$ in $\mathbb{R}^{d}$ , this similarity can be chosen as the opposite of the log-likelihood of the discrete set of points $\{x_{j}\}_{j=1,\dots,J_{i}}$ and the parameters of the Gaussian mixture can be infered thanks to the Expectation-Maximization algorithm. Observe that this log-likelihood can also be written

[TABLE]

If $\nu_{i}$ is absolutely continuous, we can instead choose ${\mu}_{i}$ which minimizes $\mathrm{KL}(\nu_{i},{\mu}_{i})$ among GMM of order $K$. The discrete and continuous formulations coincide since

[TABLE]

where $H(\nu_{i})$ is the differential entropy of $\nu_{i}$.

In both cases, the corresponding $MW_{K,2}$ does not define a distance since two different distributions may have the same corresponding Gaussian mixture. However, for $K$ large enough, their approximation by Gaussian mixtures will become different. The choice of $K$ must be a compromise between the quality of the approximation given by Gaussian mixture models and the affordable computing time. In any case, the optimal transport plan $\gamma_{K}$ involved in $MW_{2}({\mu}_{0},{\mu}_{1})$ can be used to compute an approximate transport map between $\nu_{0}$ and $\nu_{1}$.

In the experimental section, we will use this approximation for different data, generally with $K=10$.

6.2 A similarity measure mixing $MW_{2}$ and $KL$

In the previous paragraphs, we have seen how to use our Wasserstein-type distance $MW_{2}$ and its associated optimal transport plan on probability measures $\nu_{0}$ and $\nu_{1}$ that are not GMM. Instead of a two step formulation (first an approximation by two GMM, and second the computation of $MW_{2}$), we propose here a relaxed formulation combining directly $MW_{2}$ with the Kullback-Leibler divergence.

Let $\nu_{0}$ and $\nu_{1}$ be two probability measures on $\mathbb{R}^{d}$, we define

[TABLE]

where $\lambda>0$ is a parameter.

In the case where $\nu_{0}$ and $\nu_{1}$ are absolutely continuous with respect to the Lebesgue measure, we can write instead

[TABLE]

and $\widetilde{E_{K,\lambda}}(\nu_{0},\nu_{1})=E_{K,\lambda}(\nu_{0},\nu_{1})-\lambda H(\nu_{0})-\lambda H(\nu_{1}).$ Note that this formulation does not define a distance in general.

This formulation is close to the unbalanced formulation of optimal transport proposed by Chizat et al. in [10], with two differences: a) we constrain the solution $\gamma$ to be a GMM; and b) we use $\mathrm{KL}(\nu_{0},P_{0}\#\gamma)$ instead of $\mathrm{KL}(P_{0}\#\gamma,\nu_{0})$. In their case, the support of $P_{i}\#\gamma$ must be contained in the support of $\nu_{i}$. When $\nu_{i}$ has a bounded support, this constraint is quite strong and would not make sense for a GMM $\gamma$.

For discrete measures $\nu_{0}$ and $\nu_{1}$, when $\lambda$ goes to infinity, minimizing (26) becomes equivalent to approximate $\nu_{0}$ and $\nu_{1}$ by the EM algorithm and this only imposes the marginals of $\gamma$ to be as close as possible to $\nu_{0}$ and $\nu_{1}$. When $\lambda$ decreases, the first term favors solutions $\gamma$ whose marginals become closer.

Solving this problem (Equation (26)) leads to computations similar to those used in the EM iterations [4]. By differentiating with respect to the weights, means and covariances of $\gamma$, we obtain equations which are not in closed-form. For the sake of simplicity, we illustrate here what happens in one dimension.

Let $\gamma\in GMM_{2}(K)$ be a Gaussian mixture in dimension $2d=2$ with $K$ elements. We write

[TABLE]

We have that the marginals are given by the 1d Gaussian mixtures

[TABLE]

Then, to minimize, with respect to $\gamma$, the energy $E_{K,\lambda}(\nu_{0},\nu_{1})$ above, since the $\mathrm{KL}$ terms are independent of the $a_{k}$, we can directly take $a_{k}=\sigma_{0,k}\sigma_{1,k}$, and the transport cost term becomes

[TABLE]

Therefore, we have to consider the problem of minimizing the following “energy”:

[TABLE]

It can be optimized through a simple gradient descent on the parameters $\pi_{k}$, $m_{i,k}$, $\sigma_{i,k}$ for $i=0,1$ and $k=1,\ldots,K$. Indeed a simple calculus shows that we can write

[TABLE]

[TABLE]

[TABLE]

where we have introduced some auxilary empirical estimates of the variables given, for $i=0,1$ and $k=1,\ldots,K$, by

[TABLE]

[TABLE]

Automatic differenciation of $F$ can also be used in practice. At each iteration of the gradient descent, we project on the constraints $\pi_{k}\geq 0$, $\sigma_{i,k}\geq 0$ and $\sum_{k}\pi_{k}=1$.

On Figure 6, we illustrate this approach on a simple example. The distributions $\nu_{0}$ and $\nu_{1}$ are 1d discrete distributions, plotted as the red and blue histograms. On this example, we choose $K=3$ and we use automatic differenciation (with the torch.autograd Python library) for the sake of convenience. The red and blue plain curves represent the final distributions $P_{0}\#\gamma$ and $P_{1}\#\gamma$, for $\lambda$ in the set $\{10,2,1,0.5,0.1,0.01\}$. The behavior is as expected: when $\lambda$ is large, the KL terms are dominating and the distribution $\gamma$ tends to have its marginal fitting well the two distributions $\nu_{0}$ and $\nu_{1}$. Whereas, when $\lambda$ is small, the Wasserstein transport term dominates and the two marginals of $\gamma$ are almost equal.

6.3 From a GMM transport plan to a transport map

Usually, we need not only to have an optimal transport plan and its corresponding cost, but also an assignment giving for each $x\in\mathbb{R}^{d}$ a corresponding value $T(x)\in\mathbb{R}^{d}$. Let $\mu_{0}$ and $\mu_{1}$ be two GMM. Then, the optimal transport plan between $\mu_{0}$ and $\mu_{1}$ for $MW_{2}$ is given by

[TABLE]

It is not of the form $(\mathrm{Id},T)\#\mu_{0}$ (see also Figure 1 for an example), but we can however define a unique assignment of each $x$, for instance by setting

[TABLE]

where here $(X,Y)$ is distributed according to the probability distribution $\gamma$. Then, since the distribution of $Y|X=x$ is given by the discrete distribution

[TABLE]

we get that

[TABLE]

Notice that the $T_{mean}$ defined this way is an assignment that will not necessarily satisfy the properties of an optimal transport map. In particular, in dimension $d=1$, the map $T_{mean}$ may not be increasing: each $T_{k,l}$ is increasing but because of the weights that depend on $x$, their weighted sum is not necessarily increasing. Another issue is that $T_{mean}\#\mu_{0}$ may be “far” from the target distribution $\mu_{1}$. This happens for instance, in 1D, when $\mu_{0}=\mathcal{N}(0,1)$ and $\mu_{1}$ is the mixture of $\mathcal{N}(-a,1)$ and $\mathcal{N}(a,1)$, each with weight $0.5$. In this extreme case we even have that $T_{mean}$ is the identity map, and thus $T_{mean}\#\mu_{0}=\mu_{0}$, that can be very far from $\mu_{1}$ when $a$ is large.

Now, another way to define an assignment is to define it as a random assignment using the optimal plan $\gamma$. More precisely, for a fixed value $x$ we can define

[TABLE]

Observe that, from a mathematical point of view, we can define a random variable $T_{rand}(x)$ for a fixed value of $x$, or also a finite set of independent random variables $T_{rand}(x)$ for a finite set of $x$. But constructing and defining $T_{rand}$ as a stochastic process on the whole space $\mathbb{R}^{d}$ would be mathematically much more difficult (see [20] for instance).

Now, for any measurable set $A$ of $\mathbb{R}^{d}$ and any $x\in\mathbb{R}^{d}$, we can define the map $\kappa(x,A):=\mathbb{P}[T_{rand}(x)\in A],$ and we have

[TABLE]

It means that if the measure $T_{rand}$ could be defined everywhere, then “$T_{rand}\#\mu_{0}$”, would be equal in expectation to $\mu_{1}$.

Figure 7 illustrates these two possible assignments $T_{mean}$ and $T_{rand}$ on a simple example. In this example, two discrete measures $\nu_{0}$ and $\nu_{1}$ are approximated by Gaussian mixtures $\mu_{0}$ and $\mu_{1}$ of order $K$, and we compute the transport maps $T_{mean}$ and $T_{rand}$ for these two mixtures. These maps are used to displace the points of $\nu_{0}$. We show the result of these displacements for different values of $K$. We can see that depending on the configuration of points, the results provided by $T_{mean}$ and $T_{rand}$ can be quite different. As expected, the measure $T_{rand}\#\nu_{0}$ (well-defined since $\nu_{0}$ is composed of a finite set of points) looks more similar to $\nu_{1}$ than $T_{mean}\#\nu_{0}$ does. And $T_{rand}$ is also less regular than $T_{mean}$ (two close points can be easily displaced to two positions far from each other). This may not be desirable in some applications, for instance in color transfer as we will see in Figure 9 in the experimental section.

7 Two applications in image processing

We have already illustrated the behaviour of the distance $MW_{2}$ in small dimension. In the following, we investigate more involved examples in larger dimension. In the last ten years, optimal transport has been thoroughly used for various applications in image processing and computer vision, including color transfer, texture synthesis, shape matching. We focus here on two simple applications: on the one hand, color transfer, that involves to transport mass in dimension $d=3$ since color histograms are 3D histograms, and on the other hand patch-based texture synthesis, that necessitates transport in dimension $p^{2}$ for $p\times p$ patches. These two applications require to compute transport plans or barycenters between potentially millions of points. We will see that the use of $MW_{2}$ makes these computations much easier and faster than the use of classical optimal transport, while yielding excellent visual results. The codes of the different experiments are available through Jupyter notebooks on https://github.com/judelo/gmmot.

7.1 Color transfer

We start with the problem of color transfer. A discrete color image can be seen as a function $u:\Omega\rightarrow\mathbb{R}^{3}$ where $\Omega=\{0,\dots n_{r}-1\}\times\{0,\dots n_{c}-1\}$ is a discrete grid. The image size is $n_{r}\times n_{c}$ and for each $i\in\Omega$, $u(i)\in\mathbb{R}^{3}$ is a set of three values corresponding to the intensities of red, green and blue in the color of the pixel. Given two images $u_{0}$ and $u_{1}$ on grids $\Omega_{0}$ and $\Omega_{1}$, we define the discrete color distributions $\eta_{k}=\frac{1}{|\Omega_{k}|}\sum_{i\in\Omega_{k}}\delta_{u_{k}(i)}$, $k=0,1$, and we approximate these two distributions by Gaussian mixtures $\mu_{0}$ and $\mu_{1}$ thanks to the Expectation-Maximization (EM) algorithm333In practice, we use the scikit-learn implementation of EM with the kmeans initialization.. Keeping the notations used previously in the paper, we write $K_{k}$ the number of Gaussian components in the mixture $\mu_{k}$, for $k=0,1$. We compute the $MW_{2}$ map between these two mixtures and the corresponding $T_{mean}$. We use it to compute $T_{mean}(u_{0})$, an image with the same content as $u_{0}$ but with colors much closer to those of $u_{1}$. Figure 8 illustrates this process on two paintings by Renoir and Gauguin, respectively Le déjeuner des canotiers and Manhana no atua. For this experiment, we choose $K_{0}=K_{1}=10$. The corresponding transport map for $MW_{2}$ is relatively fast to compute (less than one minute with a non-optimized Python implementation, using the POT library [15] for computing the map between the discrete distributions of $10$ masses). We also show on the same figure $T_{rand}(u_{0})$, the result of the sliced optimal transport [25, 5], and the result of the separable optimal transport (i.e. on each color channel separately). Notice that the complete optimal transport on such huge discrete distributions (approximately 800000 Dirac masses for these $1024\times 768$ images) is hardly tractable in practice. As could be expected, the image $T_{rand}(u_{0})$ is much noisier than the image $T_{mean}(u_{0})$. We show on Figure 9 the discrete color distributions of these different images and the corresponding classes provided by EM (each point is assigned to its most likely class).

The value $K=10$ that we have chosen here is the result of a compromise. Indeed, when $K$ is too small, the approximation by the mixtures is generally too rough to represent the complexity of the color data properly. At the opposite, we have observed that increasing the number of components does not necessarily help since the corresponding transport map will loose regularity. For color transfer experiments, we found in practice that using around $10$ components yields the best results. We also illustrate this on Figure 10, where we show the results of the color transfer with $MW_{2}$ for different values of $K$. On the different images, one can appreciate how the color distribution gets closer and closer to the one of the target image as $K$ increases.

We end this section with a color manipulation experiment, shown on Figure 11. Four different images being given, we create barycenters for $MW_{2}$ between their four color palettes (represented again by mixtures of 10 Gaussian components), and we modify the first of the four images so that its color palette spans this space of barycenters. For this experiment (and this experiment only), a spatial regularization step is applied in post-processing [24] to remove some artifacts created by these color transformations between highly different images.

7.2 Texture synthesis

Given an exemplar texture image $u:\Omega\rightarrow R^{3}$, the goal of texture synthesis is to synthetize images with the same perceptual characteristics as $u$, while keeping some innovative content. The literature on texture synthesis is rich, and we will only focus here on a bilevel approach proposed recently in [16]. The method relies on the optimal transport between a continuous (Gaussian or Gaussian mixtures) distribution and a discrete distribution (distribution of the patches of the exemplar texture image). The first step of the method can be described as follows. For a given exemplar image $u:\Omega\rightarrow R^{3}$, the authors compute the asymptotic discrete spot noise (ADSN) associated with $u$, which is the stationary Gaussian random field $U:\mathbb{Z}^{2}\rightarrow\mathbb{R}^{3}$ with same mean and covariance as $u$, i.e.

[TABLE]

with $W$ a standard normal Gaussian white noise on $\mathbb{Z}^{2}$. Once the ADSN $U$ is computed, they extract a set $S$ of $p\times p$ sub-images (also called patches) of $u$. In our experiments, we extract one patch for each pixel of $u$ (excluding the borders), so patches are overlapping and the number of patches is approximately equal to the image size. The authors of [16] then define $\eta_{1}$ the empirical distribution of this set of patches (thus $\eta_{1}$ is in dimension $3\times p\times p$, i.e. $27$ for $p=3$) and $\eta_{0}$ the Gaussian distribution of patches of $U$, and compute the semi-discrete optimal transport map $T_{SD}$ from $\eta_{0}$ to $\eta_{1}$. This map $T_{SD}$ is then applied to each patch of a realization of $U$, and an ouput synthetized image $v$ is obtained by averaging the transported patches at each pixel. Since the semi-discrete optimal transport step is numerically very expensive in such high dimension, we propose to make use of the $MW_{2}$ distance instead. For that, we approximate the two discrete patch distributions of $u$ and $U$ by Gaussian Mixture models $\mu_{0}$ and $\mu_{1}$, and we compute the optimal map $T_{mean}$ for $MW_{2}$ between them. The rest of the algorithm is similar to the one described in [16]. Figure 12 shows the results for different choices of exemplar images $u$. In practice, we use $K_{0}=K_{1}=10$, as in color transfer, and $3\times 3$ color patches. The results obtained with our approach are visually very similar to the ones obtained with [16], for a computational time approximately 10 times smaller. More precisely, for instance for an image of size $256\times 256$, the proposed approach takes about $35$ seconds, whereas the semi-discrete approach of [16] takes about $400$ seconds. We are currently exploring a multiscale version of this approach, inspired by the recent [22].

8 Discussion and conclusion

In this paper, we have defined a Wasserstein-type distance on the set of Gaussian mixture models, by restricting the set of possible coupling measures to Gaussian mixtures. We have shown that this distance, with an explicit discrete formulation, is easy to compute and suitable to compute transport plans or barycenters in high dimensional problems where the classical Wasserstein distance remains difficult to handle. We have also discussed the fact that the distance $MW_{2}$ could be extended to other types of mixtures, as soon as we have a marginal consistency property and an identifiability property similar to the one used in the proof of Proposition 4. In practice, Gaussian mixture models are versatile enough to represent large classes of concrete and applied problems. One important question raised by the introduced framework and its generalization in Section 6.2 is how to estimate the mixtures for discrete data, since the obtained result will depend on the number $K$ of Gaussian components in the mixtures and on the parameter $\lambda$ that weights the data-fidelity terms. If the number of Gaussian components is chosen large enough, and covariances small enough, the transport plan for $MW_{2}$ will look very similar to the one of $W_{2}$, but at the price of a high computational cost. If, on the contrary, we choose a very small number of components (like in the color transfer experiments of Section 7.1), the resulting optimal transport map will be much simpler, which seems to be desirable for some applications.

Acknowledgments

We would like to thank Arthur Leclaire for his valuable assistance for the texture synthesis experiments.

Appendix: proofs

Density of $GMM_{d}(\infty)$ in $\mathcal{P}_{p}(\mathbb{R}^{d})$

Lemma** (Lemma 3.1).**

The set

[TABLE]

*is dense in $\mathcal{P}_{p}(\mathbb{R}^{d})$ for the metric $W_{p}$, for any $p\geq 1$. *

Proof 8.1**.**

The proof is adapted from the proof of Theorem 6.18 in [31] and given here for the sake of completeness.

*Let $\mu\in\mathcal{P}_{p}(\mathbb{R}^{d})$. For each $\epsilon>0$, we can find $r$ such that $\int_{B(0,r)^{c}}\|y\|^{p}d\mu(x)\leq\epsilon^{p}$, where $B(0,r)\subset\mathbb{R}^{d}$ is the ball of center [math] and radius $r$, and $B(0,r)^{c}$ denotes its complementary set in $\mathbb{R}^{d}$. The ball $B(0,r)$ can be covered by a finite number of balls $B(y_{k},\epsilon)$, $1\leq k\leq N$. Now, define $B_{k}=B(y_{k},\epsilon)\setminus\cup_{1\leq j<k}B(y_{j},\epsilon)$, all these sets are disjoint and still cover $B(0,r)$.

Define $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ on $\mathbb{R}^{d}$ such that*

[TABLE]

Then,

[TABLE]

and

[TABLE]

*which finishes the proof. *

Identifiability properties of Gaussian mixture models

Proposition** (Proposition 2).**

*The set of finite Gaussian mixtures is identifiable, in the sense that two mixtures $\mu_{0}=\sum_{k=1}^{K_{0}}\pi_{0}^{k}\mu_{0}^{k}$ and $\mu_{1}=\sum_{k=1}^{K_{1}}\pi_{1}^{k}\mu_{1}^{k}$, written such that all $\{\mu_{0}^{k}\}_{k}$ (resp. all $\{\mu_{1}^{j}\}_{j}$) are pairwise distinct, are equal if and only if $K_{0}=K_{1}$ and we can reorder the indexes such that for all $k$, $\pi_{0}^{k}=\pi_{1}^{k}$, $m_{0}^{k}=m_{1}^{k}$ and $\Sigma_{0}^{k}=\Sigma_{1}^{k}$.

This result is classical and the proof is also given here in the Appendix for the sake of completeness.

Proof 8.2**.**

This proof is an adaptation and simplification of the proof of Proposition 2 in [34]. First, assume that $d=1$ and that two Gaussian mixtures are equal:

[TABLE]

We start by identifying the Dirac masses from both sums, so only non-degenerate Gaussian components remain. Writing $\mu_{i}^{k}=\mathcal{N}(m_{i}^{k},(\sigma_{i}^{k})^{2})$, it follows that

[TABLE]

Now, define $k_{0}=\mathrm{argmax}_{k}\sigma_{0}^{k}$ and $j_{0}=\mathrm{argmax}_{j}\sigma_{1}^{j}$. If the maximum is attained for several values of $k$ (resp. $j$), we keep the one with the largest mean $m_{0}^{k}$ (resp. $m_{1}^{j}$). Then, when $x\rightarrow+\infty$, we have the equivalences

[TABLE]

Since the two sums are equal, these two terms must also be equivalent when $x\rightarrow+\infty$, which implies necessarily that $\sigma_{0}^{k_{0}}=\sigma_{1}^{j_{0}}$, $m_{0}^{k_{0}}=m_{1}^{j_{0}}$ and $\pi_{0}^{k_{0}}=\pi_{1}^{j_{0}}$. Now, we can remove these two components from the two sums and we obtain

[TABLE]

We can start over and show recursively that all components are equal.

For $d>1$, assume once again that two Gaussian mixtures $\mu_{0}$ and $\mu_{1}$ are equal, written as in Equation (28). The projection of this equality yields

[TABLE]

At this point, observe that for some values of $\xi$, some of these projected components may not be pairwise distinct anymore, so we cannot directly apply the result for $d=1$ to such mixtures. However, since the pairs $(m_{0}^{k},\Sigma_{0}^{k})$ (resp. $(m_{1}^{j},\Sigma_{1}^{j})$) are all distinct, then for $i=0,1$, the set

[TABLE]

is of Lebesgue measure [math] in $\mathbb{R}^{d}$. For any $\xi$ in $\mathbb{R}^{d}\setminus\Theta_{0}\cup\Theta_{1}$, the pairs $\{(\langle m_{0}^{k},\xi\rangle,\xi^{t}\Sigma_{0}^{k}\xi)\}_{k}$ (resp. $\{(\langle m_{1}^{j},\xi\rangle,\xi^{t}\Sigma_{1}^{j}\xi)\}_{j}$) are pairwise distinct. Consequently, using the first part of the proof (for $d=1$), we can deduce that $K_{0}=K_{1}$ and that

[TABLE]

where

[TABLE]

*Now, assume that the two sets $\{(\pi_{0}^{k},m_{0}^{k},\Sigma_{0}^{k})\}_{k}$ and $\{(\pi_{1}^{j},m_{1}^{j},\Sigma_{1}^{j})\}_{j}$ are different. Since each of these sets is composed of different triplets, it is equivalent to assume that there exists $k$ in $\{1,\dots K_{0}\}$ such that $(\pi_{0}^{k},m_{0}^{k},\Sigma_{0}^{k})$ is different from all triplets $(\pi_{1}^{j},m_{1}^{j},\Sigma_{1}^{j})$. In this case, the sets $\Xi_{k,j}$ for $j=1,\dots K_{0}$ are all of Lebesgue measure [math] in $\mathbb{R}^{d}$, which contradicts (30). We conclude that the sets $\{(\pi_{0}^{k},m_{0}^{k},\Sigma_{0}^{k})\}_{k}$ and $\{(\pi_{1}^{j},m_{1}^{j},\Sigma_{1}^{j})\}_{j}$ are equal. *