* [23] Frank S. de Boer, Dave Clarke, and Einar Broch Johnsen.  A complete guide to the future.  In Rocco De Nicola, editor, Proceedings of the 16th European Symposium on Programming, (ESOP 2007), volume 4421 of LNCS, pages 316–330, Berlin, Heidelberg, 2007. Springer-Verlag.
  * [24] Edsger W. Dijkstra.  Guarded commands, nondeterminacy and formal derivation of programs.  Commun. ACM, 18(8):453–457, 1975.
  * [25] Crystal Chang Din, Richard Bubel, and Reiner Hähnle.  KeY-ABS: A deductive verification tool for the concurrent modelling language ABS.  In Amy Felty and Aart Middeldorp, editors, Proc. 25th Intl. Conf. on Automated Deduction (CADE), Berlin, Germany, volume 9195 of LNCS, pages 517–526, Cham, 2015. Springer.
  * [26] Crystal Chang Din, Johan Dovland, Einar Broch Johnsen, and Olaf Owe.  Observable behavior of distributed systems: Component reasoning for concurrent objects.  Journal of Logic and Algebraic Programming, 81(3):227–256, 2012\.
  * [27] Crystal Chang Din, Reiner Hähnle, Einar Broch Johnsen, Violet Ka I Pun, and Silvia Lizeth Tapia Tarifa.  Locally abstract, globally concrete semantics of concurrent programming languages.  In Cláudia Nalon and Renate Schmidt, editors, Proc. 26th Intl. Conf. on Automated Reasoning with Tableaux and Related Methods, volume 10501 of LNCS, pages 22–43, Cham, Switzerland, September 2017. Springer.
  * [28] Crystal Chang Din and Olaf Owe.  A sound and complete reasoning system for asynchronous communication with shared futures.  J. Log. Algebraic Methods Program., 83(5–6):360–383, 2014.
  * [29] Crystal Chang Din and Olaf Owe.  Compositional reasoning about active objects with shared futures.  Formal Aspects of Computing, 27(3):551–572, 2015.
  * [30] Crystal Chang Din, Silvia Lizeth Tapia Tarifa, Reiner Hähnle, and Einar Broch Johnsen.  History-based specification and verification of scalable concurrent and distributed systems.  In Michael Butler, Sylvain Cochon, and Fatiha Zaïdi, editors, Proc. 17th International Conference on Formal Engineering Methods, ICFEM, Paris, volume 9407 of LNCS, pages 217–233, Cham, Switzerland, 2015\. Springer-Verlag.
  * [31] Johan Dovland, Einar Broch Johnsen, and Olaf Owe.  Verification of concurrent objects with asynchronous method calls.  In Proc. IEEE Intl. Conference on Software Science, Technology & Engineering(SwSTE’05), pages 141–150, Los Alamitos, CA, February 2005. IEEE Computer Society Press.
  * [32] Patrice Godefroid.  Model checking for programming languages using Verisoft.  In Peter Lee, Fritz Henglein, and Neil D. Jones, editors, Proc. 24th Symp. on Principles of Programming Languages (POPL), pages 174–186, New York, NY, 1997. ACM.
  * [33] Patrice Godefroid.  Test generation using symbolic execution.  In Deepak D’Souza, Telikepalli Kavitha, and Jaikumar Radhakrishnan, editors, IARCS Annual Conference on Foundations of Software Technology and Theoretical Computer Science, FSTTCS, Hyderabad, India, volume 18 of LIPIcs, pages 24–33, Schloss Dagstuhl, 2012. Leibniz-Zentrum fuer Informatik.
  * [34] Reiner Hähnle.  The Abstract Behavioral Specification language: A tutorial introduction.  In Marcello Bonsangue, Frank de Boer, Elena Giachino, and Reiner Hähnle, editors, International School on Formal Models for Components and Objects: Post Proceedings, volume 7866 of LNCS, pages 1–37, Cham, Switzerland, 2013. Springer.
  * [35] Reiner Hähnle and Marieke Huisman.  Deductive verification: from pen-and-paper proofs to industrial tools.  In Bernhard Steffen and Gerhard Woeginger, editors, Computing and Software Science: State of the Art and Perspectives, volume 10000 of LNCS, pages 345–373. Springer, Cham, Switzerland, 2019.
  * [36] Joseph Y. Halpern, Zohar Manna, and Ben C. Moszkowski.  A hardware semantics based on temporal intervals.  In Josep Díaz, editor, Automata, Languages and Programming, 10th Colloquium, Barcelona, Spain, volume 154 of LNCS, pages 278–291, Berlin, Heidelberg, 1983. Springer.
  * [37] David Harel, Dexter Kozen, and Jerzy Tiuryn.  Dynamic Logic.  Foundations of Computing. MIT Press, Boston, MA, USA, October 2000.
  * [38] Niklas Heidler.  Mechanization of LAGC semantics in Isabelle.  Bachelor thesis, Technical University of Darmstadt, Department of Computer Science, November 2021.
  * [39] Matthew Hennessy and Gordon D. Plotkin.  Full abstraction for a simple parallel programming language.  In Jirí Becvár, editor, Proc. 8th Symposium on the Mathematical Foundations of Computer Science, volume 74 of LNCS, pages 108–120, Heidelberg, 1979. Springer.
  * [40] Stefan Heule, K.Rustan M. Leino, Peter Müller, and AlexanderJ. Summers.  Abstract read permissions: Fractional permissions without the fractions.  In Roberto Giacobazzi, Josh Berdine, and Isabella Mastroeni, editors, Verification, Model Checking, and Abstract Interpretation, volume 7737 of LNCS, pages 315–334. Springer, Berlin Heidelberg, 2013.
  * [41] Carl Hewitt, Peter Boehler Bishop, and Richard Steiger.  A universal modular ACTOR formalism for artificial intelligence.  In Nils J. Nilsson, editor, Proc. 3rd Intl. Joint Conf. on Artificial Intelligence. Standford, CA, USA, pages 235–245, San Francisco, CA, USA, 1973. William Kaufmann.
  * [42] C. A. R. Hoare.  An axiomatic basis for computer programming.  Comm. of the ACM, 12(10):576–580, 583, October 1969.
  * [43] C. A. R. Hoare.  Communicating Sequential Processes.  International Series in Computer Science. Prentice Hall, Upper Saddle River, NJ, 1985.
  * [44] Gerard J. Holzmann.  The SPIN Model Checker.  Pearson Education, Boston, MA, USA, 2003.
  * [45] Michael Huth and Mark Ryan.  Logic in Computer Science: Modelling and Reasoning about Systems.  Cambridge University Press, Cambridge, UK, 2nd edition, 2004.
  * [46] Bart Jacobs, Frank Piessens, K. Rustan M. Leino, and Wolfram Schulte.  Safe concurrency for aggregate objects with invariants.  In Bernhard K. Aichernig and Bernhard Beckert, editors, Third IEEE Intl. Conf. on Software Engineering and Formal Methods (SEFM 2005), pages 137–147, Los Alamitos, CA, 2005. IEEE Computer Society.
  * [47] Alan Jeffrey and Julian Rathke.  A fully abstract may testing semantics for concurrent objects.  Theor. Comput. Sci., 338(1-3):17–63, 2005.
  * [48] Alan Jeffrey and Julian Rathke.  Java Jr: Fully abstract trace semantics for a core java language.  In Shmuel Sagiv, editor, Proc. 14th European Symposium on Programming (ESOP 2005), volume 3444 of LNCS, pages 423–438, Berlin Heidelberg, 2005. Springer.
  * [49] Einar Broch Johnsen, Reiner Hähnle, Jan Schäfer, Rudolf Schlatte, and Martin Steffen.  ABS: A core language for abstract behavioral specification.  In Bernhard K. Aichernig, Frank de Boer, and Marcello M. Bonsangue, editors, Proc. 9th International Symposium on Formal Methods for Components and Objects (FMCO 2010), volume 6957 of LNCS, pages 142–164, Heidelberg, 2011. Springer.
  * [50] C. B. Jones.  Development Methods for Computer Programmes Including a Notion of Interference.  PhD thesis, Oxford University, UK, June l981.
  * [51] Eduard Kamburjan.  Behavioral program logic.  In Serenella Cerrito and Andrei Popescu, editors, Proc. 28th International Conference on Automated Reasoning with Analytic Tableaux and Related Methods (TABLEAUX 2019), volume 11714 of LNCS, pages 391–408, Cham, Switzerland, 2019. Springer.
  * [52] Eduard Kamburjan.  Modular Verification of a Modular Specification: Behavioral Types as Program Logics.  PhD thesis, Darmstadt University of Technology, Germany, 2020.
  * [53] Eduard Kamburjan, Crystal Chang Din, Reiner Hähnle, and Einar Broch Johnsen.  Behavioral contracts for cooperative scheduling.  In Wolfgang Ahrendt, Bernhard Beckert, Richard Bubel, Reiner Hähnle, and Mattias Ulbrich, editors, Deductive Software Verification: Future Perspectives, volume 12345 of LNCS, pages 85–121. Springer, Cham, 2020.
  * [54] James C. King.  Symbolic execution and program testing.  Communications of the ACM, 19(7):385–394, July 1976.
  * [55] K.Rustan M. Leino, Peter Müller, and Jan Smans.  Verification of concurrent programs with chalice.  In Alessandro Aldini, Gilles Barthe, and Roberto Gorrieri, editors, Foundations of Security Analysis and Design V, volume 5705 of LNCS, pages 195–222. Springer, Berlin Heidelberg, 2009.
  * [56] Jayadev Misra and K. Mani Chandy.  Proofs of networks of processes.  IEEE Transactions on Software Engineering, 7(4):417–426, 1981\.
  * [57] John C. Mitchell.  Foundations for programming languages.  Foundation of computing. MIT Press, Boston, MA, USA, 1996.
  * [58] Peter D. Mosses.  Formal semantics of programming languages: — an overview —.  Electron. Notes Theor. Comput. Sci., 148(1):41–73, 2006.
  * [59] Wojciech Mostowski.  From explicit to implicit dynamic frames in concurrent reasoning for Java.  In Wolfgang Ahrendt, Bernhard Beckert, Richard Bubel, Reiner Hähnle, and Mattias Ulbrich, editors, Deductive Software Verification: Future Perspectives - Reflections on the Occasion of 20 Years of KeY, volume 12345 of LNCS, pages 177–203. Springer, Cham, 2020.
  * [60] Keiko Nakata and Tarmo Uustalu.  A Hoare logic for the coinductive trace-based big-step semantics of While.  Logical Methods in Computer Science, 11(1):1–32, 2015.
  * [61] Tobias Nipkow, Lawrence C. Paulson, and Markus Wenzel.  Isabelle/HOL — A Proof Assistant for Higher-Order Logic, volume 2283 of LNCS.  Springer, Berlin Heidelberg, 2002.
  * [62] Peter W. O’Hearn.  Resources, concurrency, and local reasoning.  Theoretical Computer Science, 375(1-3):271–307, 2007.
  * [63] Susan S. Owicki and David Gries.  An axiomatic proof technique for parallel programs I.  Acta Informatica, 6:319–340, 1976.
  * [64] Gordon D. Plotkin.  A structural approach to operational semantics.  J. Log. Algebr. Program., 60–61:17–139, 2004.
  * [65] A.P. Sistla and E.M. Clarke.  The complexity of propositional linear temporal logics.  Journal of the ACM, 32(3):733–749, 1985.
  * [66] Neelam Soundararajan.  Axiomatic semantics of communicating sequential processes.  ACM Transactions on Programming Languages and Systems, 6(4):647–662, 1984.
  * [67] Dominic Steinhöfel and Reiner Hähnle.  The trace modality.  In Alexandru Baltag and Luis S. Barbosa, editors, 2nd Intl. Workshop on Dynamic Logic: New Trends and Applications, volume 12005 of LNCS, pages 124–140, Cham, January 2020. Springer.
  * [68] Li-yao Xia, Yannick Zakowski, Paul He, Chung-Kil Hur, Gregory Malecha, Benjamin C. Pierce, and Steve Zdancewic.  Interaction trees: representing recursive and impure programs in Coq.  Proceedings of the ACM on Programming Languages, 4(POPL):1–32, Jan 2020.
  * [69] Marina Zaharieva-Stojanovski, Marieke Huisman, and Stefan Blom.  Verifying functional behaviour of concurrent programs.  In Proceedings of 16th Workshop on Formal Techniques for Java-like Programs, FTfJP’14, pages 4:1–4:6, New York, NY, USA, 2014. ACM.
  * [70] Job Zwiers.  Compositionality, Concurrency and Partial Correctness: Proof Theories for Networks of Processes, and Their Relationship, volume 321 of LNCS.  Springer-Verlag, Heidelberg, 1989.


# Optimal transport for vector Gaussian mixture models

Jiening Zhu, Kaiming Xu, Allen Tannenbaum J. Zhu and K. Xu are with the
Department of Applied Mathematics & Statistics, Stony Brook University, NY;
email<EMAIL_ADDRESS><EMAIL_ADDRESS>Tannenbaum is
with the Departments of Computer Science and Applied Mathematics & Statistics,
Stony Brook University, NY; email<EMAIL_ADDRESS>
###### Abstract

Vector Gaussian mixture models form an important special subset of vector-
valued distributions. Any physical entity that can mutate or transit among
alternative manifestations distributed in a given space falls into this
category. A key example is color imagery. In this note, we vectorize the
Gaussian mixture model and study different optimal mass transport related
problems for such models. The benefits of using vector Gaussian mixture for
optimal mass transport include computational efficiency and the ability to
preserve structure.

##  1 Introduction

Finite mixture models can describe a wide range of statistical phenomena. They
have been successfully applied to numerous fields including biology,
economics, engineering, and social sciences [15]. The first major use and
analysis of mixture models is perhaps due to the mathematician and
biostatistician Karl Pearson over 120 years ago who explicitly decomposed a
distribution into two normal distributions for characterizing non-normal
attributes of forehead to body length ratios in female shore crab populations
[16]. The literature on analyzing and applying mixture models is growing due
to their simplicity, versatility and flexibility. One of the most commonly
used mixture models is the Gaussian mixture model (GMM), which is a weighted
sum of Gaussian distributions.

Optimal mass transport (OMT) is a very important subject in mathematics,
originating with the French civil engineer and mathematician Gaspard Monge in
1781 [19, 20]. OMT considers the distance between two probability
distributions, which makes it a very powerful tool to analyze the geometry of
distributions. Its applications include but not limited to signal processing,
machine learning, computer vision, meteorology, statistical physics, quantum
mechanics, and network theory [1, 12, 17, 14]. Milestones of this subject
include Leonid Kantorovich [19, 20] who relaxed the original problem so that
it can be solved through simple linear programming, and Benamou and Brenier
[3] who introduced a computational fluid dynamics (CFD) approach to OMT. More
recent developments involve extensions of the theory to the vector-valued,
matrix-valued and unbalanced case [6, 5, 4, 9].

The problem that motivated the present work arose when the authors were
working with certain medical images. The object was to compute optimal mass
transport while preserving key structures. [7] studied OMT for GMM, which
however can only work on single layered data, i.e., gray scale images. The
need of working directly on the original color images with the potential of
capturing more information inspired us to generalize the OMT setting for one-
layer data to three-layer data, i.e., color images. More generally, in this
note, we develop optimal transport for vector Gaussian mixture models, which
can have any dimension and any general connection structures among the layers.
Furthermore, corresponding to unbalanced OMT, we also develop unbalanced
version for Gaussian mixture models.

In the present work, we will first give some background on GMM and OMT. After
that, we summarize some of the work of [7]. Based on that, we introduce two
different approaches for the vector case. We also look into unbalanced GMM
problem and conclude with some illustrative numerical results.

##  2 Gaussian mixture model

A Gaussian mixture model is one of the most important example of a mixture
model. Mathematically, a GMM is a probability distribution which is the
weighted sum of several Gaussian distributions in ${\mathbb{R}}^{N}$. Namely,
an $n$ Gaussians mixture model has the form

$\mu=p_{1}\nu_{1}+p_{2}\nu_{2}+\cdots+p_{n}\nu_{n}.$ (1)

Here

$\nu_{i}=\frac{1}{\sqrt{(2\pi)^{N}|\Sigma_{i}|}}\exp\\{-\frac{1}{2}(x-m_{i})^{T}\Sigma_{i}^{-1}(x-m_{i})\\},$
(2)

where $m_{i}\in{\mathbb{R}}^{N}$ is the mean and
$\Sigma_{i}\in{\mathbb{R}}^{N\times N}$ is the co-variant matrix. Further,

$\sum_{i=1}^{N}p_{i}=1,$ (3)

so that $\mu$ is a probability distribution.

We denote the set of all the GMMs in ${\mathbb{R}}^{N}$ by
$\mathcal{G}({\mathbb{R}}^{N})$. It is a dense subset of all the probability
distributions in the sense of the weak$*$ topology. Thus one can use GMM to
fit a distribution with arbitrarily small error, which makes GMM a very useful
model.

##  3 Optimal mass transport

In this section we sketch the basics of optimal mass transport.

The original formulation of OMT due to Monge [19, 20] may be expressed as
follows:

$\inf_{T}\\{\int_{E}c(x,T(x))\rho_{0}(x)dx\ |\ T_{\\#}\rho_{0}=\rho_{1}\\},$
(4)

where $c(x,y)$ is the cost of moving unit mass from $x$ to $y$, $T$ is the
transport map, and $\rho_{0},\rho_{1}$ are two probability distributions
defined on $E$, a subdomain of $\mathbb{R}^{n}$. $T_{\\#}$ denotes the push-
forward of $T$.

As pioneered by Leonid Kantorovich [19, 20], the Monge formulation of OMT may
be relaxed replacing transport maps $T$ by couplings $\pi$:

$\inf_{\pi\in\Pi(\rho_{0},\rho_{1})}\int_{E}c(x,y)\pi(dx,dy),$ (5)

where $\Pi(\rho_{0},\rho_{1})$ denotes the set of all the couplings between
$\rho_{0}$ and $\rho_{1}$ (joint distributions whose two marginal
distributions are $\rho_{0}$ and $\rho_{1}$).

The discrete Kantorovich form may be written as follows:

$\min_{\pi\in\Pi(\rho_{0},\rho_{1})}\sum_{i,j}c(x,y)\pi(i,j),$ (6)

where $\Pi(\rho_{0},\rho_{1})$ is the set of matrices
$\\{\pi\in{\mathbb{R}}_{+}^{m\times n}|\pi\vec{1}_{n}=\rho_{0},\
\pi^{T}\vec{1}_{m}=\rho_{1}\\}$, $\vec{1}_{m}$ and $\vec{1}_{n}$ are vectors
all of whose entries are 1 of length $m$ and $n$, respectively. $\min$ is used
instead of $\inf$ here because it is a finite dimensional minimizing problem
on a closed nonempty set
($\rho_{0}\otimes\rho_{1}^{T}\in\Pi(\rho_{0},\rho_{1})$).

One may show that for $c(x,y)=\|x-y\|^{2}$ (square of distance function), the
Kantorovich and Monge formulations are equivalent; see [19, 20] and the
references therein. Moreover for $c(x,y)=||x-y||^{2}$, the specific infimum is
called _Wasserstein-2 distance_ ($\mathcal{W}_{2}$).

##  4 Optimal mass transport for Gaussian mixture models

We can consider GMM as a general distribution and directly compute optimal
mass transport using the exact definition from last section. However, there
are two benefits of considering GMM as a special case. First, for multiple
Gaussians, there is no guarantee of preserving that structure on its geodesic
path. Second, for data with certain structure, expressing the data in the GMM
form can be considered as a way of dimension reduction, which can greatly save
computational cost.

###  4.1 OMT between Gaussian distributions

For two Gaussian distributions $\mu_{i},\ i=0,1$ whose means and covariances
are $m_{i}$ and $\Sigma_{i}$, it is standard [19, 20] that the
$\mathcal{W}_{2}$ distance between $\mu_{0}$ and $\mu_{1}$ has a closed form
solution:

$\mathcal{W}_{2}(\mu_{0},\mu_{1})^{2}=||m_{0}-m_{1}||^{2}+{\rm
trace}(\Sigma_{0}+\Sigma_{1}-2(\Sigma_{0}^{1/2}\Sigma_{1}\Sigma_{0}^{1/2})^{1/2}).$
(7)

For each $t\in[0,1]$, the distribution $\mu_{t}$ on the geodesic path is a
Gaussian whose mean and covariance matrix are defined as follows:

$\displaystyle m_{t}=(1-t)m_{0}+tm_{1}$ (8)
$\displaystyle\Sigma_{t}=\Sigma_{0}^{-1/2}((1-t)\Sigma_{0}+t(\Sigma_{0}^{1/2}\Sigma_{1}\Sigma_{0}^{1/2})^{1/2})^{2}\Sigma_{0}^{-1/2}.$
(9)

###  4.2 OMT between GMMs

Let $\mu_{0},\mu_{1}$ be two Gaussian mixture models of the form

$\mu_{i}=p_{i}^{1}\nu_{i}^{1}+p_{i}^{2}\nu_{i}^{2}+\cdots+p_{i}^{n_{i}}\nu_{i}^{n_{i}},\
i=0,1.$

The distance between $\mu_{0},\mu_{1}$ is defined as

$d(\mu_{0},\mu_{1})^{2}=\min_{\pi\in\Pi(p_{0},p_{1})}\sum_{i,j}c(i,j)\pi(i,j),$
(10)

where

$c(i,j)=\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}.$ (11)

As $\nu_{0}^{i}$ and $\nu_{1}^{j}$ are Gaussian distributions, the
$\mathcal{W}_{2}$ distance may be computed as in (7). In [7], it is proven
that $d(\cdot,\cdot)$ is indeed a metric on $\mathcal{G}({\mathbb{R}}^{N})$.
Furtehr, the geodesic on $\mathcal{G}({\mathbb{R}}^{n})$ connecting $\mu_{0}$
and $\mu_{1}$ is given by

$\mu_{t}=\sum_{i,j}\pi^{*}(i,j)\nu_{t}^{ij},$ (12)

where $\nu_{t}^{ij}$ is the displacement interpolation in (8) between
$\nu_{0}^{i}$ and $\nu_{1}^{j}$. $\pi^{*}(\cdot,\cdot)$ is the optimal
solution of (10).

##  5 Vector GMM

In this section, we extend the definition of GMM to the vector-valued case,
based on which we will formulate generalizations of the work of [7] in the
following sections.

###  5.1 Vector-valued distributions

A _vector-valued distribution_ $\rho$ on ${\mathbb{R}}^{n}$ is a map from
${\mathbb{R}}^{N}$ to ${\mathbb{R}}^{M}_{+}$ such that

$\sum_{i=1}^{M}\int_{{\mathbb{R}}^{N}}\rho_{i}(x)dx=1.$

It may represent a physical entity that can mutate or be transported among
several alternative manifestations with certain connections among its $M$
channels. This may be determined by a connected graph $G=(V,E)$ which has $M$
nodes and the edges determining the connections. Thus, $\rho$ may also be
understood as a general distribution on ${\mathbb{R}}^{N}\times G$, where
$V=\\{1,2,\cdots,M\\}$ and $E$ defines the connections among channels.

Figure 1: An example of a vector-valued distribution in ${\mathbb{R}}^{2}$
with 5 channels and a specific graph structure. The distribution takes value
only within each channel and its total sum is 1.

###  5.2 Vector GMMs as a subset of vector-valued distributions

Vector-valued GMMs are those vector-valued distributions such that the
distribution in each layer is a weighted sum of Gaussians and the weights of
the Gaussians sum up to 1. Formally,

$\rho=p^{1}\nu^{1}\vec{\delta}_{q^{1}}+p^{2}\nu^{2}\vec{\delta}_{q^{2}}+\cdots+p^{n}\nu^{n}\vec{\delta}_{q^{n}},$
(13)

where $\vec{\delta}_{k}$ is a column vector which is the $k^{\rm th}$ column
of the $M$ by $M$ identity matrix and $q^{i}$ is the index of channel where
the $i^{\rm th}$ Gaussian lies in. We will always assume that the latter is a
probability distribution, i.e.,

$\sum_{i=1}^{n}p^{i}=1.$ (14)

##  6 Generalization of the OMT framework in [7] to vector GMM

Consider two vector GMMs $\rho_{0}$ and $\rho_{1}$:

$\displaystyle\rho_{i}=$ $\displaystyle
p_{i}^{1}\nu_{i}^{1}\vec{\delta}_{q_{i}^{1}}+p_{i}^{2}\nu_{i}^{2}\vec{\delta}_{q_{i}^{2}}+\cdots+p_{i}^{n_{i}}\nu_{i}^{n_{i}}\vec{\delta}_{q_{i}^{n_{i}}},\
i=0,1.$

We want to compute an OMT based distance and a displacement interpolation
between these two vector valued distributions under the assumption that the
vector GMM structure is preserved along the interpolation path. In short, we
want to construct the same framework as in [7], but replace scalar-valued GMM
with vector-valued GMM.

As above, $\Pi(p_{0},p_{1})$ denotes the set of joint probabilities with given
marginals $p_{0}$ and $p_{1}$. Given a graph structure, the most
straightforward idea is to require only certain parts of $\Pi$ to be nonzero,
namely only when the the source and target Gaussians are in the same channel
or when they are located in adjacent channels. A more detailed description is
given in Appendix A.

However, this natural and straightforward generalization may not admit a
solution. Indeed, the newly added constraints on $\Pi$ may not work for
general graph structures. The only other choice left to be modified in (10) is
the cost matrix $c(\cdot,\cdot)$.

We propose two different approaches. Both of them modify the cost matrix
$c(\cdot,\cdot),$ but from two totally different angles.

##  7 Approach 1: modify cost matrix

In light of the limitations listed in the Appendix A.1 of the approach just
outlined, we propose a different method. As mentioned, since adding
restrictions on $\Pi$ is a delicate issue, we modify the cost matrix
$c(\cdot,\cdot)$ instead.

Since we know how to deal with OMT on ${\mathbb{R}}^{n}$, and noting that it
is also easy to compute the shortest path on a connected graph using Bellman-
Ford algorithm [2, 11] or Dijkstra’s algorithm [10], our idea is to consider
${\mathbb{R}}^{n}$ and $G$ separately, and combine them in the Kantorovich
step (10) as a new cost term.

More precisely, as the one-channel GMM case, using the Kantorovich step (10),
we take

$c_{1}(i,j)=\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})+\gamma\tilde{d}_{G}(q_{0}^{i},q_{1}^{j}),$
(15)

where $\tilde{d}_{G}(\cdot,\cdot)$ is the length of the shortest path on the
graph $G$ and $\gamma$ is a weight parameter to control how much we allow the
cross-channel transport.

With our new $c_{1}(i,j)$ defined, we can compute the distance as

$d_{V_{1}}(\rho_{0},\rho_{1})=\min_{\pi\in\Pi(p_{0},p_{1})}\sum_{i,j}c_{1}(i,j)\pi(i,j).$
(16)

As in the one-channel case, we can prove $d_{V_{1}}(\cdot,\cdot)$ is indeed a
metric on $\mathcal{G}({\mathbb{R}}^{N}\times G)$.

######  Theorem 1.

$d_{V_{1}}(\cdot,\cdot)$ defines a metric on
$\mathcal{G}({\mathbb{R}}^{N}\times G)$

###### Proof.

See Appendix B. ∎

The _displacement interpolation_ may be defined as follows:

$\rho_{t}=\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)\nu_{t}^{ij}\vec{\delta}_{path_{G}(q_{0}^{i},q_{1}^{j},t)},$
(17)

where $\tilde{\pi}_{1}^{*}(\cdot,\cdot)$ denotes the optimal solution of (16),
$\nu_{t}^{ij}$ is the displacement interpolation in (8) between $\nu_{0}^{i}$,
and $\nu_{1}^{j}$, $path_{G}(\cdot,\cdot,t)$ is the path interpolation between
two nodes of $G$. The value of $\vec{\delta}$ between two nodes is taken as
the weighted sum of the two $\vec{\delta}$ vectors of the two nodes.

For example, suppose that the shortest path between the two nodes $node_{1}$
and $node_{2}$ is $node_{1}\rightarrow node_{3}\rightarrow node_{4}\rightarrow
node_{2}$. Then at $t=0.5$ the interpolation point lies exactly in the middle
of points 3 and 4, which is to say
$path_{G}(node_{1},node_{2},0.5)=0.5*node_{3}+0.5*node_{4}$. Further,
$\vec{\delta}_{path_{G}(node_{1},node_{2},0.5)}=\vec{\delta}_{0.5*node_{3}+0.5*node_{4}}=0.5*\vec{\delta}_{3}+0.5*\vec{\delta}_{4}$.

######  Theorem 2.

The displacement interpolation (17) is a geodesic in the sense of
$d_{V_{1}}(\cdot,\cdot)$

###### Proof.

See Appendix C. ∎

##  8 Approach 2: Continuous version of $\mathcal{G}({\mathbb{R}}^{N}\times
G)$

For vector-valued GMMs, one cannot directly apply the same OMT framework as in
the scalar case [7] since the last index is taken discretely. Thus we will
generalize the framework by making the last index continuous as well. The
basic idea is to consider a continuous problem and view the vector-valued
distribution as a projection of the continuous solution onto original discrete
space.

More precisely, we propose to extend each point on the edges of the graph
instead of only taking values on vertices of graph $G$. Moreover, we extend
each edge half way from both ends, so that newly added points are centered at
original vertices of the graph. Thus, we consider the following point set of a
continuous version of the graph $G$ given as follows:

$G^{c}=\\{u+t(w-u)|u,w\in V(G),u\sim w,t\in[-0.5,0.5]\\}.$ (18)

Here $u,w\in V(G)$ are taken as abstract vertices, not as integers. In
addition, we assign a length to each edge, $\gamma$, so that we are able to
perform integration on that set. (We may consider the use of nonuniform edge
lengths, in case we are given specific edge weights.)

So the continuous structure we are going to employ is ${\mathbb{R}}^{N}\times
G^{c}$. Though it is possible to realize this whole structure in a Euclidean
space ${\mathbb{R}}^{N+K}$, one finds an unnecessarily complicated topology
space in general, as the two examples shown in Figures 2 3 indicate.

Figure 2: Consider a 3-vector distribution in ${\mathbb{R}}^{2}$ (may be used
for color image data). There are essentially two ways to connect the R,G and B
channels. If the graph is given by $R\rightarrow G\rightarrow B,$ then its
continuous version becomes a cube. On the other hand if additionally, R is
also connected to B, then the continuous version is like a hollow cylinder,
which is already a rather complicated space to employ in an OMT framework.
Figure 3: Consider a 4-vector distribution in ${\mathbb{R}}$ given by four 1D
edges ($[0,1]$) with each edge connected to to the other three edges making a
Klein bottle with a hole on its surface, illustrated by the shape shown on the
right. Computing optimal transport on its surface in general would be quite
challenging.

In fact, we do not need to realize the global structure of the complicated
space ${\mathbb{R}}^{N}\times G^{c}$ as a whole. Instead, we can just consider
the local structure. A natural and simple way to do that is to impose a
manifold structure that we will now elucidate. (For the convenience of the
reader we have included some background definitions concerning manifolds in
Appendix D ). We denote the manifold as $\mathcal{M}$.

In order to define $\mathcal{M}$, we need to specify its atlas:

$A=\\{{\mathbb{R}}^{N}\times p|p\in[G^{c}]_{0}\\},$ (19)

where $[G^{c}]_{0}$ is a subset of all continuous paths on $G^{c}$ which have
no cycles (no recurring vertices of $G$ on the path). We can characterize the
charts as we stack layers like bricks where we follow the order of the path on
$G^{c}$. It’s clear to see that each $p$ is homeomorphic to ${\mathbb{R}}$, so
that each chart is homeomorphic to ${\mathbb{R}}^{N+1}$.

We want to define a distribution on $\mathcal{M}$ in a way that the original
distribution is the projection of each layer’s range. The projection is
defined as the integral of the last index:

$P_{u}(f(x,z))=\sum_{w\sim u}\int_{-0.5}^{0.5}f(x,u+t(w-u))dt,$ (20)

where $P_{u}(\cdot)$ is the projection of the range of layer $u$, and $f$ is a
distribution on the manifold. The integral range of the last index is the
intersection of a ball centered at $u$ which has half-edge radius with $G^{c}$
(layer $u$’s range). Note that for different $w$’s that are connected to $u$,
the ranges are like different orbits centered at $u$.

One of the simplest choices for lifting the original distribution to the
manifold is a “Gaussian cylinder,” i.e., a product of a Gaussian distribution
and a uniform distribution within the range of the layer. Thus, we thicken
each Gaussian.

Figure 4: Left hand side is one of the layers of vector GMM. Right hand side
is the chart centered at that layer. Gaussians in the original layer become
”Gaussian cylinders” on the manifold.

If a layer has more than one edge connected to it, then the original
distribution may be lifted to multiple “Gaussian cylinders” (located at all
the possible orbits that are centered at the given layer) with combined
weights. Notice that though those “Gaussian cylinders” project to be the same
vector-valued distribution within that layer, they may have different
potentials to transport to different directions on the graph.

Let us briefly summarize the optimal transport problem we are going to solve
on the manifold $\mathcal{M}$ with the notions we just introduced. Given the
projection of each layer’s range for source (starting) and target (terminal)
distributions, we want to find corresponding source and target distributions
on $\mathcal{M}$ such that the transport cost is optimally low.

As before, we first consider the sub-problem where the starting and terminal
vector-valued distributions are two Gaussian-distributions which may be
located on different layers.

Figure 5: When we consider the transport map from the red Gaussian
distribution to the green Gaussian, we consider the transport problem on all
the charts that covers both Gaussian distributions. The above figure gives two
of the charts.

######  Theorem 3.

For any two Gaussian-cylinder-shaped distributions whose projections of each
layer’s range are simple Gaussian distributions denoted by $\nu_{0}$ and
$\nu_{1},$ located on layers $u$ and $w$, respectively, the optimal transport
$\mathcal{W}_{2}$ distance between them on $\mathcal{M}$ is given by
$d_{\mathcal{M}}=\mathcal{W}_{2}(\nu_{0},\nu_{1})^{2}+\gamma\tilde{d}_{G}(u,w)^{2}$

###### Proof.

We consider all the couplings on the manifold $\mathcal{M}$ denoted by
$\Pi(\mathcal{M})$. More precisely, we consider all the possible transports on
the charts in $A$ which can cover the supports of both Gaussian cylinders
lifted from the two original Gaussians, that is, all the charts in
$\\{{\mathbb{R}}^{N}\times p|p\in[G^{c}]_{0}^{uw}\\}$ where $[G^{c}]_{0}^{uw}$
denotes the subset of $[G^{c}]_{0}$ of those paths contain both layer $u$ and
layer $w$. With this definition, we write down the optimization problem
explicitly:

$\displaystyle d_{\mathcal{M}}=$
$\displaystyle\inf_{\pi\in\Pi(\mathcal{M})}\int_{\mathcal{M}\times\mathcal{M}}||\tilde{x}-\tilde{y}||^{2}\pi(d\tilde{x},d\tilde{y})$
$\displaystyle=$
$\displaystyle\inf_{p\in[G^{c}]_{0}^{uw}}\inf_{\pi\in\Pi^{p}({\mathbb{R}}^{N+1})}\int_{{\mathbb{R}}\times{\mathbb{R}}}\int_{{\mathbb{R}}^{N}\times{\mathbb{R}}^{N}}||x-y||^{2}+|z_{1}-z_{2}|^{2}\pi(dxdz_{1},dydz_{2})$

Here, $\Pi^{p}({\mathbb{R}}^{N+1})$ denotes the couplings in
${\mathbb{R}}^{N+1}$ (which is homeomorphic to ${\mathbb{R}}^{N}\times p$)
whose marginals are the source and target Gaussian cylinders, respectively.

Further, because the special structure of “Gaussian cylinders,” the first $N$
indices and the last index may be treated separately. If we denote by
$\Pi_{1}({\mathbb{R}}^{N})$ the set of couplings in ${\mathbb{R}}^{N}$ for
which the two marginals are the original source and target Gaussians (which
does not depend on the path $p$), and denote by $\Pi^{p}_{2}({\mathbb{R}})$
the set of couplings whose two marginals are two uniform distributions located
in their corresponding layers, the distance expression can be divided into two
parts:

$\displaystyle\inf_{p\in[G^{c}]_{0}^{uw}}\inf_{\pi_{1}\in\Pi_{1}({\mathbb{R}}^{N})}\int_{{\mathbb{R}}^{N}\times{\mathbb{R}}^{N}}||x-y||^{2}\pi_{1}(dx,dy)+\inf_{p\in[G^{c}]_{0}^{uw}}\inf_{\pi_{2}\in\Pi^{p}_{2}({\mathbb{R}})}\int_{{\mathbb{R}}\times{\mathbb{R}}}|z_{1}-z_{2}|^{2}\pi_{2}(dz_{1},dz_{2})$

The second term is a simple 1D optimal transport problem between two uniform
distributions which are centered at $u$ and $w$, respectively, with the same
radius of thickness of each layer. The optimal transport distance between them
is simply the distance between their respective centers, which is easy to
calculate since both two centers are located on the path $p$. To be specific,
the distance is the length of the path connecting $u$ and $w$ times the
thickness of each layer. Hence,

$\displaystyle d_{\mathcal{M}}=$
$\displaystyle\inf_{\pi_{1}\in\Pi_{1}({\mathbb{R}}^{N})}\int_{{\mathbb{R}}^{N}\times{\mathbb{R}}^{N}}||x-y||^{2}\pi_{1}(dx,dy)+\inf_{p\in[G^{c}]_{0}^{uw}}\Delta_{p}z^{2}$
$\displaystyle=$ $\displaystyle\quad\
\mathcal{W}_{2}(\nu_{0},\nu_{1})^{2}+\gamma\tilde{d}_{G}(u,w)^{2}.$

Here the relative distance $\Delta_{p}z$ is determined by the path $p$.
Moreover, $\gamma$ is introduced as a parameter for the thickness of each
layer’s range. We assume that the thickness of each layer is $\sqrt{\gamma}$.
The minimum among all the possible paths is just $\tilde{d}_{G}(u,w)$, the
shortest distance on the graph $G$ between vertices $u$ and $w$. ∎

Using the theorem above, we can compute the minimum $\mathcal{W}_{2}$ cost of
moving a source Gaussian distribution to a Gaussian target distribution.
Indeed, for the $i^{\rm th}$ and $j^{\rm th}$ Gaussian cylinders on
$\mathcal{M}$ we set

$c_{2}(i,j)=\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}+\gamma\tilde{d}_{G}(q_{0}^{i},q_{1}^{j})^{2}$
(21)

If we take $c_{2}(\cdot,\cdot)$ in (21) as the cost matrix and compute the
Kantorovich formulation of OMT, we can derive a distance:

$d_{V_{2}}(\rho_{0},\rho_{1})^{2}=\min_{\pi\in\Pi(p_{0},p_{1})}\sum_{i,j}c_{2}(i,j)\pi(i,j).$
(22)

This distance is derived from the continuous manifold but actually it can be
shown to be a metric for our original vector-valued distributions. See the
proof in Appendix E.

In addition to the latter distance, we can get the optimal transport plan
$\tilde{\pi}_{2}^{*}(\cdot,\cdot)$ from the optimal solution of (22). The
transport plan gives the combination of weights for how the “Gaussian
cylinders” are arranged at different orbits within each layer’s range.

Based on the optimal transport plan, a geodesic (proof in Appendix F) on the
manifold $\mathcal{M}$ may be expressed in the following manner:

$\rho_{t}=\sum_{i,j}\tilde{\pi}_{2}^{*}(i,j)\nu_{t}^{ij}U_{z}(path_{G}(q_{0}^{i},q_{1}^{j},t)),$
(23)

where $U_{z}(z_{0})$ is the 1D uniform distribution density function centered
at $z_{0}$. The distribution
$\nu_{t}^{ij}U_{z}(path_{G}(q_{0}^{i},q_{1}^{j},t))$ at time $t$ is considered
to be the zero expansion on the whole manifold, which is supported on the
chart defined by the shortest path that connects $q_{0}^{i}$ to $q_{1}^{j}$ on
the graph, expressed in its own local continuous coordinates for the $z$ index
along that path.

For each pair of Gaussians, source and target, a deformation Gaussian cylinder
moves across layers following the shortest path on the graph. When it moves
across a layer boundary, the Gaussian cylinder is cut into two parts belonging
to the ranges of two adjacent layers. Each part is still a Gaussian cylinder.
Hence after projection of each layer’s range, the projected distribution is
still a vector GMM distribution. Now if we project (23) onto the range of each
layer, we get the following displacement interpolation in the original space:

$\rho_{t}=\sum_{i,j}\tilde{\pi}_{2}^{*}(i,j)\nu_{t}^{ij}\vec{\delta}_{path_{G}(q_{0}^{i},q_{1}^{j},t)}.$
(24)

Note this form of displacement interpolation is very similar to (17) just with
slightly different weights.

Remark 1: Actually the way that we define vector GMM makes it already a
manifold (each layer is its own chart). However, it is impossible to define
charts that contain all the possible paths, which is to say the atlas contains
only the geometric information within each layer. Our manifold $\mathcal{M}$
on the other hand, has an atlas that contains all the possible paths which
encode global geometric information for which we can solve the OMT problem.

Remark 2: This approach gives a geometric intuitive understanding of the
optimal transport for vector-valued GMMs. Comparing to the cost matrix (15) in
Approach 1, (21) just uses the sum of squares instead of the direct sum, in
analogy to the difference between $\mathcal{W}_{1}$ and $\mathcal{W}_{2}$. It
may seem complicated to consider all the possible paths on the graph in that
approach. But in actual computation, we only need to consider the shortest
path on the $G$ that connects two layers. Moreover, from the parameter
$\gamma$ that also appears in Approach 1, we find a clear geometric meaning:
it represents the thickness of each layer.

##  9 OMT for unbalanced distributions

Unbalanced OMT problem seeks to remove the restriction that the source and
target have the same total mass. The total density of the spaces may be
different, which means that mass can be created or destroyed during the
process of transport. Very commonly, we encounter unbalanced data typically
related to imaging problems in which intensity is taken to be mass. In
medicine, mass may be created with cell proliferation and destroyed with cell
apoptosis. Normalizing both inputs to be probability distributions is not
natural and may cause numerical issues and in some cases loss of information.
Thus we want to find a way to directly employ unnormalized data for the OMT on
Gaussian mixtures. In short, for GMMs, we want to use the general form (1) but
without the restriction (3) i.e., source and target distributions can have
different total densities.

##  10 Vector-valued GMM OMT and unbalanced GMM

The idea comes from [21], where we found that by adding a source layer to the
original scalar problem and using a special set of weight parameters, we can
reformulate the unbalanced OMT problem as a version of the vector-valued OMT
problem. But that work was based on Benamou-Brenier’s computational fluid
dynamics (CFD) form of OMT [3] and a Fisher-Rao source term [8]. Here we will
show that we can use the same source layer trick in the Kantorovich
formulation of OMT problem, so that can be used to solve unbalanced GMM
problem as well.

The idea is to add a source layer to both the source and target distributions
and place the mass difference in that layer, so that the problem becomes
balanced but vector-valued. If the starting distribution has more total mass,
then the mass difference is added to the source layer of target distribution.
On the other hand, if the target distribution has more total mass, then the
mass difference is added to the source layer of starting distribution. So
instead of considering the $n_{0}$ Gaussians to $n_{1}$ Gaussians one-layer
problem, now we consider $n_{0}$ Gaussians to $(n_{1}+1)$ Gaussians (or
$(n_{0}+1)$ Gaussians to $n_{1}$ Gaussians) two-layer vector problem.

We should note that we do not need to specify the distribution (a Gaussian
mixture) in the source layer in the first place, but just to consider a single
implicit Gaussian distribution in that layer. In the $c(\cdot,\cdot)$ matrix
of size $n_{0}\times(n_{1}+1)$ or $(n_{0}+1)\times n_{1}$, the sub-matrix
$c(1:n_{0},1:n_{1})$ may be computed as before (11). We put $\gamma$ in the
last row (or column), which is to say the cost of using source for unit mass
is $\gamma$. Note that $\gamma$ is a parameter to control the weight of the
source in the total distance formula.

After we have solved this extended Kantorovich problem, we can compute a
displacement interpolation, since we have the optimal solution of
$\pi^{*}(\cdot,\cdot)$. Note that the meaning of the last row (or the last
column) is the amount of source the $j^{\rm th}$ Gaussian uses in the target
distribution (or $i^{\rm th}$ Gaussian use in the starting distribution). So
we can add exactly the necessary amount of mass as a weighted same shaped
Gaussian in the source layer for each corresponding Gaussian in the original
space. Then for the two explicit 2-vector GMMs, we can compute the
displacement interpolations between them. The interpolation of only the
original channel is the displacement interpolation for the two unbalanced
GMMs.

Remark Actually we can use this approach more generally than the GMM or source
layer case. Indeed, for any two unbalanced discrete distributions, one can add
an implicit “source node” to the source or target distribution in order to
make the two extended distributions balanced. Then if we consider the extended
Kantorovich problem, we can get an OMT related “distance” between these two
unbalanced distributions.

##  11 Vector GMM OMT and unbalanced vector GMM problem

We can treat the unbalanced vector GMM problem in this framework as well. The
problem of computing the OMT distance and the associated interpolation for
color image data falls into this category, which makes it of great practical
usefulness.

The idea is to add a new source layer for both the source and target vector
GMM structures. The source layer is made to connect with each of the other
layers. If we again put an implicit Gaussian function in that layer, we can
get a balanced but extended version of vector GMM OMT problem, which can be
solved through the scheme from Section 7. The interpolation of the original
part of two extended structures (the part other than the source layers) is
then the displacement interpolation for the two unbalanced vector GMMs.

##  12 Numerical results

###  12.1 1D examples: approach 1 vs approach 2

We first look at the simplest 1D 2-channel vector GMM problem. Figure 6 shows
the distributions of starting and targeting vector distributions.

Figure 6: 1D 2-channel example: starting and target vector distributions

We test the same example via both our approach 1 (see Figure 7) and approach 2
(see Figure 8).

i͡n 1,3,5,7,9

Figure 7: 1D 2-channel Example 1 - using approach 1: vector distributions over
time

i͡n 1,3,5,7,9

Figure 8: 1D 2-channel example 2 - using approach 2: vector distributions over
time

Two sets of results look very similar. For most examples that we have tested,
we cannot find discernible difference between these two approaches. This is
due to the similar and correlated formulations of two cost matrices. For later
examples, we will use approach 2 formulation.

###  12.2 1D examples: graph structure

We consider the simplest 1D 3-channel problem, which is the simplest example
with nontrivial graph structure. Figure 9 shows the distributions of the
starting and target vector distributions. The distributions in their 3
channels are explicitly and separately plotted. All the mass of the starting
vector distribution locates in its first channel (one Gaussian), and all the
mass of the target vector distribution is two Gaussians which locate in the
third channel.

We tested two different graph structures for the vector distribution problem.
In Figure 10, we show the interpolation results whose second channel is
connected to first and third ones, but whose first and third channels are not
directly connected. In Figure 11, we show the interpolation results which
three layers are fully connected.

We can see from the results that the mass transfer in 10 goes through the
second layer while in 11 mass in layer one can directly transfer to layer
three.

Figure 9: 1D 3 channels example: starting and target vector distributions

i͡n 1,3,5,7,9

Figure 10: 1D 3-channel Example 1 - not fully connected graph: vector
distributions over time

i͡n 1,3,5,7,9

Figure 11: 1D 3-channel example 2 - fully connected graph: vector
distributions over time

###  12.3 2D examples: parameter $\gamma$

Now we consider 2D 3-vector distributions, and we take a graph which is not
completely connected: the green channel is connected to red channel and blue
channel, but red and blue channels are not directly connected. This is exactly
the model for RGB three-channel images. So the following vector-valued
distributions are plotted as color images. We still use approach 2 method for
this numerical test.

There are two different kinds of mass transfer: within layers and between
layers. Parameter $\gamma$ controls the amount of mass that moves across
layers because with $\gamma$ large, the cost of moving mass between layers is
large. By using different values of $\gamma$, we may get very different
results.

We use a two-Gaussian to two-Gaussian example to show the effect of different
values of $\gamma$. Both the starting and target vector-valued distributions
have one Gaussian in the red layer and one in the blue layer. But the location
of the balls are interchanged (see Figure 12). With $\gamma$ large, the two
Gaussians just move within its own layer (see Figure 13). With $\gamma$ small,
red ball goes through green layer to blue layer while blue ball takes the
inverse path (see Figure 14).

Figure 12: 2D 3 channels example: starting and target vector distributions

i͡n 1,3,5,7,9

Figure 13: 2D Example 1 - $\gamma$ large: vector distributions over time

i͡n 1,3,5,7,9

Figure 14: 2D Example 2 - $\gamma$ small: vector distributions over time

###  12.4 Fitted image data

The above examples use GMMs as inputs. For data that is represented as a GMM,
we can first fit the data to GMM form using method like EM algorithm [15] and
use the fitted GMMs as inputs.

Figure 15 is the image data we use. We fit each layer with 20 Gaussians. The
fitting results are shown in Figure 16. The interpolations between them are
computed in Figure 17.

Figure 15: image data example

Figure 16: Fitted GMMs as starting and target

i͡n 1,3,5,7,9

Figure 17: Image Example: vector distributions over time

###  12.5 Unbalanced GMM OMT

We also tested on unbalanced GMM data. The total mass of two Gaussians on the
right in the starting distribution is equal to the mass of the large Gaussian
at the center of the target distribution. Two small Gaussians on the left of
the target GMM are the mass difference. See Figure 18.

As expected, the two Gaussians in the starting distribution move together and
combine as the large Gaussian at the center of target distribution while the
presence of two small Gaussians on the left is due to the injection of source
(see Figure 19). We also plotted the distribution function of both original
layer and source layer as two surfaces in 3D (See Figure 20). It is clear from
the surfaces plots that there are two Gaussian distributions in the source
layer as the source of the mass difference in the original layer.

Figure 18: Fitted GMMs as starting and target

i͡n 1,3,5,7,9

Figure 19: Unbalanced Example: vector distributions over time

i͡n 1,3,5,7,9

Figure 20: Unbalanced Example: vector distributions and source layer plotted
as surface

##  13 Conclusion

This work focuses on the optimal transport for vector GMM which is a special
case of vector OMT. As an extension of [7], we define a distance and geodesic
path for vector GMM distributions. The proposed transport is of great value
because of its speed advantage and unique ability of preserving structure. We
also consider the unbalanced GMM problems and extend our source layer idea in
[21] to the Kantorovich formulation. We are going to use this framework to
analyze medical images and other appropriate vector-valued distribution data.
And potentially we can extract more information by using this.

##  14 Acknowledgments

This study was supported by AFOSR grants (FA9550-17-1-0435, FA9550-20-1-0029),
NIH grant (R01-AG048769), MSK Cancer Center Support Grant/Core Grant (P30
CA008748), and a grant from Breast Cancer Research Foundation (grant
BCRF-17-193).

##  Appendix A Approach 0: restrictions on transport matrix

Specifically, we add the following constraints on $\Pi(p_{0},p_{1})$ based on
the graph structure:

$\Pi_{G}(p_{0},p_{1})=\left\\{\pi\in{\mathbb{R}}_{+}^{n_{0}\times
n_{1}}|\pi\vec{1}_{n_{1}}=p_{0},\ \pi^{T}\vec{1}_{n_{0}}=p_{1},\pi(i,j)=0\
\text{if}\ (q^{i}_{0},q^{j}_{1})\notin E(G)\cup\\{(v,v)|v\in
V(G)\\}\right\\}.$ (25)

Thus transport is allowed between two Gaussians only when they are in the same
channel or they located in adjacent channels.

Accordingly, we can define a “distance:”

$d_{V_{0}}(\rho_{0},\rho_{1})^{2}=\inf_{\pi\in\Pi_{G}(p_{0},p_{1})}\sum_{i,j}\pi(i,j)\mathcal{W}_{2}(\nu^{i}_{0},\nu^{j}_{1})^{2}.$
(26)

This “distance” is a pseudo-metric, which is to say, $d_{V_{0}}(\cdot,\cdot)$
satisfies all the other metric conditions, but may be zero for two distinct
distributions. Indeed, $d_{V_{0}}(\cdot,\cdot)$ is zero if one vector GMM is a
permutation of channels of the other GMM.

######  Theorem 4.

$d_{V_{0}}(\cdot,\cdot)\geq 0$ and satisfies triangle inequality.

###### Proof.

Clearly, $d_{V_{0}}(\rho_{0},\rho_{1})\geq 0$ for any
$\rho_{0},\rho_{1}\in\mathcal{G}({\mathbb{R}}^{N})$. To prove the triangular
inequality $d_{G}(\rho_{0},\rho_{1})+d_{G}(\rho_{1},\rho_{2})\geq
d_{G}(\rho_{0},\rho_{2})$, let $\pi_{01}$($\pi_{12}$) is the optimal solution
with marginal $\rho_{0},\rho_{1}$($\rho_{1},\rho_{2}$).

Denote

$\pi_{02}(i,k)=\sum_{j=1}^{n_{1}}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}.$

Clearly $\pi_{02}$ is a joint distribution between $\rho_{0}$ and $\rho_{2}$,
namely, $\pi_{02}\in\Pi_{G}(\rho_{0},\rho_{2})$. It follows from direct
calculation:

$\displaystyle\sum_{i}\pi_{02}(i,k)$
$\displaystyle=\sum_{i,j}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}$
$\displaystyle=\sum_{j}\frac{p_{1}^{j}\pi_{12}(j,k)}{p_{1}^{j}}$
$\displaystyle=p_{2}^{k}.$

Similarly we have $\sum_{k}\pi_{02}(i,k)=p_{0}^{i}$. Therefore,

$\displaystyle d_{V_{0}}(\rho_{0},\rho_{2})$
$\displaystyle\leq\sqrt{\sum_{i,k}\pi_{02}(i,k)\mathcal{W}_{2}(\nu_{0}^{i},\nu_{2}^{k})^{2}}$
$\displaystyle\leq\sqrt{\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}(\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})+\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k}))^{2}}$
$\displaystyle\leq\sqrt{\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}}+\sqrt{\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k})^{2}}$
$\displaystyle=\sqrt{\sum_{i,j}\pi_{01}(i,j)\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}}+\sqrt{\sum_{j,k}\pi_{12}(j,k)\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k})^{2}}$
$\displaystyle=d_{V_{0}}(\rho_{0},\rho_{1})+d_{V_{0}}(\rho_{1},\rho_{2}).$

∎

A notion of _geodesic_ may be defined as follows:

$\rho_{t}=\sum_{i,j}\pi^{*}(i,j)\nu_{t}^{ij}\left((1-t)\vec{\delta}_{q_{0}^{i}}+t\vec{\delta}_{q_{1}^{j}}\right),$
(27)

where $\pi^{*}(i,j)$ is the optimal solution of (26) and $\nu_{t}^{ij}$ is the
displacement interpolation in (8).

######  Theorem 5.

$d_{V_{0}}(\rho_{s},\rho_{t})=(t-s)d_{V_{0}}(\rho_{0},\rho_{1}),\qquad 0\leq
s<t\leq 1$

###### Proof.

$\displaystyle d_{V_{0}}(\rho_{s},\rho_{t})$
$\displaystyle\leq\sqrt{\sum_{i,j}\pi^{*}(i,j)\mathcal{W}_{2}(\nu_{s}^{i},\nu_{t}^{j})^{2}}$
$\displaystyle=(t-s)\sqrt{\sum_{i,j}\pi^{*}(i,j)\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}}$
$\displaystyle=(t-s)d_{V_{0}}(\rho_{0},\rho_{1})$

By the triangle inequality,

$\displaystyle d_{V_{0}}(\rho_{0},\rho_{1})$ $\displaystyle\leq
d_{V_{0}}(\rho_{0},\rho_{s})+d_{V_{0}}(\rho_{s},\rho_{t})+d_{V_{0}}(\rho_{t},\rho_{1})$
$\displaystyle\leq
sd_{V_{0}}(\rho_{0},\rho_{1})+(t-s)d_{V_{0}}(\rho_{0},\rho_{1})+(1-t)d_{V_{0}}(\rho_{0},\rho_{1})$
$\displaystyle=d_{V_{0}}(\rho_{0},\rho_{1}).$

These two inequalities give the result
$d_{V_{0}}(\mu_{s},\mu_{t})=(t-s)d_{V_{0}}(\mu_{0},\mu_{1})$. ∎

###  A.1 Limitations of approach 0

There are several problems with this approach. So that we propose two other
approaches in the main text trying to ameliorate in what follows:

  * •

There may be no solution, because $\Pi_{G}(p_{0},p_{1})$ may be an empty set.
For example, in a 3-channel structure, with channel 2 connected to channels 1
and 3, but 1 and 3 are not directly connected. And source and target are
described as follows:
Source: all the mass locates in channel 1, $p_{0}=(\nu_{0},0,0)^{T}$.
Target: all the mass locates in channel 3, $p_{1}=(0,0,\nu_{1})^{T}$.
From the definition of joint distribution, we have
$\sum_{i}\pi(i,j)=p_{1}^{j}$, $\sum_{j}\pi(i,j)=p_{0}^{i}$. From the zeros in
$p_{0}$ and $p_{1}$, we can see $\pi(i,j)=0$ when $i=2,3$ or $j=1,2$.
Meanwhile, $\pi(1,3)=0$ because of the restriction, $\pi(i,j)=0$ if channel
$i$ and channel $j$ are not directly connected. In a word, all $\pi(i,j)$ is
$0$ here, which contradicts to the definition of joint distribution. So such
joint distribution does not exist.

  * •

The transport can only go through two connected channels.
This is a problem because even two channels are not connected but there is a
path between them, the transport can go along the path.

  * •

The cost of mass transport between channels does not include.
In this approach, a Gaussian can go to the different channels freely. This
does not make sense because it is a nontrivial mass transport.

  * •

Not a metric.
Like we discussed above, the distance can be $0$ when two marginals are
different.

Despite these latter limitations, in practice this approach is very useful
when employed for vector-valued GMM with a fully connected graph. This
approach gives reasonable results for color images as we will show below.

###  A.2 Numerical results of approach 0

We consider 2D 3-vector distributions with fully connected graph structure.
Like what we did in Section 12.3, we use a two-Gaussian to two-Gaussian
example. Both starting and target vector distribution has one Gaussian in the
red layer and one in the blue layer, but the location are interchanged (see
Figure 21). As we mentioned above, the Gaussian can go between channel freely,
so the blue ball just go to the red layer and the red ball go to the blue
layer (see Figure 22).

Figure 21: 2D 3 channels example of approach 0: starting and target vector
distributions

i͡n 1,3,5,7,9

Figure 22: 2D Example of approach 0: vector distributions over time

##  Appendix B Proof of Theorem 1

This section focus on a new distance and its properties. The distance is
defined as:

$d_{V_{1}}(\rho_{0},\rho_{1})=\inf_{\pi\in\Pi(p_{0},p_{1})}\sum_{i,j}\pi(i,j)[\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})+\gamma\tilde{d}_{G}(q_{0}(i),q_{1}(j))]$

######  Theorem 6.

$d_{V_{1}}(\cdot,\cdot)$ defines a metric on
$\mathcal{G}({\mathbb{R}}^{N}\times G)$.

###### Proof.

Clearly, $d_{V_{1}}(\rho_{0},\rho_{1})\geq 0$,
$\forall\rho_{0},\rho_{1}\in\mathcal{G}(\mathbb{R}^{N})$ and
$d_{V_{1}}(\rho_{0},\rho_{1})=0$ if and only if $\rho_{0}=\rho_{1}$.

We now prove the triangle inequality $d_{V_{1}}(\rho_{0},\rho_{2})\leq
d_{V_{1}}(\rho_{0},\rho_{1})+d_{V_{1}}(\rho_{1},\rho_{2})$. Similar to our
argument in Theorem 4, denote

$\pi_{02}(i,k)=\sum_{j=1}^{n_{1}}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}.$

As $\pi_{02}$ is a joint distribution between $\rho_{0}$ and $\rho_{2}$, we
have

$\displaystyle d_{V_{1}}(\rho_{0},\rho_{2})$
$\displaystyle\leq\sum_{i,k}\pi_{02}(i,k)[\mathcal{W}_{2}(\nu_{0}^{i},\nu_{2}^{k})+\gamma\tilde{d}_{G}(q_{0}(i),q_{2}(k))]$
$\displaystyle\leq\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}[\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})+\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k})+\gamma\tilde{d}_{G}(q_{0}(i),q_{1}(j))+\gamma\tilde{d}_{G}(q_{1}(j),q_{2}(k))]$
$\displaystyle=\sum_{i,j}\frac{\pi_{01}(i,j)p^{j}_{1}}{p^{j}_{1}}[\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})+\gamma\tilde{d}_{G}(q_{0}(i),q_{1}(j))]$
$\displaystyle+\sum_{j,k}\frac{p^{j}_{1}\pi_{12}(j,k)}{p^{j}_{1}}[\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k})+\gamma\tilde{d}_{G}(q_{1}(j),q_{2}(k))]$
$\displaystyle=d_{V_{1}}(\rho_{0},\rho_{1})+d_{V_{1}}(\rho_{1},\rho_{2})$

∎

##  Appendix C Proof of theorem 2

###### Proof.

$\rho_{t}=\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)\nu_{t}^{ij}\vec{\delta}_{path_{G}(q_{0}^{i},q_{1}^{j},t)}$

Locally, let $\bar{d}=\tilde{d}_{G}(q_{0}(i),q_{i}(j))$. Rename the vertices
on this shortest path in order $0-1-\cdots-\bar{d}$. Denote
$a=t\bar{d}-\lfloor t\bar{d}\rfloor$, $b=s\bar{d}-\lfloor s\bar{d}\rfloor$,
$r=\lfloor t\bar{d}\rfloor-\lfloor s\bar{d}\rfloor$. Noticing that, all above
notations depend on $i,j$. Using these notations, the displace interpolation
can be expressed as:

$\rho_{t}=\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)\nu_{t}^{ij}((1-a)\vec{\delta}_{\lfloor
t\bar{d}\rfloor}+a\vec{\delta}_{\lfloor t\bar{d}\rfloor+1})$

$\rho_{s}=\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)\nu_{t}^{ij}((1-b)\vec{\delta}_{\lfloor
s\bar{d}\rfloor}+b\vec{\delta}_{\lfloor s\bar{d}\rfloor+1})$

Then we have

$\displaystyle d_{V_{1}}(\rho_{s},\rho_{t})$
$\displaystyle\leq\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)[W_{2}(\nu_{s}^{ij},\nu_{t}^{ij})+\gamma(r(1-a)(1-b)+rab+(r-1)(1-a)b+(r+1)(1-b)a)]$
$\displaystyle=\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)[W_{2}(\nu_{s}^{ij},\nu_{t}^{ij})+\gamma(r+a-b)]$
$\displaystyle=\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)[(t-s)W_{2}(\nu_{0}^{i},\nu_{1}^{j})+\gamma(t-s)\bar{d}$
$\displaystyle=(t-s)\sum_{i,j}\tilde{\pi}_{1}^{*}(i,j)[W_{2}(\nu_{0}^{i},\nu_{1}^{j})+\gamma\tilde{d}_{G}(q_{0}(i),q_{i}(j))]$
$\displaystyle=(t-s)d_{V_{1}}(\rho_{0},\rho_{1})$

By the triangle inequality,

$\displaystyle d_{V_{1}}(\rho_{0},\rho_{1})$ $\displaystyle\leq
d_{V_{1}}(\rho_{0},\rho_{s})+d_{V_{1}}(\rho_{s},\rho_{t})+d_{V_{1}}(\rho_{t},\rho_{1})$
$\displaystyle\leq
sd_{V_{1}}(\rho_{0},\rho_{1})+(t-s)d_{V_{1}}(\rho_{0},\rho_{1})+(1-t)d_{V_{1}}(\rho_{0},\rho_{1})$
$\displaystyle=d_{V_{1}}(\rho_{0},\rho_{1}).$

These two inequalities give the result
$d_{V_{1}}(\mu_{s},\mu_{t})=(t-s)d_{V_{1}}(\mu_{0},\mu_{1})$. ∎

##  Appendix D Background on manifolds

The books [13] and [18] introduce the basic concepts of manifolds. Generally,
manifold may be described as follows:

######  Definition 1.

An n-dimensional manifold is a topological space with the property that each
point has a neighborhood that is homeomorphic to $\mathbb{R}^{n}$.

######  Definition 2.

A homeomorphism is a continuous function between topological spaces that has a
continuous inverse function. Two spaces with a homeomorphism between them are
called homeomorphic.

The simplest n-dimensional manifold is $\mathbb{R}^{n}$. There are two
important concepts here, charts and atlas.

######  Definition 3.

A chart is an invertible map between a subset of the manifold and a subset of
$\mathbb{R}^{n}$. A chart can be denoted by $(U,\phi:U\to\mathbb{R}^{n})$.

######  Definition 4.

Two charts $(U,\phi:U\to\mathbb{R}^{n}),(V,\psi:V\to\mathbb{R}^{n})$ of a
manifold are $C^{\infty}$-compatible if the two maps

$\phi\circ\psi^{-1}:\psi(U\cap V)\to\phi(U\cap
V)\qquad\psi\circ\phi^{-1}:\phi(U\cap V)\to\psi(U\cap V)$

are $C^{\infty}$. These two maps are called the transition functions between
the charts.

######  Definition 5.

An atlas is an specific collection
$\mathcal{U}=\\{(U_{\alpha},\phi_{\alpha})\\}$ of pairwise
$C^{\infty}$-compatible charts which covers a manifold $\mathcal{M}$, i.e.,
such that $\mathcal{M}=\bigcup_{\alpha}U_{\alpha}$.

##  Appendix E metric $d_{V_{2}}(\cdot,\cdot)$

######  Theorem 7.

$d_{V_{2}}(\cdot,\cdot)$ defines a metric on
$\mathcal{G}({\mathbb{R}}^{N}\times G)$

Clearly, $d_{V_{2}}(\rho_{0},\rho_{1})\geq 0$ for any
$\rho_{0},\rho_{1}\in\mathcal{G}({\mathbb{R}}^{N}\times G)$ and
$d_{V_{2}}(\rho_{0},\rho_{1})=0$ if and only if $\rho_{0}=\rho_{1}$. We next
prove the triangular inequality, namely,

$d_{V_{2}}(\rho_{0},\rho_{1})+d_{V_{2}}(\rho_{1},\rho_{2})\geq
d_{V_{2}}(\rho_{0},\rho_{2})$

Before proving the triangle inequality, we first need the following lemma:

######  Lemma 1.

$\sqrt{\sum_{i=1}^{n}(a_{i}+b_{i})^{2}+(c_{i}+d_{i})^{2}}\leq\sqrt{\sum_{i=1}^{n}a_{i}^{2}+c_{i}^{2}}+\sqrt{\sum_{i=1}^{n}b_{i}^{2}+d_{i}^{2}},$
(28)

$\forall\ a_{i},b_{i},c_{i},d_{i}\in{\mathbb{R}},\ i=1,...,n$ and $n\geq 1$.

###### Proof.

$\displaystyle\mbox{LHS}^{2}$
$\displaystyle=\sum(a_{i}^{2}+b_{i}^{2}+c_{i}^{2}+d_{i}^{2}+2a_{i}b_{i}+2c_{i}d_{i})$
$\displaystyle\leq\sum(a_{i}^{2}+b_{i}^{2}+c_{i}^{2}+d_{i}^{2})+2(\sqrt{\sum
a_{i}^{2}\sum b_{i}^{2}}+\sqrt{\sum c_{i}^{2}\sum d_{i}^{2}}),$
$\displaystyle(\sqrt{\sum a_{i}^{2}\sum b_{i}^{2}}+\sqrt{\sum c_{i}^{2}\sum
d_{i}^{2}})^{2}$ $\displaystyle=\sum a_{i}^{2}\sum b_{i}^{2}+\sum
c_{i}^{2}\sum d_{i}^{2}+2\sqrt{\sum a_{i}^{2}\sum b_{i}^{2}\sum c_{i}^{2}\sum
d_{i}^{2}}$ $\displaystyle\leq\sum a_{i}^{2}\sum b_{i}^{2}+\sum c_{i}^{2}\sum
d_{i}^{2}+\sum a_{i}^{2}\sum d_{i}^{2}+\sum b_{i}^{2}\sum c_{i}^{2}$
$\displaystyle=(\sum a_{i}^{2}+\sum c_{i}^{2})(\sum b_{i}^{2}+\sum
d_{i}^{2}).$

By combining the above two inequalities, we have

$\displaystyle\mbox{LHS}^{2}$
$\displaystyle\leq\sum(a_{i}^{2}+b_{i}^{2}+c_{i}^{2}+d_{i}^{2})+2\sqrt{(\sum
a_{i}^{2}+\sum c_{i}^{2})(\sum b_{i}^{2}+\sum d_{i}^{2})}$
$\displaystyle=(\sqrt{\sum a_{i}^{2}+c_{i}^{2}}+\sqrt{\sum
b_{i}^{2}+d_{i}^{2}})^{2}$ $\displaystyle=\mbox{RHS}^{2}.$

∎

Next we prove the triangle inequality. Similar to the proof in Appendix 4, we
define

$\pi_{02}=\sum_{j}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p_{1}^{j}}.$

Because $\pi_{02}\in\Pi(p_{0},p_{2})$, we have

$\displaystyle
d_{V_{2}}(\rho_{0},\rho_{2})\leq\sqrt{\sum_{i,k}\pi_{02}(i,k)[W_{2}(\nu_{0}^{i},\nu_{2}^{k})^{2}+\gamma\tilde{d}_{G}(q_{0}^{i},q_{2}^{k})^{2}}]$
$\displaystyle\leq\sqrt{\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}[(\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})+\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k}))^{2}+\gamma(\tilde{d}_{G}(q_{0}^{i},q_{1}^{j})+\tilde{d}_{G}(q_{1}^{j},q_{2}^{k}))^{2}]}$
$\displaystyle\leq\sqrt{\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}[\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}+\gamma\tilde{d}_{G}(q_{0}^{i},q_{1}^{j})^{2}]}+\sqrt{\sum_{i,j,k}\frac{\pi_{01}(i,j)\pi_{12}(j,k)}{p^{j}_{1}}[\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k})^{2}+\gamma\tilde{d}_{G}(q_{1}^{j},q_{2}^{k})^{2}]}$
$\displaystyle=\sqrt{\sum_{i,j}\pi_{01}(i,j)[\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}+\gamma\tilde{d}_{G}(q_{0}^{i},q_{1}^{j})^{2}]}+\sqrt{\sum_{j,k}\pi_{12}(j,k)[\mathcal{W}_{2}(\nu_{1}^{j},\nu_{2}^{k})^{2}+\gamma\tilde{d}_{G}(q_{1}^{j},q_{2}^{k})^{2}]}$
$\displaystyle=d_{V_{2}}(\rho_{0},\rho_{1})+d_{V_{2}}(\rho_{1},\rho_{2}).$

The last inequality used the lemma above.

$\hfill\square$

##  Appendix F Geodesic of $d_{V_{2}}(\cdot,\cdot)$

######  Theorem 8.

$d_{V_{2}}(\rho_{s},\rho_{t})=(t-s)d_{V_{2}}(\rho_{0},\rho_{1}),\qquad 0\leq
s<t\leq 1$

###### Proof.

$\displaystyle d_{V_{2}}(\rho_{s},\rho_{t})$
$\displaystyle\leq\sqrt{\sum_{i,j}\pi^{*}(i,j)(\mathcal{W}_{2}(\nu_{s}^{i},\nu_{t}^{j})^{2}+\gamma\tilde{d}_{G}(q_{s}^{i},q_{t}^{j})^{2})}$
$\displaystyle=(t-s)\sqrt{\sum_{i,j}\pi^{*}(i,j)\mathcal{W}_{2}(\nu_{0}^{i},\nu_{1}^{j})^{2}+\gamma\tilde{d}_{G}(q_{0}^{i},q_{1}^{j})^{2}}$
$\displaystyle=(t-s)d_{V_{2}}(\rho_{0},\rho_{1})$

By the triangle inequality,

$\displaystyle d_{V_{2}}(\rho_{0},\rho_{1})$ $\displaystyle\leq
d_{V_{2}}(\rho_{0},\rho_{s})+d_{V_{2}}(\rho_{s},\rho_{t})+d_{V_{2}}(\rho_{t},\rho_{1})$
$\displaystyle\leq
sd_{V_{2}}(\rho_{0},\rho_{1})+(t-s)d_{V_{2}}(\rho_{0},\rho_{1})+(1-t)d_{V_{2}}(\rho_{0},\rho_{1})$
$\displaystyle=d_{V_{2}}(\rho_{0},\rho_{1}).$

This two inequalities give the result
$d_{V_{2}}(\mu_{s},\mu_{t})=(t-s)d_{V_{2}}(\mu_{0},\mu_{1})$. ∎

## References

  * [1] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein GAN.  arxiv.org, 1701.07875, 2017.
  * [2] Richard Bellman.  On a routing problem.  Quarterly of Applied Mathematics, 16(1):87–90, 1958.
  * [3] Jean-David Benamou and Yann Brenier.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  Numerische Mathematik, 84(3):375–393, 2000.
  * [4] Yongxin Chen, Tryphon T Georgiou, and Allen Tannenbaum.  Interpolation of density matrices and matrix-valued measures: The unbalanced case.  Euro. Jnl of Applied Mathematics, 30(3):458–480, 2018.
  * [5] Yongxin Chen, Tryphon T Georgiou, and Allen Tannenbaum.  Matrix optimal mass transport: a quantum mechanical approach.  IEEE Trans. Automatic Control, 63(8):2612 – 2619, 2018.
  * [6] Yongxin Chen, Tryphon T Georgiou, and Allen Tannenbaum.  Vector-valued optimal mass transport.  SIAM Journal Applied Mathematics, 78(3):1682–1696, 2018.
  * [7] Yongxin Chen, Tryphon T. Georgiou, and Allen Tannenbaum.  Optimal Transport for Gaussian Mixture Models.  IEEE Access, 7(c):6269–6278, 2019.
  * [8] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Unbalanced Optimal Transport: Geometry and Kantorovich Formulation.  working paper or preprint, August 2015.
  * [9] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and Francois-Xavier Vialard.  An interpolating distance between optimal transport and Fisher-Rao metrics.  Foundations of Computational Mathematics, 10:1–44, 2016.
  * [10] E. W. Dijkstra.  A Note on Two Problems in Connexion with Graphs.  Numerische Mathematlk, 1:269–271, 1959.
  * [11] Lester Ford.  Network Flow Theory.  Rand Corporation, 1956.
  * [12] Steven Haker, Lei Zhu, Allen Tannenbaum, and Sigurd Angenent.  Optimal mass transport for registration and warping.  International Journal of Computer Vision, 60(3):225–240, 2004.
  * [13] Jürgen Jost and Jèurgen Jost.  Riemannian geometry and geometric analysis, volume 42005.  Springer, 2008.
  * [14] James C. Mathews, Saad Nadeem, Maryam Pouryahya, Zehor Belkhatir, Joseph O. Deasy, Arnold J. Levine, and Allen R. Tannenbaum.  Functional network analysis reveals an immune tolerance mechanism in cancer.  Proceedings of the National Academy of Sciences, 117(28):16339–16345, jul 2020.
  * [15] Geoffrey McLachlan and David Peel.  Finite Mixture Models.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, USA, sep 2000.
  * [16] Karl Pearson.  III. Contributions to the mathematical theory of evolution.  Philosophical Transactions of the Royal Society of London. (A.), 185:71–110, dec 1894.
  * [17] Svetlozar T Rachev and Ludger Rüschendorf.  Mass Transportation Problems: Volumes I and II.  Springer Science & Business Media, 1998.
  * [18] Loring W.. Tu.  An introduction to manifolds.  Springer., 2011.
  * [19] Cédric Villani.  Topics in Optimal Transportation.  American Mathematical Soc., 2003.
  * [20] Cédric Villani.  Optimal Transport: Old and New, volume 338.  Springer Science & Business Media, 2008.
  * [21] J. Zhu, R. Elkin, J. H. Oh, J. O. Deasy, and A. Tannenbaum.  A vectorial approach to unbalanced optimal mass transport.  IEEE Access, 8:209224–209231, 2020.