1

Optimal transport for Gaussian mixture models
Yongxin Chen, Tryphon T. Georgiou and Allen Tannenbaum

arXiv:1710.07876v2 [math.PR] 30 Jan 2018

Abstract
We present an optimal mass transport framework on the space of Gaussian mixture models. These models are
widely used in statistical inference. We treat such models as elements on the submanifold of probability densities
with Gaussian mixture structure and embed them to the density spaces equipped with Wasserstein metric. An
equivalent view relates to discrete measures on the space of Gaussian densities. Our method leads to a natural
way to compare, interpolate and average Gaussian mixture models, with low computational complexity. Different
aspects of our framework are discussed and several examples are presented for illustration. The method represents a
first attempt to study optimal transport problems for probability densities with specific structure that can be suitably
exploited.

I. I NTRODUCTION
A mixture model is a probabilistic model describing properties of populations with subpopulations.
Formally, it is a mixture distribution with each component representing a subpopulation. Mixture models
are widely used in statistics in detecting subgroups, inferring properties of subpopulations, and many other
areas [1]. An important case of mixture models is the so-called Gaussian mixture model (GMM), which
is simply a weighted average of several Gaussian distributions. Each Gaussian component stands for a
subpopulation. The Gaussian mixture model is commonly used in applications due to its mathematical
simplicity as well as efficient algorithms in inference (e.g., Expectation Maximization algorithm).
Optimal mass transport (OMT) is an old but very active research area dealing with probability densities.
Starting from the original formulation of Monge [2], and the relaxation of Kantorovich [3], and following
up with a sequence of important works [4], [5], [6], [7], [8], [9], the subject of OMT has become a
powerful tool in mathematics, physics, economics, engineering and biology etc [10], [11], [12], [13],
[14], [15], [16]. Benefitting from the development of alternative algorithms [17], [18], [19], [20], [21],
[22], [23], OMT has recently found applications in data science [24], [25]. Briefly, OMT deals with
problems of transporting masses from an initial distribution to a terminal distribution in a mass preserving
manner with minimum cost. When the unit cost is the square of the Euclidean distance, the OMT problem
induces an extremely rich geometry for probability densities. It endows a Riemannian metric on the space
of probability densities [26], [27], [28]. This geometry enables us to compare, interpolate and average
probability densities in a very natural way, which is in line with the needs in a range of applications.
Despite the inherent elegance and beauty OMT geometry, transport on the entire manifold of probability
densities is computationally expensive. However, in many applications, the probability densities often
have specific structure and may be parameterized [29]. Thus, we are motivated to study OMT on certain
submanifolds of probability densities. To retain the nice properties of OMT, herein, we seek an explicit
OMT framework on Gaussian mixture models. The extension to more general structured densities will be
a future research topic.
Y. Chen is with the Department of Electrical and Computer Engineering, Iowa State University, IA; email: yongchen@iastate.edu
T. T. Georgiou is with the Department of Mechanical and Aerospace Engineering, University of California, Irvine, CA; email:
tryphon@uci.edu
A. Tannenbaum is with the Departments of Computer Science and Applied Mathematics & Statistics, Stony Brook University, NY; email:
allen.tannenbaum@stonybrook.edu

2

This work is also partially motivated by problems in data science. Meaningful data are often of high
dimension and always have some structure. Thus, they are not densely distributed in the high dimensional
space. Instead, they usually live in a low dimensional sub-manifold. Besides, in many cases, the data
are sparsely distributed among subgroups. The difference between data within a subgroup is way less
significant than that between subgroups. For example, the differences between two dogs of the same
species are expected to be smaller than that of different species. In such applications, mixture models are
suitable, and therefore it is of importance to develop a mathematical framework that respects such data
structure.
II. BACKGROUND ON OMT
We now give a very brief overview of OMT theory. We only cover materials that are related to the present
work. We refer the reader to [27] for more details.
Consider two measures µ0 , µ1 on Rn with equal total mass. Without loss of generality, we take µ0 and
µ1 to be probability distributions. In the original formulation of OMT, a transport map
T : Rn → Rn : x 7→ T (x)
is sought that specifies where mass µ0 (dx) at x should be transported so as to match the final distribution
in the sense that T] µ0 = µ1 , i.e. µ1 is the “push-forward” of µ0 under T , meaning
µ1 (B) = µ0 (T −1 (B))
for every Borel set B in Rn . Moreover, the map should achieve a minimum cost of transportation
Z
c(x, T (x))µ0 (dx).
Rn

Here, c(x, y) represents the transportation cost per unit mass from point x to y. In this paper we focus
on the case when c(x, y) = kx − yk2 . To ensure finite cost, it is standard to assume that µ0 and µ1 live
in the space of probability densities with finite second moments, denoted by P2 (Rn ).
The dependence of the transportation cost on T is highly nonlinear and a minimum may not exist in general.
This fact complicated early analyses of the problem [27]. To circumvent this difficulty, Kantorovich
presented a relaxed formulation in 1942. In this, instead of seeking a transport map, one seeks a joint
distribution Π(µ0 , µ1 ) on Rn × Rn , referred to as “coupling” of µ0 and µ1 , so that the marginals along the
two coordinate directions coincide with µ0 and µ1 , respectively. Thus, in the Kantorovich formulation,
we solve
Z
inf
kx − yk2 π(dxdy).
(1)
π∈Π(µ0 ,µ1 )

Rn ×Rn

For the case where µ0 , µ1 are absolutely continuous with corresponding densities ρ0 and ρ1 , it is a standard
result that OMT (1) has a unique solution [4], [27], [28]. Moreover, the unique optimal transport T is the
gradient of a convex function φ, i.e.,
y = T (x) = ∇φ(x).
(2)
Having the optimal mass transport map T , as in (2), the optimal coupling is
π = (Id × T )] µ0 ,
where Id stands for the identity map. The square root of the minimum of the cost defines a Riemannian
metric on P2 (Rn ), known as the Wasserstein metric W2 [7], [26], [27], [28]. On this Riemannian-type
manifold, the geodesic curve connecting µ0 and µ1 is given by
µt = (Tt )] µ0 ,

Tt (x) = (1 − t)x + tT (x),

(3)

which is called displacement interpolation. It satisfies
W2 (µs , µt ) = (t − s)W2 (µ0 , µ1 ),

0 ≤ s < t ≤ 1.

(4)

3

A. Gaussian marginal distributions
When both of the marginals µ0 , µ1 are Gaussian distributions, the problem can be greatly simplified [30].
In fact, a closed-form solution exists. Denote the mean and covariance of µi , i = 0, 1 by mi and Σi ,
respectively. Let X, Y be two Gaussian random vectors associated with µ0 , µ1 , respectively. Then the cost
in (1) becomes
E{kX − Y k2 } = E{kX̃ − Ỹ k2 } + km0 − m1 k2 ,
(5)
where X̃ = X − m0 , Ỹ = Y − m1 are zero mean versions of X and Y . We minimize (5) over all the
possible Gaussian joint distributions between X and Y . This gives




Σ0 S
2
≥0 ,
(6)
min km0 − m1 k + trace(Σ0 + Σ1 − 2S) |
S T Σ1
S
with S = E{X̃ Ỹ T }. The constraint is semidefinite constraint, so the above problem is a semidefinite
programming (SDP). It turns out that the minimum is achieved by the unique minimizer in closed-form
1/2

1/2

1/2

−1/2

S = Σ0 (Σ0 Σ1 Σ0 )1/2 Σ0
with minimum value

1/2

1/2

W2 (µ0 , µ1 )2 = km0 − m1 k2 + trace(Σ0 + Σ1 − 2(Σ0 Σ1 Σ0 )1/2 ).
The consequent displacement interpolation µt is a Gaussian distribution with mean mt = (1 − t)m0 + tm1
and covariance

2
1/2
1/2 1/2
−1/2
−1/2
(1 − t)Σ0 + t(Σ0 Σ1 Σ0 )
Σ0 .
(7)
Σt = Σ0
The Wasserstein distance can be extended to singular Gaussian distributions by replacing the inverse by
the pseudoinverse †, which leads to
1/2

1/2

1/2

1/2

W2 (µ0 , µ1 )2 = km0 − m1 k2 + trace(Σ0 + Σ1 − 2Σ0 ((Σ0 )† Σ1 (Σ0 )† )1/2 Σ0 ).

(8)

In particular, when Σ0 = Σ1 = 0, we have
W2 (µ0 , µ1 ) = km0 − m1 k,
implying that the Wasserstein space of Gaussian distributions, denoted by G(Rn ), is an extension, at least
formally, of the Euclidean space Rn .
III. OMT FOR G AUSSIAN MIXTURE MODELS
A Gaussian mixture model is an important instance of mixture models, which are commonly used to
study properties of populations with several subgroups. Mathematically, a Gaussian mixture model is a
probability density consisting of several Gaussian components. Namely, it has the form
µ = p1 ν 1 + p2 ν 2 + · · · + pN ν N ,
where each ν k is a Gaussian distribution and p = (p1 , p2 , . . . , pN )T is a probability vector. Here the
finite number N stands for the number of components of µ. We denote the space of Gaussian mixture
distributions by M (Rn ).
As we have already seen in Section II-A, the displacement interpolation of two Gaussian distributions
remains Gaussian. This invariance, however, no longer holds for Gaussian mixtures. Yet, the mixture
models may contain some physical or statistical features that we may want to retain. This gives rise to

4

the following question we would like to address. How do we establish a geometry that inherits the nice
properties of OMT and in the meantime keeps the Gaussian mixture structure?
Our approach relies on a different way of looking at Gaussian mixture models. Instead of treating the
given mixture as a distribution on the Euclidean space Rn , we view it as a discrete distribution on the
Wasserstein space of Gaussian distributions G(Rn ). A Gaussian mixture distribution is equivalent to a
discrete measure, and therefore we can apply OMT theory to such discrete measures. We will see next
that this strategy retains the Gaussian mixture structure.
Let µ0 , µ1 be two Gaussian mixture models of the form
i Ni
µi = p1i νi1 + p2i νi2 + · · · + pN
i νi , i = 0, 1.

Here N0 maybe different to N1 . The distribution µi is equivalent to a discrete measure pi with supports
νi1 , νi2 , . . . , νiNi for each i = 0, 1. Our framework is built on the discrete OMT problem
X
min
c(i, j)π(i, j)
(9)
π∈Π(p0 ,p1 )

i,j

for these two discrete measures. Here Π(p0 , p1 ) denote the space of joint distributions between p0 and p1 .
The cost c(i, j) is taken to be the square of the Wasserstein metric on G(Rn ), that is,
c(i, j) = W2 (ν0i , ν1j )2 .
By standard linear programming theory, the discrete OMT problem (9) always has at least one solution.
Let π ∗ be a minimizer, and define
sX
c(i, j)π ∗ (i, j).
(10)
d(µ0 , µ1 ) =
i,j

Theorem 1: d(·, ·) defines a metric on M (Rn ).
Proof 1: Apparently, d(µ0 , µ1 ) ≥ 0 for any µ0 , µ1 ∈ M (Rn ) and d(µ0 , µ1 ) = 0 if and only if µ0 = µ1 .
We next prove the triangular inequality, namely,
d(µ0 , µ1 ) + d(µ1 , µ2 ) ≥ d(µ0 , µ2 )
for any µ0 , µ1 , µ2 ∈ M (Rn ). Denote the probability vector associated with µ0 , µ1 , µ2 by p0 , p1 , p2 respectively. The Gaussian components of µi is denoted by νij . Let π01 (π12 ) be the solution to (9) with marginals
µ0 , µ1 (µ1 , µ2 ). Define π02 by
X π01 (i, j)π12 (j, k)
π02 (i, k) =
.
pj1
j
Clearly, π02 is a joint distribution between p0 and p2 , namely, π02 ∈ Π(p0 , p2 ). It follows from direct
calculation
X
X π01 (i, j)π12 (j, k)
π02 (i, k) =
pj1
i
i,j
=
=

X pj π12 (j, k)
2

j
k
p2 .

pj1

5

π02 (i, k) = pi0 . Therefore,
sX
d(µ0 , µ2 ) ≤
π02 (i, k)W2 (ν0i , ν2k )2

Similarly, we have

P

k

i,k

=

≤

≤

=

v
uX
u
π01 (i, j)π12 (j, k)
t
W2 (ν0i , ν2k )2
j
p1
i,j,k
v
uX
u
π01 (i, j)π12 (j, k)
t
(W2 (ν0i , ν1j ) + W2 (ν1i , ν2k ))2
j
p
1
i,j,k
v
v
uX
uX
u
u
π01 (i, j)π12 (j, k)
π01 (i, j)π12 (j, k)
j 2
i
t
W2 (ν0 , ν1 ) + t
W2 (ν1j , ν2k )2
j
j
p1
p1
i,j,k
i,j,k
sX
sX
j 2
i
π01 (i, j)W2 (ν0 , ν1 ) +
π12 (j, k)W2 (ν1j , ν2k )2
i,j

j,k

= d(µ0 , µ1 ) + d(µ1 , µ2 ).
In the above, the second inequality is due to the fact W2 is a metric, and the third inequality is an
application of the Minkowski inequality.
A. Geodesic
A geodesic on M (Rn ) connecting µ0 and µ1 is given by
X
µt =
π ∗ (i, j)νtij ,

(11)

i,j

where νtij is the displacement interpolation (see (7)) between ν0i and ν1j .
Theorem 2:
d(µs , µt ) = (t − s)d(µ0 , µ1 ),

0 ≤ s < t ≤ 1.

(12)

Proof 2: For any 0 ≤ s ≤ t ≤ 1, we have
sX
d(µs , µt ) ≤
π ∗ (i, j)W2 (νsij , νtij )2
i,j

sX
π ∗ (i, j)W2 (ν0i , ν1j )2 = (t − s)d(µ0 , µ1 )
= (t − s)
i,j

where we have used the property (4) of W2 . It follows that
d(µ0 , µs ) + d(µs , µt ) + d(µt , µ1 ) ≤ sd(µ0 , µ1 ) + (t − s)d(µ0 , µ1 ) + (1 − t)d(µ0 , µ1 ) = d(µ0 , µ1 ).
On the other hand, by Theorem 1, we have
d(µ0 , µs ) + d(µs , µt ) + d(µt , µ1 ) ≥ d(µ0 , µ1 ).
Combining these two, we obtain (12).
We remark that µt is a Gaussian mixture model since it is a weighted average of the Gaussian distributions
νtij . Even though the solution to (9) is not unique in some instances, it is unique for generic µ0 , µ1 ∈
M (Rn ). Therefore, in most real applications, we need not worry about the uniqueness.

6

B. Relation between d and W2
We first note that we have
d(µ0 , µ1 ) ≥ W2 (µ0 , µ1 )
for any µ0 , µ1 ∈ M (Rn ). Equality holds when both µ0 and µ1 have only one Gaussian component. In
general, d > W2 . This is due to the fact that the restriction to the submanifold M (Rn ) induces suboptimality in the transport plan. Let γ(t), 0 ≤ t ≤ 1 be any piecewise smooth curve on M (Rn ) connecting
µ0 and µ1 . Define the Wasserstein length of γ by
X
LW (γ) =
sup
W2 (γtk , γtk+1 ),
0=t0 <t1 <···<ts =1

k

and natural length by
L(γ) =

X

sup
0=t0 <t1 <···<ts =1

d(γtk , γtk+1 ).

k

Then LW (γ) ≤ L(γ).
Using the metric property of d we get
d(µ0 , µ1 ) ≤ inf L(γ),
γ

where the minimization is taken over all the piecewise smooth curve on M (Rn ) connecting µ0 and µ1 .
In view of (12), we conclude
d(µ0 , µ1 ) = inf L(γ) ≥ inf LW (γ).
γ

γ

Therefore, it is unclear whetherd is the restriction of W2 to M (Rn ).
In general, d is a very good approximation of W2 if the variances of the Gaussian components are
small compared with the differences between the means. This may lead to an efficient algorithm to
approximate Wasserstein distance between two distributions with such properties. If we want to compute
the Wasserstein distance W2 (µ0 , µ1 ) between two distributions µ0 , µ1 ∈ M (Rn ), a standard procedure is
discretizing the densities first, and then solving a discrete OMT problem. Depending upon the resolution
of the discretization, the second step may become very costly. In contrast, to compute our new distance
d(µ0 , µ1 ), we need only to solve (9). When the number of Gaussian components of µ0 , µ1 are small, this
is extremely efficient.
IV. BARYCENTER OF G AUSSIAN MIXTURES
The barycenter [31] of L distributions µ0 , µ1 , . . . , µL is defined to be the minimizer of
L

1X
W2 (µ, µk )2 .
J(µ) =
L k=1
This resembles the average

1
(x1
L

(13)

+ x2 + · · · + xL ) of L points in the Euclidean space, which minimizes
1X
kx − xk k2 .
J(x) =
L k=1

The above definition can be generalized to the cost
min

µ∈P2 (Rn )

L
X
k=1

λk W2 (µ, µk )2 .

(14)

7

where λ = [λ1 , λ2 , . . . , λL ] is a probability vector. The existence and uniqueness of (14) has been
extensively studied in [31] where it is shown that under some mild assumptions, the solution exists
and is unique.
In the special case when all µk are Gaussian distributions, the barycenter remains Gaussian. In particular,
denoting the mean and covariance of µk as mk , Σk , then the barycenter has mean
m=

L
X

λk mk

(15)

λk (Σ1/2 Σk Σ1/2 )1/2 .

(16)

k=1

and covariance Σ solving
Σ=

L
X
k=1

A fast algorithm to get the solution of (16) is through the fixed point iteration [32]
!2
L
X
λk (Σ1/2 Σk Σ1/2 )1/2 Σ−1/2 .
(Σ)next = Σ−1/2
k=1

In practice, the iteration
(Σ)next =

L
X

λk (Σ1/2 Σk Σ1/2 )1/2

k=1

appears to also work. However, no convergence proof for the latter is known at present [31], [32].
For general distributions, the barycenter problem (14) is difficult to solve. It can be reformulated as a
multi-marginal optimal transport problem and is therefore convex. Recently several algorithms have been
proposed to solve (14) through entropic regularization [22]. However, due to the curse of dimensionality,
solving such a problem for dimension greater than 3 is still unrealistic. This is the case even for Gaussian
mixture models. What’s more, the Gaussian mixture structure is often lost when solving problem (14).
To overcome this issue for Gaussian mixtures, herein, we propose to solve a modified barycenter problem
min

µ∈M (Rn )

L
X

λk d(µ, µk )2 .

(17)

k=1

The optimization variable is restricted to be Gaussian mixture distribution and the Wasserstein distance
W2 is replaced by its relaxed version (10). Let µk be a Gaussian mixture distribution with Nk components,
n
k Nk
namely, µk = p1k νk1 + p2k νk2 + · · · + pN
k νk . If we view µ as a discrete measure on G(R ), then clearly, it
can only have support at the points (Gaussian distributions) of the form
argminν

L
X

λk W2 (ν, νkik )2

(18)

k=1

with νkik being any component of µk . As we discussed before, the optimal ν is Gaussian. Denote the set
of all such minimizers as {ν 1 , ν 2 , . . . , ν N }, then µ is equal to
µ = p1 ν 1 + p2 ν 2 + · · · + pN ν N ,

8

for some probability vector p = (p1 , p2 , . . . , pN )T . The number of element N is bounded above by
N1 N2 · · · NL . Finally, utilizing the definition of d(·, ·) we obtain an equivalent formulation of (17), which
reads as
min

π1 ≥0,··· ,πL ≥0
N
X
i=1
N1
X
j1 =1

Nk
L X
N X
X

λk ck (i, jk )πk (i, jk )

(19a)

k=1 i=1 jk =1

πk (i, jk ) = pjkk ,
π1 (i, j1 ) =

N2
X

∀1 ≤ k ≤ L, 1 ≤ jk ≤ Nk
π2 (i, j2 ) = · · · =

j2 =1

NL
X

πL (i, jL ),

(19b)
∀1 ≤ i ≤ N.

(19c)

jL =1

The cost
ck (i, j) = W2 (ν i , νkj )2

(20)

is the optimal transport cost from ν i to νkj . After solving the above linear programming problem (19), we
get the barycenter µ = p1 ν 1 + p2 ν 2 + · · · + pN ν N with
i

p =

N1
X

π1 (i, j)

j=1

for each 1 ≤ i ≤ N . We remark that our formulation is independent of the dimension of the underlying
space Rn . The dimension affects only the computation of the cost function (20) where a closed-form
(8) is available. The complexity of (19) relies on the numbers of components of the Gaussian mixtures
distributions {µk }. Therefore, our formulation is extremely efficiently for high dimensional Gaussian
mixtures with small number of components.
The difficulty of formulation (19) lies on the number N of components of the barycenter µ, which is
usually of order N1 N2 · · · NL . To overcome this issue, we can consider the barycenter problem for Gaussian
mixture with specified components. More specifically, given N Gaussian components ν 1 , ν 2 , . . . , ν N , we
would like to find a minimizer of the optimization problem (14) subject to the structure constraint that
µ = p1 ν 1 + p2 ν 2 + · · · + pN ν N
for some probability vector p = (p1 , p2 , . . . , pN )T . Note that ν k here doesn’t have to be of the form (18).
It can be any Gaussian distribution. Moreover, the number N can be chosen to be small. It turns out
this problem can be solved in exactly the same way. Clearly, a linear programming reformulation (19) is
straightforward.
V. N UMERICAL E XAMPLES
Several examples are provided to illustrate our framework in computing distance, geodesic and barycenter.
A. d vs W2
To demonstrate the difference between d and W2 , we investigate a simple example here. We choose µ0
to be a zero-mean Gaussian distribution with unit variance. The terminal distribution µ1 is set to be the
average of two unit variance Gaussian distributions, one with mean ∆ and the other one with mean −∆.
Clearly, d(µ0 , µ1 ) = ∆. Figure 1 depicts d(µ0 , µ1 ) and W2 (µ0 , µ1 ) for different ∆ values. As can be seen,
these two distances are not equivalent and d is always bounded below by W2 .

9

Fig. 1: d vs W2
B. Geodesic
We compare the displacement interpolation in standard OMT theory, and our proposed geodesic interpolation (see (11)) on M (Rn ). Consider the two Gaussian mixture models in Figure 2. Both of
them have two components; one in red, one in blue and the mixture model is in black. For both
marginals, the masses are equally distributed among the components. The means and covariances are
m10 = 0.5, m20 = 0.1, Σ10 = 0.01, Σ20 = 0.05 for µ0 , and m11 = 0, m21 = −0.35, Σ11 = 0.02, Σ21 = 0.02
for µ1 . Figure 3 depicts the interpolation results based on standard OMT and our method. As we can
see from the figures, the intermediate densities based on standard OMT interpolation lose the Gaussian
mixture structure. This is not the case for our method. The two Gaussian components of the interpolation
based on proposed method described are shown in Figure 4.
We have similar observations for an example on 2-dimensional spaces; see Figures 5-7. The two marginal
distributions are Gaussian mixtures shown in Figure 5. Figure 6 and Figure 7 are two interpolations, based
on OMT and our method, respectively. We can see that the Gaussian mixture structure is undermined in
Figure 6.

(a) µ0

(b) µ1

Fig. 2: Marginal distributions

10

(a) OMT

(b) our framework

Fig. 3: Interpolations

Fig. 4: Two Gaussian components of the interpolation
C. Barycenter
As shown in Figure 8, three Gaussian mixture distributions are given. The masses are equally distributed
among the components. The statistics for µ1 , µ2 , µ3 are (m11 = 0.5, m21 = 0.1, Σ11 = 0.01, Σ21 = 0.05),
(m12 = 0, m22 = −0.35, Σ12 = 0.02, Σ22 = 0.02) and (m13 = 0.4, m23 = −0.45, Σ13 = 0.025, Σ23 = 0.021)
respectively. We apply both our method and the traditional OMT theory to compute the barycenter. Two
sets of weights are considered and the results are displayed in Figure 9 and 10. Apparently, our method
gives better average. The Gaussian mixtures structure is damaged if traditional OMT is used.
VI. C ONCLUSION
In this note, we have defined a new optimal mass transport distance for Gaussian mixture models by
restricting ourselves to the submanifold of Gaussian mixture distributions. Consequently, the geodesic
interpolation utilizing this metric remains on the submanifold of Gaussian mixture distributions. On the
numerical side, computing this distance between two densities is equivalent to solving a linear programming problem whose number of variables grows linearly as the number of Gaussian components. This is
a huge reduction in computational cost compared with traditional OMT. Finally, when the covariances of
the components are small, our distance is a very good approximation of the standard OMT distance. The
extension to general mixture models or structural models will be an interesting direction in the future.

11

(a) µ0

(b) µ1

Fig. 5: Marginal distributions

Fig. 6: OMT Interpolation
ACKNOWLEDGEMENTS
This project was supported by AFOSR grants (FA9550-15-1-0045 and FA9550-17-1-0435), grants from
the National Center for Research Resources (P41- RR-013218) and the National Institute of Biomedical
Imaging and Bioengineering (P41-EB-015902), NCI grant (1U24CA18092401A1,) and NIA grant (R01
AG053991), and the Breast Cancer Research Foundation.
R EFERENCES
[1] G. McLachlan and D. Peel, Finite mixture models.

John Wiley & Sons, 2004.

[2] G. Monge, Mémoire sur la théorie des déblais et des remblais.

De l’Imprimerie Royale, 1781.

[3] L. V. Kantorovich, “On the transfer of masses,” in Dokl. Akad. Nauk. SSSR, vol. 37, no. 7-8, 1942, pp. 227–229.
[4] Y. Brenier, “Polar factorization and monotone rearrangement of vector-valued functions,” Communications on pure and applied
mathematics, vol. 44, no. 4, pp. 375–417, 1991.
[5] W. Gangbo and R. J. McCann, “The geometry of optimal transportation,” Acta Mathematica, vol. 177, no. 2, pp. 113–161, 1996.
[6] R. J. McCann, “A convexity principle for interacting gases,” Advances in mathematics, vol. 128, no. 1, pp. 153–179, 1997.

12

Fig. 7: Our Interpolation

(a) µ1

(b) µ2

(c) µ3

Fig. 8: Marginal distributions

[7] R. Jordan, D. Kinderlehrer, and F. Otto, “The variational formulation of the Fokker–Planck equation,” SIAM journal on mathematical
analysis, vol. 29, no. 1, pp. 1–17, 1998.
[8] J. D. Benamou and Y. Brenier, “A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem,” Numerische
Mathematik, vol. 84, no. 3, pp. 375–393, 2000.
[9] F. Otto and C. Villani, “Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality,” Journal of
Functional Analysis, vol. 173, no. 2, pp. 361–400, 2000.
[10] L. C. Evans and W. Gangbo, Differential equations methods for the Monge-Kantorovich mass transfer problem. American Mathematical
Soc., 1999, vol. 653.
[11] S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, “Optimal mass transport for registration and warping,” International Journal of
Computer Vision, vol. 60, no. 3, pp. 225–240, 2004.
[12] M. Mueller, P. Karasev, I. Kolesov, and A. Tannenbaum, “Optical flow estimation for flame detection in videos,” IEEE Transactions
on image processing, vol. 22, no. 7, pp. 2786–2797, 2013.
[13] Y. Chen, T. T. Georgiou, and M. Pavon, “On the relation between optimal transport and Schrödinger bridges: A stochastic control
viewpoint,” Journal of Optimization Theory and Applications, vol. 169, no. 2, pp. 671–691, 2016.
[14] Y. Chen, T. T. Georgiou, and M. Pavon, “Optimal transport over a linear dynamical system,” IEEE Transactions on Automatic Control,
vol. 62, no. 5, pp. 2137–2152, 2017.
[15] A. Galichon, Optimal Transport Methods in Economics.

Princeton University Press, 2016.

[16] Y. Chen, “Modeling and control of collective dynamics: From Schrödinger bridges to optimal mass transport,” Ph.D. dissertation,
University of Minnesota, 2016.
[17] S. Angenent, S. Haker, and A. Tannenbaum, “Minimizing flows for the Monge–Kantorovich problem,” SIAM journal on mathematical
analysis, vol. 35, no. 1, pp. 61–97, 2003.

13

(a) our method

(b) optimal transport

Fig. 9: Barycenters with λ = (1/3, 1/3, 1/3)

(a) our method

(b) optimal transport

Fig. 10: Barycenters with λ = (1/4, 1/4, 1/2)

[18] M. Cuturi, “Sinkhorn distances: Lightspeed computation of optimal transport,” in Advances in Neural Information Processing Systems,
2013, pp. 2292–2300.
[19] J. D. Benamou, B. D. Froese, and A. M. Oberman, “Numerical solution of the optimal transportation problem using the Monge-Ampere
equation,” Journal of Computational Physics, vol. 260, pp. 107–126, 2014.
[20] E. G. Tabak and G. Trigila, “Data-driven optimal transport,” Commun. Pure. Appl. Math. doi, vol. 10, p. 1002, 2014.
[21] E. Haber and R. Horesh, “A multilevel method for the solution of time dependent optimal transport,” Numerical Mathematics: Theory,
Methods and Applications, vol. 8, no. 01, pp. 97–111, 2015.
[22] J. D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré, “Iterative Bregman projections for regularized transportation problems,”
SIAM Journal on Scientific Computing, vol. 37, no. 2, pp. A1111–A1138, 2015.
[23] Y. Chen, T. Georgiou, and M. Pavon, “Entropic and displacement interpolation: a computational approach using the Hilbert metric,”
SIAM Journal on Applied Mathematics, vol. 76, no. 6, pp. 2375–2396, 2016.
[24] G. Montavon, K.-R. Müller, and M. Cuturi, “Wasserstein training of restricted boltzmann machines,” in Advances in Neural Information
Processing Systems, 2016, pp. 3718–3726.
[25] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein gan,” arXiv preprint arXiv:1701.07875, 2017.
[26] F. Otto, “The geometry of dissipative evolution equations: the porous medium equation,” Communications in Partial Differential
Equations, 2001.
[27] C. Villani, Topics in Optimal Transportation.

American Mathematical Soc., 2003, no. 58.

14

[28] C. Villani, Optimal Transport: Old and New.

Springer, 2008, vol. 338.

[29] S. I. Amari and H. Nagaoka, Methods of information geometry.

American Mathematical Soc., 2007, vol. 191.

[30] A. Takatsu, “Wasserstein geometry of gaussian measures,” Osaka Journal of Mathematics, vol. 48, no. 4, pp. 1005–1026, 2011.
[31] M. Agueh and G. Carlier, “Barycenters in the Wasserstein space,” SIAM Journal on Mathematical Analysis, vol. 43, no. 2, pp. 904–924,
2011.
[32] P. C. Álvarez-Esteban, E. del Barrio, J. Cuesta-Albertos, and C. Matrán, “A fixed-point approach to barycenters in Wasserstein space,”
Journal of Mathematical Analysis and Applications, vol. 441, no. 2, pp. 744–762, 2016.

