\section{Introduction}
\label{sec:intro}
For a metric space $(\mathcal{X},d)$ and a weighting exponent
$p \in [1,\infty)$ optimal transport induces the $p$-Wasserstein
distances $W_p$ on the probability measures over $\mathcal{X}$.  A
remarkable property of Wasserstein distances is that they form a
length space if the base space $(\mathcal{X},d)$ is a length space,
inducing the so-called displacement interpolation between probability
measures \cite{McCannConvexity1997}. The celebrated Benamou--Brenier
formula for $W_2$ over $\mathbb{R}^n$ \cite{BeBr00} can be interpreted as an
explicit search for the shortest path between two probability
measures. In the last two decades the geometry of metric spaces has
extensively been studied by means of optimal transport. In explicit it
has been observed that the 2-Wasserstein metric over probability
densities in $\mathbb{R}^n$ formally resembles a Riemannian manifold
\cite{Ot01} and that various diffusion-type equations can be
interpreted as gradient flows for entropy-type functionals with
respect to this metric \cite{JKO1998}. For a comprehensive
introduction we refer to the monographs
\cite{Villani-OptimalTransport-09,Santambrogio-OTAM}.

Unfortunately, this rich geometry is not directly available when the
base space $\mathcal{X}$ is discrete, since $W_2$ degenerates and does
not admit geodesics.
Maas \cite{Ma11} introduced a transport-type Riemannian metric $\mathcal{W}$
on probability measures over a discrete space $\mathcal{X}$ equipped
with a reversible Markov kernel $Q$, based on an adaption of the
Benamou--Brenier formula.  A key ingredient in the construction is the
choice of a \emph{`mass averaging'} function $\theta$ that
interpolates the amount of mass on neighbouring graph vertices. For
the particular choice of $\theta$ being the logarithmic mean, the heat
equation (with respect to the underlying Markov kernel) arises as
gradient flow of the entropy with respect to this metric \cite{Ma11,
  Mielke-Nonlinearity-2011}, yielding a discrete analogue of Otto's
interpretation of diffusive PDEs, see also \cite{MaasDiscretePME2014}
for a generalization to non-linear evolution equations on discrete
spaces. In analogy to the Lott--Sturm--Villani theory the displacement
interpolation on graphs has been used to introduce a notion of Ricci
curvature lower bounds for discrete spaces equipped with Markov
kernels \cite{ErbarMaas-Ricci2012} that implies a variety of
functional inequalities in analogy to the theory of
Lott--Sturm--Villani. The study of transport-type distances on
discrete domains has various connections to the original Wasserstein
distances on continuous domains.  Approximating a torus with an
increasingly finer toroidal graph, the discrete transport metric $\mathcal{W}$
has been shown to converge to the continuous underlying 2-Wasserstein
distance on the torus in the sense of Gromov--Hausdorff \cite{GiMa13}.
Conversely, the introduction of a mass averaging function for discrete
spaces has in turn inspired the design of new non-local transport-type
metrics in continuous domains \cite{Erbar-Jump2014}.
\medskip

Computing classical Wasserstein distances $W_2$ numerically is often a
challenge.  While the classical Kantorovich formulation via transport
couplings is a standard linear program, its na\"ive dense form
requires $(\mathrm{card} \mathcal{X})^2$ variables which may quickly become
computationally unfeasible as $\mathcal{X}$ increases in size.  On
arbitrary metric graphs $(\mathcal{X},d)$ an additional problem arises:
only local edge lengths are usually prescribed and the full distance
function $d : \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is in general unknown
a priori. On large graphs, computing $d$ from local edge lengths may
be computationally prohibitive or even storing $d$ may exceed the
memory capacities.

Owing to its particular structure, the 1-Wasserstein distance over a
discrete graph can be reformulated as a min cost flow problem along
its edges, thus drastically reducing the number of required variables
if the graph is sparse and requiring no pre-computation of $d$, see
for instance \cite{NetworkFlows1993}. On continuous domains this
corresponds to \emph{Beckmann's problem} \cite{Santambrogio-OTAM}.  A
numerical scheme tailored to application on meshed surfaces is
presented in \cite{SolomonEMDSurfaces2014}.  A computational approach
that uses quadratic regularization to break the non-uniqueness of the
optimal flow is described in \cite{EssidSolomon-QuadraticRegOT-2017}.

For the 2-Wasserstein distance on continuous domains the
Benamou--Brenier formula serves a similar purpose, see for instance
\cite{PaPe14} for a numerical scheme based on proximal point
algorithms. However, this does not immediately carry over to discrete
graphs, as the mass averaging function $\theta$ introduces a
non-trivial coupling of the mass variables along graph edges.  In
\cite{SolomonTransportGraphs2016} a Benamou--Brenier-type transport
distance on discrete metric graphs is developed, similar to the
construction of Maas, and a corresponding numerical scheme is
developed. A crucial design choice is that $\theta$ is picked to be
the harmonic mean which allows the application of second-order convex
cone programs for numerical optimization.  This does not extend to
other choices of $\theta$ and thus, for instance, hinders the
numerical study of the gradient flow when $\theta$ is the logarithmic
mean.

\subsubsection*{Contribution}  In this article we present a scheme
for the numerical approximation of the distance $\mathcal{W}$ on discrete sets
$\mathcal{X}$ equipped with irreducible Markov kernels $Q$ as introduced
by Maas.  We pick up the Benamou--Brenier-type formulation and provide
a temporal discretization of the action functional to obtain a
finite-dimensional convex problem and prove $\Gamma$-convergence of
the discretized functional to the original problem, as well as strong
convergence of the discrete geodesics to the continuous geodesics.  To
overcome the strong coupling of mass variables along graph edges
caused by the mass averaging function we introduce a set of slack
variables to remedy this entanglement.  This allows us to apply a
robust proximal point algorithm for the optimization. Due to the slack
variables, all involved proximal mappings can be computed efficiently
by either solving a sparse linear program (if $Q$ is sparse) or by
decomposing them into independent low-dimensional sub-problems.

In particular this numerical scheme does not depend critically on the
choice of $\theta$ and can be quickly adapted to different
variants. We provide formulas for the logarithmic and geometric mean.
For a series of numerical test cases we visualize and discuss the
behaviour of the interpolating flow.
Finally, we adopt the algorithm to approximate gradient flows with
respect to the discrete transportation distance $\mathcal{W}$. In particular,
we test the algorithm against a classical backward Euler
discretization of the heat equation on a graph which coincides with
the gradient flow of the entropy.

\subsubsection*{Organization} The paper is organized as follows. At first we review the
construction of the $L^2$-Wasserstein metric on discrete spaces by
Maas \cite{Ma11} in Section \ref{sec:review}.  Then, in Section
\ref{sec:Discretization} we will derive the time discretization and
establish $\Gamma$-convergence of the time discrete action functional
and the convergence of time discrete geodesics to a continuous
geodesic.  Next, the proximal splitting algorithm with suitably chosen
slack variables is presented in detail in Section
\ref{sec:ProximalSplitting}.  Numerical results are discussed in
Section \ref{sec:numerics} and the experimental comparison of
solutions of a JKO scheme for the entropy and solutions of the Markov
semigroup are presented in Section \ref{sec:gradflow}.

\section{Optimal transport on graphs}
\label{sec:review}
In this section we briefly review the discrete transportation metric
on the space of probability measures over a graph and in particular
recall the basic definitions and discuss the analogy to the
$L^2$-Wasserstein metric on probability measures over $\mathbb{R}^n$. Then we
derive a priori bounds on feasible curves of measures.

\subsection{The discrete transportation distance}
Let $\mathcal{X}$ be a finite set and let
$Q : \mathcal{X} \times \mathcal{X} \to [0,\infty)$ be the transition rate
matrix of a continuous time Markov chain on $\mathcal{X}$.
I.e.~we have $Q(x,y)\geq0$ for
$x\neq y$ and make the convention that $Q(x,x)=0$ for all
$x\in\mathcal{X}$. Then $\mathcal{X}$ can be interpreted as the set of
vertices of a graph with directed edges $(x,y)$ for those
$(x,y) \in \mathcal{X} \times \mathcal{X}$ with positive weight $Q(x,y)$.
We assume the Markov chain to be irreducible or equivalently the
corresponding graph to be strongly connected. Thus, there exists a
unique stationary distribution $\pi : \mathcal{X} \to (0,1]$ of the
Markov chain with $\sum_{x \in \mathcal{X}} \pi(x) = 1$. We further
assume that the Markov chain is reversible with respect to $\pi$,
i.e.~the detailed balance condition $\pi(x) Q(x,y) = \pi(y) Q(y,x)$
holds for all $x,y \in \mathcal{X}$.  Now, the set of probability
densities on $\mathcal{X}$ with respect to $\pi$ is given by
\begin{align*}
 \mathcal{P}(\mathcal{X}) := \left\{ \rho: \mathcal{X} \to \mathbb{R}^+_0 : \sum_{x \in \mathcal{X}} \pi(x) \rho(x) = 1 \right\} \, .
\end{align*}
For brevity, in the following we will write $\mathbb{R}^{\mathcal{X}}$ and $\mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ for the spaces of real functions over $\mathcal{X}$ and $\mathcal{X} \times \mathcal{X}$ respectively.

Next, we define the following inner products on $\mathbb{R}^{\mathcal{X}}$ and $\mathbb{R}^{\mathcal{X} \times \mathcal{X}}$
\begin{align}
 \label{def:innerProducts}
 \InProdV{\phi}{\psi} & := \sum_{x \in \mathcal{X}} \phi(x) \psi(x) \pi(x), &
 \InProdE{\Phi}{\Psi} & := \frac{1}{2} \sum_{x,y \in \mathcal{X}} \Phi(x,y) \Psi(x,y) Q(x,y) \pi(x)
 \end{align}
 for $\phi$, $\psi \in \mathbb{R}^{\mathcal{X}}$ and $\Phi$, $\Psi \in \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$.
The corresponding induced norms are denoted by $\NormV{\cdot}$ and $\NormE{\cdot}$.
A discrete gradient $\nabla_{\mathcal{X}} : \mathbb{R}^{\mathcal{X}} \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ and a discrete divergence 
$\div_{\mathcal{X}} : \mathbb{R}^{\mathcal{X} \times \mathcal{X}} \to \mathbb{R}^{\mathcal{X}}$ are given by
\begin{align}
 \label{def:graddiv}
 (\nabla_{\mathcal{X}} \psi)(x,y) & := \psi(x) - \psi(y), &
 (\div_{\mathcal{X}} \Psi)(x) & := \frac{1}{2} \sum_{y \in \mathcal{X}} Q(x,y) (\Psi(y,x) - \Psi(x,y)).
 \end{align}
Then the duality between these two operators formulated as  the discrete integration by parts formula
\begin{align*}
  \InProdV{\phi}{\div_{\mathcal{X}} \Psi} = -\InProdE{\nabla_\mathcal{X} \phi}{\Psi}
\end{align*}
can easily be verified.
The associated discrete Laplace-operator $\Delta_{\mathcal{X}} : \mathbb{R}^{\mathcal{X}} \to \mathbb{R}^{\mathcal{X}}$ is given by
\begin{align*}
 \Delta_{\mathcal{X}} \psi (x) 
 := \div_{\mathcal{X}} (\nabla_{\mathcal{X}} \psi) (x)
 = \sum_{y \in \mathcal{X}} Q(x,y)\left[ \psi(y)-\psi(x) \right]
 = (Q-D)\psi(x) \, ,
\end{align*}
where $D=\text{diag}(\sum_yQ(x,y))_{x\in\mathcal{X}}$.  The graph divergence allows
to formulate a continuity equation for time-dependent probability
densities $\rho : [0,1] \to \mathbb{R}^{\mathcal{X}}$ and momenta
$m : [0,1] \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ describing the
flow of mass along the graph edges. In explicit, we consider the
following definition of solutions to the continuity equation with
boundary values at time $t=0$ and $t=1$.
\begin{Definition}[Continuity equation]\label{def:CEStrong}
The set $\mathcal{CE}(\rho_A, \rho_B)$ of solutions of the continuity equations for given 
boundary data $\rho_A$, $\rho_B \in \mathcal{P}(\mathcal{X})$ is defined as the set of all pairs $(\rho,m)$ 
with $\rho:[0,1] \times \mathbb{R}^\mathcal{X} \to \mathbb{R}$ and $m:[0,1] \times \mathbb{R}^{\mathcal{X} \times \mathcal{X}} \to \mathbb{R}$ measurable, such that  
\begin{equation}
\int_0^1 \InProdV{\partial_t \varphi(t,\cdot)}{\rho(t,\cdot)} +
 \InProdE{\nabla_\mathcal{X} \varphi(t,\cdot)}{m(t,\cdot)}   \d t= 
\InProdV{\varphi(1,\cdot)}{\rho_B} - \InProdV{\varphi(0,\cdot)}{\rho_A}
\label{eq:CEStrong}
\end{equation}
for all $\varphi \in C^1([0,1],\mathbb{R}^\mathcal{X})$. 
\end{Definition}

For $m\in L^2((0,1),\mathbb{R}^{\mathcal{X} \times \mathcal{X}})$ (see Lemma \ref{Lemma:L2}) 
one gets $\rho \in H^{1,2}((0,1),\mathbb{R}^{\mathcal{X}})$ and thus 
$\partial_t \rho + \div_\mathcal{X} m =0$ holds a.e.~. 
Furthermore, $\rho \in C^{0,\frac12}([0,1],\mathbb{R}^{\mathcal{X} \times \mathcal{X}})$ and 
$\rho(1,\cdot) =  \rho_B$, $\rho(0,\cdot) =  \rho_A$.
If $\rho(t,\cdot) \geq 0$ is ensured for all $t\in (0,1)$ via a finite energy property (see \eqref{def:localAction} below), then 
testing with $\varphi(t,x) = \zeta(t)$ implies that $\rho(t,\cdot) \in \mathcal{P}(\mathcal{X})$.

The Benamou--Brenier formula \cite{BeBr00} asserts that the squared
$L^2$-Wasserstein distance for probability measures in $\mathbb{R}^n$ is the
minimum of an action functional over solutions to the corresponding
continuity equation. Formally the action functional can be interpreted
as a Riemannian path length \cite{Ot01}. To construct an analogous
action functional for solutions
$(\rho,m) \in \mathcal{CE}(\rho_A,\rho_B)$ a mass
density on edges has to be deduced from the the mass densities on the
edge nodes. To this end, one defines an averaging function
$\theta :(\mathbb{R}^+_0)^2 \to \mathbb{R}^+_0$ which satisfies: \smallskip

$\theta$ is continuous, concave, 1-homogeneous, and symmetric, $\theta$ is $C^\infty$ on $(0,+\infty)^2$,
$\theta(0,s) = \theta(s,0) = 0$ and $\theta(s,s)=s$ for $s \in \mathbb{R}^+_0$, $\theta(s,t) > 0$ if $s>0$ and $t>0$, and
$s \mapsto \theta(t,s)$ is monotone increasing on $\mathbb{R}^+_0$ for fixed $t \in \mathbb{R}^+_0$.
\smallskip

It will be useful to consider $\theta$ as a concave function
$\mathbb{R}^2 \to \R \cup \{-\infty\}$. Therefore, we will set $\theta(s,t) = -\infty$
when $\min\{s,t\}<0$.  Possible choices for $\theta$ are for example
the logarithmic mean $\theta_{\tn{log}}$ or the geometric mean
$\theta_{\tn{geo}}$ for $s$, $t \in \mathbb{R}^+_0$:
\begin{align}
\label{def:means}
 \theta_{\tn{log}}(s,t) =
\begin{cases} 0, & \text{ if } s = 0 \text{ or } t = 0 \\
              s, & \text{ if } s=t \\
              \frac{t-s}{\log(t) - \log(s)} & \text{ otherwise }
\end{cases}
\quad  , \qquad \qquad
 \theta_{\tn{geo}}(s,t) = \sqrt{st} \, .
\end{align}
Note that the arithmetic mean is not admissible.  Based on this
averaging function one can define the discrete transportation distance
on $\mathcal{P}(\mathcal{X})$.
\begin{Definition}[Action functional and distance]\label{defi:definition_of_metric}
The action functional for measurable functions $\rho : [0,1] \to \mathbb{R}^\mathcal{X}$
and $m : [0,1] \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ is defined as
\begin{multline} 
	\label{def:localAction}
	\mathcal{A}(\rho,m) = \frac12 \int_0^1 \sum_{x,y \in \mathcal{X}} \alpha\big( \rho(t,x), \rho(t,y), m(t,x,y) \big) \, Q(x,y) \, \pi(x) \d t
	\\
\tn{with } \alpha  : \mathbb{R}^3 \to \R \cup \{\infty\};  \quad
 (s,t,m)  \mapsto 
 \begin{cases}
 \frac{m^2}{\theta(s,t)} & \tn{if } \theta(s,t) > 0, \\
 0 & \tn{if } \theta(s,t) = 0 \tn{ and } m = 0, \\
 + \infty & \tn{else.}
 \end{cases}
\end{multline}
The energy is then given by
\begin{align*}
	\mathcal{E}(\rho, m) &
	= \mathcal{A}(\rho, m)
	+ {\mathcal{I}}_{\mathcal{CE}(\rho_A,\rho_B)}(\rho, m)\,,
\end{align*}
where ${\mathcal{I}}_{\mathcal{CE}(\rho_A,\rho_B)}$ is the
indicator functional, which is zero for $(\rho,m)$ in
$\mathcal{CE}(\rho_A,\rho_B)$ and $\infty$ otherwise.  The induced
discrete transportation distance is obtained by
\begin{equation}\label{def:metric}
 \mathcal{W}(\rho_A, \rho_B) = \sqrt{\inf \mathcal{E}(\rho, m)}\,.
\end{equation}
\end{Definition}

Note that $\alpha$ is convex and lower semi-continuous and
$\mathcal{CE}(\rho_A,\rho_B)$ is a convex set. Hence,
\eqref{def:metric} is a convex optimization problem.  In is shown in
\cite[Theorem 3.8]{Ma11} that the mapping
$\mathcal{W}: \mathcal{P}(\mathcal{X}) \times \mathcal{P}(\mathcal{X}) \to \mathbb{R}$ defines a metric on
$\mathcal{P}(\mathcal{X})$, provided
\begin{align*}
  \int_0^1\frac{1}{\sqrt{\theta(1-r,1+r)}}\d r <\infty\;.
\end{align*}
This is the case for the logarithmic mean $\theta_{\log}$ and the
geometric mean $\theta_{\rm geo}$. In \cite[Theorem
3.2]{ErbarMaas-Ricci2012} it is shown that the infimum in
\eqref{def:metric} is attained by an optimal pair $(\rho,\psi)$. The
curve $(\rho_t)_{t\in[0,1]}$ is a constant speed geodesic for the
distance $\mathcal{W}$, i.e. it holds
$\mathcal{W}(\rho_t,\rho_s)=|t-s|\mathcal{W}(\rho_A,\rho_B)$ for all $s,t\in[0,1]$.

\subsection{A priori bounds}
In what follows we will investigate the numerical approximation of
$\mathcal{W}$ using a suitable Galerkin discretization in time and solving the
resulting discrete convex optimization problem.  Here the the
nonlinear averaging function $\theta$ and the resulting coupling of
the values of the probability desity on neighbouring nodes will
require special treatment in order to obtain a robust and effective
solution scheme.
To this end, we first discuss a few simplifications of the optimization problem \eqref{def:metric} that will help to reduce the computational complexity.
\begin{Remark}[Sparsity of kernel $Q$]
\label{rem:Sparsity}
Let $\mathcal{S} = \{ (x,y) \in \mathcal{X}^2 \,\colon\, Q(x,y)>0 \}$ be the set of `edges' indicated by non-zero transition probability. 
As $Q$ is reversible, one finds $(x,y) \in \mathcal{S}$ iff $(y,x) \in \mathcal{S}$. 
Furthermore,  
$\div_\mathcal{X} m(t,\cdot)$ and $\mathcal{A}(\rho,m)$ 
for $m : [0,1] \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ only depend on values of $m(t,x,y)$ where $(x,y) \in \mathcal{S}$.
Hence,if the kernel $Q$ is sparse, i.e.~if $\mathcal{S}$ is only a small subset of $\mathcal{X} \times \mathcal{X}$ this implies a considerable reduction of computational complexity.
\end{Remark}
In addition, the following Lemma allows to replace the two variables $m(t,x,y)$ and $m(t,y,x)$ by one effective variable, further reducing the problem size.
\begin{Lemma}[Antisymmetry of optimal
  momentum] \label{Lem:AntisymmetryM} If
  $\mathcal{W}(\rho_A,\rho_B)$ is finite and if
  $\rho:[0,1]\to\mathbb{R}^\mathcal{X}$ and
  $m : [0,1] \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ are optimal
  for \eqref{def:metric} then $m(t,x,y) = -m(t,y,x)$
  $t$-almost everywhere, whenever $(x,y) \in \mathcal{S}$ (see above
  remark for definition of $\mathcal{S}$).
\end{Lemma}
\begin{proof}
Let $\rho : [0,1] \to \mathbb{R}^\mathcal{X}$ and $m : [0,1] \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ be given such that $\mathcal{E}(\rho,m) < \infty$. Now set
\begin{align*}
\hat{m}(t,x,y) := -m(t,y,x)\,.
\end{align*}
One quickly verifies that $\div_\mathcal{X} \hat{m} = \div_\mathcal{X} m$ and that thus $(\rho,\hat{m}) \in \mathcal{CE}(\rho_A,\rho_B)$ as well.
Besides, by using that $Q(x,y)\,\pi(x) = Q(y,x)\,\pi(y)$ and $\alpha(s,t,m)=\alpha(t,s,-m)$ one finds that $\mathcal{A}(\rho,\hat{m})=\mathcal{A}(\rho,m)$.
Let now $\ol{m}=\frac12 (m + \hat{m})$. Note that $\ol{m}(t,x,y)$ is anti-symmetric in $x$ and $y$. By convexity of $\mathcal{CE}(\rho_A,\rho_B)$ one gets $(\rho,\ol{m}) \in \mathcal{CE}(\rho_A,\rho_B)$ and by convexity of $\mathcal{A}$ one finds
\begin{align*}
\mathcal{A}(\rho,\ol{m}) \leq \frac12 \left(
\mathcal{A}(\rho,m) + \mathcal{A}(\rho,\hat{m})\right)
= \mathcal{A}(\rho,m)\,.
\end{align*}
Further, the finiteness of $\mathcal{A}(\rho,m)$ implies that $m(t,x,y)=0$ when $\theta(\rho(t,x),\rho(t,y))=0$ 
and $(x,y) \in \mathcal{S}$ $t$-almost everywhere , values of $m(t,x,y)$ for $(x,y) \notin \mathcal{S}$ will have no impact on $\mathcal{A}$,
and the function $\mathbb{R} \ni z \mapsto \alpha(s,t,z)$ is even strictly convex for fixed $s,\,t > 0$.
Hence, we observe that $\mathcal{A}(\rho,\ol{m})<\mathcal{A}(\rho,m)$ unless $\ol{m}$ already coincides with $m$ for almost every $t$ and all $(x,y) \in \mathcal{S}$.
\end{proof}
In the $\Gamma$-convergence analysis we will make use on the following $L^2$ bound for the momentum.
Let us introduce the constants
\begin{align}\nonumber
  C^*:=&\max_{x\in\mathcal{X}} \sum_yQ(x,y)\;,\\\label{eq:constantsC}
  C_*:=&\min_{x,y\in\mathcal{X}, Q(x,y)>0}Q(x,y)\pi(x)\;.
\end{align}
\begin{Lemma}[$L^2$ bound for the momentum]\label{Lemma:L2}
Let $(\rho,m): [0,1] \to \mathbb{R}^\mathcal{X} \times \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ 
be a measurable path with energy  $\mathcal{E}(\rho,m) \leq \widebar E<\infty$. Then, $m$ and $\rho$ are uniformly bounded in 
$L^2((0,1), \mathbb{R}^{\mathcal{X} \times \mathcal{X}})$ and 
$H^{1,2}((0,1),\mathbb{R}^\mathcal{X}) \cap C^{0,\frac12}([0,1],\mathbb{R}^\mathcal{X})$, respectively, with 
bounds solely depending on $\mathcal{X}$ and $\widebar E$.
\end{Lemma}
\begin{proof}
   Since $\mathcal{E}(\rho,m) < \infty$, 
   we have $(\rho,m) \in \mathcal{CE}(\rho_A, \rho_B)$, 
   and thus for a.e.~$t \in (0,1)$ the mass is preserved, \emph{i.e. }  $\sum_{x \in \mathcal{X}} \rho(t,x) \pi(x) = \sum_{x \in \mathcal{X}} \rho_A(x) \pi(x) = 1$. 
   In addition, $\rho(t,x)$ is non-negative for all $x \in \mathcal{X}$ and a.e.~$t \in (0,1)$.
   By symmetry and concavity of $\theta$ and since $\theta(s,s)=s$, we can estimate 
   \begin{align*}
     \theta(\rho(t,x), \rho(t,y) ) &= \frac12\theta( \rho(t,x),\rho(t,y))+\frac12\theta(\rho(t,y),\rho(t,x)) \\
     &\leq \theta\Big(\frac{\rho(t,x)+\rho(t,y)}{2},\frac{\rho(t,x)+\rho(t,y)}{2}\Big) = \frac{\rho(t,x) + \rho(t,y)}{2}
   \end{align*}
   and get
  \begin{align}\nonumber 
   \sum_{x,y \in \mathcal{X}} \theta( \rho(t,x), \rho(t,y)) Q(x,y) \pi(x) 
   \leq \frac12\sum_{x,y \in \mathcal{X}} (\rho(t,x) Q(x,y) \pi(x) + \rho(t,y) Q(y,x) \pi(y) ) \\
=  \frac12\sum_{x,y \in \mathcal{X}} (\rho(t,x) Q(y,x) \pi(y) + \rho(t,y) Q(x,y) \pi(x) )   =   C^*\, \sum_{x \in \mathcal{X}} \rho(t,x) \pi(x) =  C^*\, . \label{eq:ThetaMonotonictyArg}
  \end{align}
  Thus, using the Cauchy--Schwarz inequality we obtain
  \begin{align}\nonumber
  \Big( \sum_{x,y \in \mathcal{X}} |m(t,x,y)| Q(x,y) \pi(x) \Big)^2
  \leq & \Big( \sum_{x,y \in \mathcal{X}} \alpha( \rho(t,x),\rho(t,y), m(t,x,y) ) Q(x,y) \pi(x) \Big) \\
  & \cdot\Big( \sum_{x,y \in \mathcal{X}} \theta(\rho(t,x), \rho(t,y) Q(x,y) \pi(x) \Big) \, . \label{eq:L2BoundM}
  \end{align}
Integrating in time we obtain 
  \begin{align*}
  \int_0^1||m(t,\cdot,\cdot)||^2_Q\d t = \int_0^1 \sum_{x,y \in \mathcal{X}} m(t,x,y)^2 Q(x,y) \pi(x) \d t \leq  \frac{C^*}{C_*} \, \widebar E \,.
  \end{align*}
Finally,  using the continuity equation \eqref{eq:CEStrong} and $m$ in $L^2((0,1),\mathbb{R}^{\mathcal{X} \times \mathcal{X}})$
we obtain that 
\begin{align*}
  \int_0^1||\partial_t\rho||^2_\pi\d t \leq \int_0^1\sum_x\Big|\sum_ym(t,x,y)Q(x,y)\Big|^2\pi(x)
d t \leq C^*\int_0^1\sum_{x,y}m(t,x,y)^2Q(x,y)\pi(x)d t\;.
\end{align*}

This implies that $\rho \in H^{1,2}((0,1),\mathbb{R}^\mathcal{X})$ and via
the Sobolev embedding theorem we obtain that also
$\rho \in C^{0,\frac12}((0,1),\mathbb{R}^\mathcal{X})$.
\end{proof}

\section{Discretization}\label{sec:Discretization}
\subsection{Galerkin discretization}\label{sec:Galerkin}
To approximate the minimizers of \eqref{def:metric} numerically we
choose a Galerkin discretization in time.  The time interval $[0,1]$
is divided into $N$ subintervals $\interval{i}=[t_i,t_{i+1})$ for
$i=0, \ldots, N-1$ with uniform step size $h=\frac{1}{N}$ and
$t_i=i\,h$.  Then, we define discrete spaces
\begin{align*}
 V_{n,h}^{1} & = \{ \psi_h \in C^0([0,1],\mathbb{R}^{\mathcal{X}})  \,\colon\,
 	\psi_h(\cdot)|_{\interval{i}} \tn{ is affine } \forall i=0,\ldots,N-1
 	\} \, , \\
 V_{n,h}^0 & = \{ \psi_h : [0,1] \to \mathbb{R}^{\mathcal{X}} \,\colon\,
 	\psi_h(\cdot)|_{\interval{i}} \tn{ is constant } \forall i=0, \ldots, N-1
 	\} \, , \\
 V_{e,h}^0 & = \{ \psi_h : [0,1] \to \mathbb{R}^{\mathcal{X} \times \mathcal{X}} \, \colon \,
 	\psi_h(\cdot)|_{\interval{i}} \tn{ is constant } \forall i=0, \ldots, N-1 \}\,.
\end{align*}
For a function $\psi_h \in V_{n,h}^0$ or $V_{e,h}^0$ we will often write $\psi_h(t_i)$ to refer to its value on the interval $\interval{i}=[t_i,t_{i+1})$.
For a function $\psi_h \in V_{n,h}^{1}$ the time-derivative can be interpreted as map
\begin{align*}
	\partial_t : V_{n,h}^{1} & \to V_{n,h}^0\;, &
	(\partial_t \psi_h)(t_i) & =\frac{1}{h}(\psi_h(t_{i+1})-\psi_h(t_i))
	\tn{ for } i=0,\ldots,N-1\,.
\end{align*}


We pick $V_{n,h}^{1} \times V_{e,h}^0$ as the space for discretized masses and momenta $(\rho_h,m_h)$. That is, discrete masses $\rho_h$ are continuous and piecewise affine and the corresponding momenta $m_h$ will be piecewise constant. $\partial_t \rho_h$ and $\div_\mathcal{X} m_h$ then lie in $V_{n,h}^0$.
In analogy to Definition \ref{def:CEStrong} we define discrete solutions of the continuity equation.
\begin{Definition}
The set of solutions to the discretized continuity equation for given boundary values $\rho_A, \rho_B \in \mathbb{R}^\mathcal{X}$ is given by
 \begin{multline}\label{eq:CEDisc}
   \mathcal{CE}_h(\rho_A, \rho_B)
   = \Big\{ (\rho_h, m_h) \in V_{n,h}^{1} \times V_{e,h}^0 \; : \; 
      h \sum_{i=0}^{N-1} \InProdV{\partial_t \rho_h(t_i,\cdot) + \div_\mathcal{X} m_h(t_i,\cdot)}{\varphi_h(t_i,\cdot)} = 0 \; \forall \varphi_h \in V_{n,h}^0 \, , \\
      \rho_h(t_0, x) = \rho_A(x) \, , \rho_h(t_N, x) = \rho_B(x)
      \Big\} \, .
\end{multline}
\end{Definition}
One can quickly verify that
$\mathcal{CE}_h(\rho_A,\rho_B) = \mathcal{CE}(\rho_A,\rho_B)
\cap (V_{n,h}^{1} \times V_{e,h}^0)$
and that $\partial_t \rho_h + \div_\mathcal{X} m_h=0$
holds for a.e.~$t$ when
$(\rho_h,m_h) \in \mathcal{CE}_h(\rho_A,\rho_B)$.
Next, we define a fully discrete action functional in analogy to
Definition~\ref{defi:definition_of_metric} and subsequently a discrete
version of the transport metric $\mathcal{W}$.

\begin{Definition}[Time-discrete action and transportation distance]
The averaging operator $\avg_h$ takes a measure $\psi \in \mathcal{M}([0,1],\mathbb{R}^{\mathcal{X}})$ to its average values on time intervals $\interval{i}$:
\begin{align*}
	\avg_h : \mathcal{M}([0,1],\mathbb{R}^{\mathcal{X}}) & \to V_{n,h}^0\,, \quad
	(\avg_h \psi)(t_i)  = \psi(\interval{i})\;
	\tn{ for } i=0,\ldots,N-1\,.
\end{align*}
Analogously we declare the $\avg_h$ operator for $\mathbb{R}^{\mathcal{X} \times \mathcal{X}}$-valued measures. Note that for $\psi_h \in V_{n,h}^{1}$ one finds
$(\avg_h \psi_h)(t_i)=\tfrac12(\psi_h(t_i)+\psi_h(t_{i+1}))\,$.
For $(\rho,m) \in \mathcal{M}([0,1],\mathbb{R}^{\mathcal{X}}) \times \mathcal{M}([0,1],\mathbb{R}^{\mathcal{X} \times \mathcal{X}}) $ the discrete approximation for the action is given by
\begin{align*}
	\mathcal{A}_h(\rho, m) & =
	\mathcal{A}(\avg_h \rho, \avg_h m) \\
	&=\frac{h}{2} \sum_{i=0}^{N-1} \sum_{x,y \in \mathcal{X}} \alpha \left( \avg_h \rho(t_i,x) , \avg_h \rho(t_i,y), \avg_h m(t_i,x,y) \right) Q(x,y) \pi(x) \, .
\end{align*}
Finally, the time discrete energy functional is defined by
$ \mathcal{E}_h(\rho, m) =
\mathcal{A}_h(\rho,m) +
{\mathcal{I}}_{\mathcal{CE}_h(\rho_A,\rho_B)}(\rho,
m) $
and for the associated time discrete approximation of the
transportation distance one obtains
\begin{align}\label{def:metricDiscrete}
 \mathcal{W}_h(\rho_A, \rho_B) = \sqrt{\inf \mathcal{E}_{h}(\rho, m)}\,.
\end{align}
\end{Definition}

Note that the indicator function of the discrete continuity equation entails the constraint $(\rho,m) \in V_{n,h}^{1} \times V_{e,h}^0$.
These spaces can be represented by finite-dimensional vectors, the operators $\partial_t$ and $\avg_h$ can be represented as finite-dimensional matrices and the continuity equation becomes a finite-dimensional affine constraint. 
Thus, \eqref{def:metricDiscrete} is indeed a finite-dimensional convex optimization problem. Its numerical solution by using proximal mappings will be detailed in Section~\ref{sec:ProximalSplitting}.

\subsection[Gamma-convergence]{$\Gamma$-convergence}

In the following, we will prove a $\Gamma$-convergence result of the discrete energy functional, which will justify our discretization.
First, we construct explicitly continuous and discrete trajectories between an arbitrary probability distribution on $\mathcal{X}$ and the uniform probability density $\odens \in \mathcal{P}(\mathcal{X})$ given by $\odens(x) = 1$. We show that these trajectories have uniformly bounded energy, which will be essential in the $\Gamma$-$\limsup$ inequality in Theorem~\ref{thm:GammaConvergence}. Let us define the 
Lagrange interpolation operator $\lindiscr_h: C^0([0,1],\mathbb{R}^\mathcal{X}) \to V_{n,h}^{1}; \rho \mapsto \lindiscr_h(\rho)$ given by
\begin{align*}
\left(\lindiscr_h \rho \right)(t_i,x) := \rho(t_i,x) \quad \forall i=0,\ldots,N \, .
\end{align*}

\begin{Proposition}
	\label{prop:FiniteCostToFlatMeasure}
	There is some constant $C(\mathcal{X}) < \infty$ such that for
        any $\rho_A \in \mathcal{P}(\mathcal{X})$ there is a trajectory
        $(\rho,m) \in \mathcal{CE}(\rho_A,\odens)$ with
        $\mathcal{A}(\rho,m) \leq C(\mathcal{X})$ and
        $(\lindiscr_h\rho,\avg_hm) \in
        \mathcal{CE}_h(\rho_A,\odens)$
        with
        $\mathcal{A}_h(\lindiscr_h\rho,\avg_hm)) \leq
        C(\mathcal{X})$ for every $h=1/N$.
\end{Proposition}
\begin{proof}
For $x \in \mathcal{X}$ let $\rho^x_A\in \mathcal{P}(\mathcal{X})$ be the probability density on $\mathcal{X}$ with all mass concentrated on $x$. That is, 
$\rho^x_A =  \tfrac1{\pi(x)} \delta_x$, where $\delta_x$ is the usual Kronecker symbol with $\delta_x(y) = 1$ if $x=y$ and $0$ else.
\medskip

\emph{Construction of elementary flows:}
For $(x,y) \in \mathcal{X} \times \mathcal{X}$, $x\neq y$, with $Q(x,y)>0$ we define $L[x,y] \in \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ as follows:
\begin{align*}
L[x,y](a,b) = \begin{cases}
\frac{1}{Q(x,y)\,\pi(x)} & \tn{if } (a,b) = (x,y), \\
\frac{-1}{Q(x,y)\,\pi(x)} & \tn{if } (a,b) = (y,x), \\
0 & \tn{else.}
\end{cases}
\end{align*}
Then $\div_\mathcal{X} L[x,y]=\rho^y_A -\rho^x_A$.
Now, for any $(x,y) \in \mathcal{X} \times \mathcal{X}$, $x \neq y$, there exists a path $(x=x_0,x_1,\ldots,x_K=y)$ with $K<\mathrm{card}\, \mathcal{X}$ with $Q(x_k,x_{k+1})>0$ for $k=0,\ldots, K-1$. We can add the corresponding $L(x_k,x_{k+1})$ along these edges to construct a flow $M[x,y]$ with $\div_\mathcal{X} M[x,y]=\rho^y_A -\rho^x_A$.
All entries of all $M[x,y]$ are bounded (in absolute value) by $\widetilde C(\mathcal{X}):=\mathrm{card}\mathcal{X}/C_*$, where $C_*$ is defined in \eqref{eq:constantsC}. For $x=y$, $M[x,x]$ is simply zero.

Now assume $\rho_A=\rho_A^x$ for some $x \in \mathcal{X}$. Let
$m_0 = \sum_{y \in \mathcal{X}} M[x,y]\,\pi(y)\,$. One finds
\begin{align*}
\div_\mathcal{X} m_0 = \sum_{y \in \mathcal{X}} \left(\tfrac1{\pi(y)}  \delta_y -\tfrac1{\pi(x)}  \delta_x\right) \pi(y) = \odens-\rho_A^x\,.
\end{align*}
Again, every entry of $m_0$ is bounded in absolute value by $\widetilde C(\mathcal{X})$. Now let
$m(t) =  2\,m_0\,t$,
$\rho(t) =  \rho_A^x + (\div_\mathcal{X} m_0)\,t^2 = (1-t^2) \cdot \rho_A^x  + t^2 \cdot \odens\,$.
We find $(\rho,m) \in \mathcal{CE}(\rho_A^x,\odens)$.
One has $|m(t,x,y)| \leq t \cdot 2\widetilde C(\mathcal{X})$ and $\rho(t,x)\geq t^2$ and using the monotonicity of $\alpha$ for the action $\mathcal{A}$ we get
\begin{align*}
\mathcal{A}(\rho,m) & \leq \frac12 \int_0^1 \sum_{x,y \in \mathcal{X}} \frac{(t \cdot 2 \widetilde C(\mathcal{X}))^2}{t^2} Q(x,y)\,\pi(x)\, \d t = 2\widetilde C(\mathcal{X})^2C^*\,.
\end{align*}
\medskip

\emph{Construction of discrete counterparts:}
For fixed $h=1/N$ let $\rho_h=\lindiscr_h \rho$ and $m_h = \avg_h m$.
By construction $(\rho_h,m_h) \in \mathcal{CE}_h(\rho_A^x,\odens)$. Then, one finds
$m_h(t_i,x,y)  \leq (i+\tfrac12)\,h\,2\widetilde C(\mathcal{X})$,
$\rho_h(t_i,x)  \geq i^2\,h^2$, 
$(\avg_h \rho_h)(t_i,x)  \geq (i^2 + i + \tfrac12)\,h^2\,$,
and thus
\begin{align*}
\mathcal{A}_h(\rho_h, m_h) & = \mathcal{A}(\avg_h \rho_h,m_h) \leq \frac12 \sum_{i=0}^{N-1}
h \frac{h^2\,4\widetilde C(\mathcal{X})^2\,(i+\tfrac12)^2}{h^2 (i^2 + i + \tfrac12)}
\sum_{x,y \in \mathcal{X}} Q(x,y)\,\pi(x)
\leq 2\widetilde C(\mathcal{X})^2C^*\,.
\end{align*}
\medskip

\emph{Extension for arbitrary initial data:}
For given $x \in \mathcal{X}$ let $(\rho^x,m^x)$ be the (continuous) trajectory between $\rho_A^x$ and $\odens$ as constructed above. Any $\rho_A$ is a superposition of various $\rho^x_A$:
\begin{align*}
\rho_A = \sum_{x \in \mathcal{X}} \rho_A(x) \,\delta_x = \sum_{x \in \mathcal{X}} \rho_A(x) \pi(x) \,\rho^x_A
 \end{align*}
By linearity of the continuity equation the trajectory 
$(\rho,m) = \sum_{x \in \mathcal{X}} \rho_A(x)\,\pi(x) \cdot (\rho^x,m^x)
$
then lies in $\mathcal{CE}(\rho_A,\odens)$.
Since $\mathcal{A}$ is convex and 1-homogeneous, it is sub-additive. Therefore,
\begin{align*}
\mathcal{A}(\rho,m) & \leq \sum_{x \in \mathcal{X}} \rho_A(x)\,\pi(x)\,\mathcal{A}(\rho^x,m^x) \leq \sum_{x \in \mathcal{X}} \rho_A(x)\,\pi(x)\,2\,\widetilde C(\mathcal{X})^2C^*=2\,\widetilde C(\mathcal{X})^2C^*\,.
\end{align*}
For the discrete trajectory the reasoning is completely
analogous. Thus the claim follows with
$C(\mathcal{X})=2\,\widetilde C(\mathcal{X})^2C^*$.
\end{proof}

\begin{Corollary}
	\label{cor:UniformWcBound}
	The above strategy can be used to construct trajectories between arbitrary $\rho_A$, $\rho_B$ via $\odens$ as intermediate state. This establishes that $\mathcal{W}$ and $\mathcal{W}_h$ are uniformly bounded on $\mathcal{P}(\mathcal{X})^2$.
\end{Corollary}
\begin{Remark}
	In \cite{Ma11} it is shown that $\mathcal{W}$ is bounded if the constant
	$
		C_\theta := \int_0^1 \frac{1}{\sqrt{\theta(1-r,1+r)}} \d r
	$
	is finite. 
	Here, we assumed that $\theta(s,s)=s$ for $s \in \mathbb{R}^+_0$ and that  $s \mapsto \theta(s,t)$  is increasing on $\mathbb{R}^+_0$ for fixed $t \in \mathbb{R}^+_0$ 
	which implies that $\theta(s,t) \geq \min\{s,t\}$ for $s,\,t \in \mathbb{R}^+_0$. This is sufficient for $C_\theta<\infty$.
\end{Remark}
\begin{Theorem}[$\Gamma$-convergence of time discrete energies]
	\label{thm:GammaConvergence}
	Let $\rho_A$, $\rho_B$ be fixed temporal boundary conditions. Then, the sequence of functionals $(\mathcal{E}_h)_h$ $\Gamma$-converges for $h \to 0$ to the functional $\mathcal{E}$ with respect to the weak$\ast$ topology in $\mathcal{M}([0,1],\mathbb{R}^{\mathcal{X}} \times \mathbb{R}^{\mathcal{X} \times \mathcal{X}})$.
\end{Theorem}
\begin{proof}
To establish $\Gamma$-convergence, we have to verify the $\Gamma$-$\liminf$ and $\Gamma$-$\limsup$ properties.
\medskip

For the $\Gamma$-$\liminf$ property, we have to demonstrate that  the inequality
\begin{align} \label{eq:liminf}
\mathcal{A}(\rho,m)
+ {\mathcal{I}}_{\mathcal{CE}(\rho_A,\rho_B)}(\rho,m)
\leq \liminf_{h \to 0}
\mathcal{A}_h(\rho_h,m_h)
+ {\mathcal{I}}_{\mathcal{CE}_h(\rho_A,\rho_B)}(\rho_h,m_h)
\end{align}
holds for all sequences $(\rho_h,m_h) \stackrel{*}{\rightharpoonup} (\rho,m)$ in $\mathcal{M}([0,1],\mathbb{R}^{\mathcal{X}} \times \mathbb{R}^{\mathcal{X} \times \mathcal{X}})$.
As $\mathcal{CE}(\rho_A,\rho_B)$ is weak-$\ast$ closed and $\mathcal{CE}_h(\rho_A,\rho_B) \subset \mathcal{CE}(\rho_A,\rho_B)$ the statement is trivial if there is no subsequence with $(\rho_h,m_h) \in \mathcal{CE}_h(\rho_A,\rho_B)$. 
Thus, we may assume that all $(\rho_h,m_h)$ fulfill the discrete continuity equation,  that $(\rho,m)$
fulfills the continuous continuity equation, and  all $\rho_h$ are non-negative.

Now, $\rho_h \stackrel{*}{\rightharpoonup} \rho$ implies $\avg_h \rho_h \stackrel{*}{\rightharpoonup} \rho$ and
$\mathcal{A}_h(\rho_h,m_h) = \mathcal{A}(\avg_h \rho_h,m_h)\,$.
Since $\alpha$ is jointly convex and lower semi-continuous in $\rho$ and $m$, the action functional $\mathcal{A}$ is weak-$\ast$ lower semi-continuous and \eqref{eq:liminf} holds. \medskip


To verify the $\Gamma$-$\limsup$ property we need to show that for any $(\rho,m) \in \mathcal{M}([0,1],\mathbb{R}^{\mathcal{X}} \times \mathbb{R}^{\mathcal{X} \times \mathcal{X}})$ there exists a recovery sequence $(\rho_h,m_h) \stackrel{*}{\rightharpoonup} (\rho,m)$ with
\begin{align}\label{eq:limsup}
\limsup_{h \to 0}
\mathcal{A}_h(\rho_h,m_h)
+ {\mathcal{I}}_{\mathcal{CE}_h(\rho_A,\rho_B)}(\rho_h,m_h)
\leq \mathcal{A}(\rho,m)
+ {\mathcal{I}}_{\mathcal{CE}(\rho_A,\rho_B)}(\rho,m)\,.
\end{align}
We may assume $\mathcal{A}(\rho,m)<\infty$ and $(\rho,m) \in \mathcal{CE}(\rho_A,\rho_B)$. 
Using Lemma \ref{Lemma:L2} this implies in particular that $\rho \in C^{0,\frac12}([0,1],\mathbb{R}^{\mathcal{X}})$.
For such a trajectory $(\rho,m)$ we will construct a recovery sequence in two steps:
First, the continuous trajectory $(\rho,m)$ is regularized, then, the regularized still time continuous trajectory is discretized using local averaging in time.
The regularization is necessary to control the effect of the discontinuity of $\alpha$ at the origin, see \eqref{def:localAction}.

\medskip

Let $(\rho_{A,\odens},m_{A,\odens}) \in \mathcal{CE}(\rho_A,\odens)$ be the trajectory from $\rho_A$ to $\odens$ as constructed in Proposition \ref{prop:FiniteCostToFlatMeasure}, analogously let $(\rho_{\odens,B},m_{\odens,B}) \in \mathcal{CE}(\odens,\rho_B)$ be the corresponding trajectory from $\odens$ to $\rho_B$ with $(\rho_{\odens,B},m_{\odens,B})(t,\cdot)  := (\rho_{B,\odens},-m_{B,\odens})(1-t,\cdot)$. Then, 
for $\delta \in (0,\frac12)$ and $\epsilon=\delta^2$ we define
\begin{align*}
\rho_\delta(t) & = \begin{cases}
(1-\epsilon) \cdot \rho_A + \epsilon \cdot \rho_{A,\odens}(t/\delta) &
\tn{for } t \in [0,\delta), \\
(1-\epsilon) \cdot \rho((t-\delta)/(1-2\delta)) + \epsilon \cdot \odens & 
\tn{for } t \in [\delta,1-\,\delta), \\
(1-\epsilon) \cdot \rho_B + \epsilon \cdot \rho_{\odens,B}((t-(1-\delta))/\delta) & \tn{for } t \in [1-\delta,1]
\end{cases} \\
\intertext{and}
m_\delta(t) & = \begin{cases}
\epsilon \cdot \delta^{-1} \cdot m_{A,\odens}(t/\delta) & 
\tn{for } t \in [0,\delta), \\
(1-\epsilon) \cdot (1-2\delta)^{-1} \cdot m((t-\delta)/(1-2\delta)) & 
\tn{for } t \in [\delta,1-\,\delta), \\
\epsilon \cdot \delta^{-1} \cdot m_{\odens,B}((t-(1-\delta))/\delta) & \tn{for } t \in [1-\delta,1]\,.
\end{cases}
\end{align*}
One finds that $(\rho_\delta,m_\delta) \in \mathcal{CE}(\rho_A,\rho_B)$.
To evaluate the action of $(\rho_\delta,m_\delta)$ we decompose it into the contributions of the time intervals $\interval{l}=[0,\delta]$, $\interval{m}=[\delta,1-\delta]$ and $\interval{r}=[1-\delta,1]$:
\newcommand{\Ac^{\tn{int}}}{\mathcal{A}^{\tn{int}}}
\begin{align*}
\mathcal{A}(\rho_\delta,m_\delta)=\mathcal{A}_l + \mathcal{A}_m + \mathcal{A}_r
\qquad \tn{with} \qquad
\mathcal{A}_\chi = \int_{\interval{\chi}} \Ac^{\tn{int}}(\rho_\delta(t),m_\delta(t))\d t
\tn{ for } \chi \in \{l,m,r\}\,.
\end{align*}
where
\begin{align*}
\Ac^{\tn{int}} & : \mathbb{R}^{\mathcal{X}} \times \mathbb{R}^{\mathcal{X} \times \mathcal{X}} \to \R \cup \{\infty\}, &
(\rho,m) & \mapsto \frac12 \sum_{x,y \in \mathcal{X}} \alpha\big( \rho(x), \rho(y), m(x,y) \big) \, Q(x,y) \, \pi(x)\,.
\end{align*}
$\Ac^{\tn{int}}$ is jointly convex and 1-homogeneous and therefore sub-additive. Moreover, it is 2-homogeneous in the second argument. Therefore we obtain
\begin{align*}
\nonumber
\mathcal{A}_m & \leq \frac{(1-\epsilon)}{(1-2\delta)^{2}} \int_{\interval{m}}
\Ac^{\tn{int}}\big(
\rho((t-\delta)/(1-2\delta)),
m((t-\delta)/(1-2\delta))
\big)
\d t \\\label{eq:LimSupEMBound}
&= \frac{(1-\epsilon)}{(1-2\delta)}
\int_0^1 \Ac^{\tn{int}}\big(\rho(t),m(t)\big) \d t
= \frac{(1-\epsilon)}{(1-2\delta)}
\mathcal{A}(\rho,m)\,.
\end{align*}
Further, using Proposition  \ref{prop:FiniteCostToFlatMeasure} we obtain $\mathcal{A}_{l} + \mathcal{A}_{r} \leq 2 \,C(\mathcal{X}) \, \delta$.

\medskip
Next, we discretize in time. Since
$\rho \in C^{0,\frac12}([0,1],\mathbb{R}^{\mathcal{X}})$ we have
$|\rho(t,x)-\rho(t',x)|\leq g(|t-t'|)$ with
$g(s):= C \cdot s^{\frac12}$ for all $x \in \mathcal{X}$. Now let
$\Delta=g(2h)$ and choose the regularization parameter
$\delta=\min\{i \cdot h \,\colon\, i \in \mathbb{N},\, i\cdot h \geq
\Delta^{\frac14} \}$
and as before $\epsilon=\delta^2$. Obviously $\Delta$, $\delta$ and
$\epsilon \to 0$ as $h \to 0$. In particular, for $h$ sufficiently
small $2 \geq 1/(1-2\delta)$ and thus
$\Delta=g(2h)\geq g(h/(1-2\delta))$. Therefore, $\Delta$ is a uniform
upper bound for the variation of $\rho_\delta$ on any interval
of the size $h$. We now set
\begin{align*}
\rho_h  = \lindiscr_h \rho_\delta, \quad
m_h  = \avg_h m_\delta\,,
\end{align*}
and note that
$(\rho_h,m_h) \in \mathcal{CE}_h(\rho_A,\rho_B)$.  As
$\delta \to 0$ one finds
$(\rho_\delta,m_\delta) \stackrel{*}{\rightharpoonup}
(\rho,m)$
and for $h \to 0$ we obtain
$(\rho_\delta-\rho_h,m_\delta-m_h)
\stackrel{*}{\rightharpoonup} 0$.
This implies that
$(\rho_h,m_h) \stackrel{*}{\rightharpoonup} (\rho,m)$.
\medskip

Note that $\delta$ was chosen to be an integer multiple of $h$. So the division of $[0,1]$ into the three intervals $[0,\delta]$, $[\delta,1-\delta]$ and $[1-\delta,1]$ in the construction of $(\rho_\delta,m_\delta)$ is compatible with the grid discretization of step size $h$.
Therefore, as above, the discrete action decomposes into three contributions which we denote 
$\mathcal{A}_h(\rho_h,m_h)=\mathcal{A}_{l,h} + \mathcal{A}_{m,h} + \mathcal{A}_{r,h}$.
Again, using joint 1-homogeneity and sub-additivity of $\alpha$, as well as the 2-homogeneity in the second argument one obtains
\begin{align*}
\mathcal{A}_{l,h} \leq \frac{\epsilon}{\delta} \cdot \mathcal{A}_h(\lindiscr_h \rho_{A,\odens}, \avg_h m_{A,\odens}),\quad 
\mathcal{A}_{r,h} \leq \frac{\epsilon}{\delta} \cdot \mathcal{A}_h(\lindiscr_h \rho_{\odens, B}, \avg_h m_{\odens, B})\,.
\end{align*}
Using Proposition \ref{prop:FiniteCostToFlatMeasure} we observe that
\begin{align*}
\mathcal{A}_{l,h} + \mathcal{A}_{r,h} \leq 2 \,C(\mathcal{X}) \, \delta\,.
\end{align*}
In view of \eqref{eq:LimSupEMBound}, it remains to estimate
$\mathcal{A}_{m,h}$ by a suitable constant times $\mathcal{A}_m$.
\newcommand{\intervalsM}{S_m} To this end, let
$\intervalsM \subset \{0,\ldots,N-1\}$ the the set of discrete indices
such that $\interval{i} \subset \interval{m}$ for $i \in \intervalsM$.
Then $\mathcal{A}_m$ is given by
\begin{align*}
\mathcal{A}_m = \frac12 \sum_{i\in\intervalsM} \sum_{x,y \in \mathcal{X}}
\left[ \int_{\interval{i}}
\alpha\big(\rho_\delta(t,x),\rho_\delta(t,y),m_\delta(t,x,y)\big)\,\d t
\right] Q(x,y)\,\pi(x)\,.
\end{align*}
Since $\alpha$ is convex, by Jensen's inequality one finds
\begin{align*}
\int_{\interval{i}} \alpha\big(\rho_\delta(t,x),\rho_\delta(t,y),m(t,x,y)\big)\,\d t
&\geq h \cdot
\alpha\big( (\avg_h \rho_\delta)(t_i,x), (\avg_h \rho_\delta)(t_i,y), (\avg_h m_\delta)(t_i,x,y) \big) \,.
\end{align*}
The discretized action $\mathcal{A}_{m,h}$ is a weighted sum of the form
\begin{align*}
\mathcal{A}_{m,h} & =
\frac12 \sum_{i \in \intervalsM} \sum_{x,y \in \mathcal{X}} 
h \cdot \alpha \big((\avg_h \lindiscr_h \rho_\delta)(t_i,x), (\avg_h \lindiscr_h \rho_\delta)(i,y), \avg_h m_\delta)(t_i,x,y) Q(x,y)\,\pi(x) \,.
\end{align*}
By construction $\rho_\delta$ is bounded from below by $\epsilon$ on $\interval{m}$ on all nodes and its variation within each discretization interval is bounded by $\Delta$. Therefore, for any $i \in \intervalsM$, $z \in \mathcal{X}$ one finds
\begin{align*}
(\avg_h \rho_\delta)(t_i,z) & \leq (\avg_h \lindiscr_h \rho_\delta)(t_i,z) + \Delta\,, \quad 
(\avg_h \lindiscr_h \rho_\delta)(t_i,z)  \geq \epsilon\,.
\end{align*}
Due to the monotonicity of $s\to \tfrac{s}{s+\Delta}$ we obtain 
\begin{align*}
\frac{(\avg_h \lindiscr_h \rho_\delta)(t_i,z)}{(\avg_h \rho_\delta)(t_i,z)}
\geq \frac{(\avg_h \lindiscr_h \rho_\delta)(t_i,z)}{(\avg_h \lindiscr_h \rho_\delta)(t_i,z) + \Delta}
\geq \frac{\epsilon}{\epsilon + \Delta}.
\end{align*}
Taking into account the joint 1-homogeneity of $\theta$ and the monotonicity of $\theta$ in each single argument this implies for all $x,y \in \mathcal{X}$ that
\begin{align*}
\frac{\theta\big(
(\avg_h \lindiscr_h \rho_\delta)(t_i,x),
(\avg_h \lindiscr_h \rho_\delta)(t_i,y)\big)
}{
\theta\big(
(\avg_h \rho_\delta)(t_i,x),
(\avg_h \rho_\delta)(t_i,y)\big)
}
\geq \frac{\epsilon}{\epsilon + \Delta} = \frac{1}{1 + \Delta/\epsilon}\,.
\end{align*}
Hence, 
\begin{align*}
\mathcal{A}_{m,h} & =
\frac12 \sum_{i \in \intervalsM} \sum_{x,y \in \mathcal{X}} 
h \cdot \frac{(\avg_h m_\delta)^2(t_i,x,y)}{\theta((\avg_h \lindiscr_h \rho_\delta)(t_i,x), (\avg_h \lindiscr_h \rho_\delta)(t_i,y))} Q(x,y)\,\pi(x) \\
&\leq \frac12 (1 + \Delta / \epsilon) \sum_{i \in \intervalsM} \sum_{x,y \in \mathcal{X}} 
h \cdot \frac{(\avg_h m_\delta)^2(t_i,x,y)}{\theta((\avg_h  \rho_\delta)(t_i,x), (\avg_h  \rho_\delta)(t_i,y))} Q(x,y)\,\pi(x)  = (1 + \Delta / \epsilon) \mathcal{A}_{m}\,.
\end{align*}
Our choice of $\delta$ implies that $\epsilon = \delta^2 \geq \Delta^\frac12$ and thus $\Delta/\epsilon \leq \epsilon$. 
Altogether, we obtain for $h$ sufficiently small 
\begin{align*}
\mathcal{A}_h(\rho_h,m_h)  =\mathcal{A}_{l,h} + \mathcal{A}_{m,h} + \mathcal{A}_{r,h} \leq 2\, C(\mathcal{X}) \, \delta + 
(1 +  \epsilon) \, \frac{1-\epsilon}{1-2\delta} \, \mathcal{A}(\rho,m)\,.
\end{align*}
Since $\delta \to 0$, $\epsilon \to 0$ as $h \to 0$, this establishes the $\Gamma$-$\limsup$ property.
\end{proof}

Next, we establish convergence of the discrete optimizers to a
continuous solution.  To establish compactness we first show a uniform
bound for the $L^2$ norm of the discrete momenta, in analogy to Lemma
\ref{Lemma:L2}.

\begin{Lemma}[$L^2$ bound for the discrete momentum]
\label{Lemma:L2Disc}
Let $(\rho_h,m_h) \in V_{n,h}^{1} \times V_{e,h}^0$ with discrete energy $\mathcal{E}_h(\rho_h,m_h) \leq \ol{E} < \infty$.
Then, there exists a constant $\ol{M}<\infty$ only depending on $(\mathcal{X},Q,\pi)$ and $\ol{E}$ (and not on $h$), 
such that $\|m_h\|_{L^2([0,1],\mathbb{R}^{\mathcal{X} \times \mathcal{X}})} \leq \ol{M}$.
\end{Lemma}
\begin{proof}
The proof works in complete analogy to Lemma \ref{Lemma:L2}. We bound
\begin{multline*}
\Big( \sum_{x,y \in \mathcal{X}} |m_h(t_i,x,y)| Q(x,y) \pi(x) \Big)^2
\leq  \Big( \sum_{x,y \in \mathcal{X}} \alpha( \avg_h \rho_h(t_i,x),\avg_h
\rho_h(t_i,y), m(t_i,x,y) ) Q(x,y) \pi(x) \Big) \\
 	{} \cdot  \Big( \sum_{x,y \in \mathcal{X}} \theta(\avg_h \rho_h(t_i,x), \avg_h
 		\rho(t_i,y)) Q(x,y) \pi(x) \Big)
\end{multline*}
and
\begin{align*}
\sum_{x,y \in \mathcal{X}} \theta(\avg_h \rho_h(t_i,x), \avg_h
 		\rho(t_i,y))\, Q(x,y) \pi(x)
\leq  C^*\,,
\end{align*}
where $C^*$ is defined in \eqref{eq:constantsC}. Here, we have used
that $(\rho_h,m_h) \in \mathcal{CE}(\rho_A, \rho_B)$
which implies that mass is preserved, \emph{i.e. }
$\sum_{x \in \mathcal{X}} \avg_h \rho_h(t_i,x) \pi(x) = \sum_{x
  \in \mathcal{X}} \rho_h(t_i+\frac{h}{2},x) \pi(x) = \sum_{x \in
  \mathcal{X}} \rho_A(x) \pi(x) = 1$
for all $i=0,\ldots, N-1$, and that since
$\mathcal{A}_h(\rho_h,m_h)<\infty$ one has
$\avg_h \rho_h \geq 0$.  Now, once more using that $\mathcal{X}$
is finite and integrating (or summing) in time establishes the bound.
\end{proof}

\begin{Theorem}[Convergence of discrete geodesics]
\label{thm:Convergence} 
For fixed temporal boundary conditions $\rho_A$, $\rho_B$
any sequence $(\rho_h,m_h)$ of minimizers of
$\mathcal{E}_h$ is uniformly bounded in
$C^{0,\frac12}([0,1],\mathbb{R}^\mathcal{X}) \times L^2((0,1),\mathbb{R}^{\mathcal{X}
  \times \mathcal{X}})$
for $h \to 0$. Up to selection of a subsequence,
$\rho_h \to \rho$ strongly in
$C^{0,\alpha}([0,1],\mathbb{R}^\mathcal{X})$ for any $\alpha\in0,\frac12)$ and
$m_h \to m$ weakly in $L^2$ with
$( \rho,m)$ being a minimizer of the energy
$\mathcal{E}$.
\end{Theorem}
\begin{proof} 
For a sequence of minimizers $(\rho_h,m_h)$ the discrete energy $\mathcal{E}_h(\rho_h,m_h)$ is uniformly bounded by Corollary \ref{cor:UniformWcBound} .
Since $(\rho_h,m_h) \in \mathcal{CE}(\rho_A,\rho_B)$ the total variation of all $\rho_h$ is uniformly bounded. Further, by Lemma \ref{Lemma:L2Disc} the $L^2$ norm $\|m_h\|_{L^2([0,1],\mathbb{R}^{\mathcal{X} \times \mathcal{X}})}$ is uniformly bounded. 
Hence, the sequence $(\rho_h,m_h)_h$ has a weakly$\ast$ (in the sense of measures) convergent subsequence, which by Theorem \ref{thm:GammaConvergence} and a standard  consequence of  $\Gamma$ convergence theory converges weakly$\ast$ to some minimizer $(\rho,m)$ of $\mathcal{E}$).
\medskip

Using the continuity equation this convergence can be strengthened.
We already know that $(\rho_h,m_h)$ solves the continuity equation $\partial_t \rho_h = - \div_\mathcal{X} m_h$. 
Thus, the uniform bound for $m_h$ in $L^2((0,1),\mathbb{R}^{\mathcal{X} \times \mathcal{X}})$ implies that $\rho_h$ is uniformly bounded in $H^{1,2}(\mathbb{R}^\mathcal{X})$.
From this we obtain by the Sobolev embedding theorem that  $(\rho_h)_h$ is uniformly bounded in $C^{0,\frac12}(\mathbb{R}^\mathcal{X})$ and compact in $C^{0,\alpha}(\mathbb{R}^\mathcal{X})$ for all $\alpha \in [0,\frac12)$.
\end{proof}

\section{Optimization with Proximal Splitting}
\label{sec:ProximalSplitting}
\subsection{Slack Variables and Proximal Splitting}
The computation of the discrete transportation distance
\eqref{def:metricDiscrete} and the associated transport path require
the solution of a finite-dimensional non-smooth convex optimization
problem.  To this end, we apply a proximal splitting approach with
suitably choosen slack variables.  The proximal mapping of a convex
and lower semi-continuous function $f:H \to \mathbb{R} \cup \{ \infty \}$ on a
Hilbert space $H$ with norm $\NormH{\cdot}$ is defined as
 \begin{align}
  \mathrm{prox}_f(x) = \argmin_{y \in H} \frac{1}{2} \NormH{x-y}^2 + f(y)\,.
 \end{align}
Furthermore, the indicator function of a closed convex set $K \subset H$ is given by 
${\mathcal{I}}_K(x) = 0$ for $x\in K$ and $\infty$ elsewise.
In particular, $\mathrm{prox}_{{\mathcal{I}}_K} =  \mathrm{proj}_K$, where $\mathrm{proj}_K$ is the projection onto $K$. For a function $f : H \mapsto \R \cup \{\infty\}$ its Fenchel conjugate is given by
\begin{align}
	f^\ast(y) = \sup_{x \in H} \, \InProdH{y}{x} - f(x)\,.
\end{align}
If $f(x)<\infty$ for some $x \in H$, then $f^\ast$ is convex and lower semi-continuous.
For more details and an introduction to convex analysis see e.g.~\cite{ConvexFunctionalAnalysis-11}.
The practical applicability of proximal splitting schemes depends on whether the objective can be split into terms such that the proximal mapping for each term can be computed efficiently.
In \cite{PaPe14} a spatiotemporal discretization with staggered grids of the classical Benamou--Brenier formulation \cite{BeBr00} of optimal transport of Lebesgue densities on $\mathbb{R}^n$ was presented and several proximal splitting methods were considered to solve the discrete problem.
However, this approach can not directly be transfered to problem \eqref{def:metricDiscrete} since the action $\mathcal{A}$ couples the variables $\rho$ and $m$ in a non-linear way via the terms $\alpha(m(t_i,x,y), \avg_h \rho(t_i,x), \avg_h \rho(t_i,y))$ spatially over the whole graph according to the transition kernel $Q$ 
and temporally via the averaging operator $\avg_h$.
Thus, the proximal mapping of the $\mathcal{A}$-term is not separable in space or time and thus requires the solution of a fully coupled, nonlinear minimization problem.
As a remedy, we propose to introduce auxiliary variables to decouple the variables and rewrite the action $\mathcal{A}$ with terms where variables only interact locally, thus leading to separable, hence simpler, proximal mappings.

\begin{Lemma}
For $(\rho,m) \in V_{n,h}^{1} \times V_{e,h}^0$ one finds
\begin{align}
\label{eq:actionReform1}
\mathcal{A}_h(\rho,m)  = \mathcal{A}(\avg_h\rho,m) =
\inf\left\{\widehat{\mathcal{A}}(\vartheta,m) + {\mathcal{I}}_{\mathcal{K}_{\tn{pre}}}(\avg_h \rho,\vartheta)
\colon \vartheta \in V_{e,h}^0 \right\}
\end{align}
with the convex set
\begin{align}
 \mathcal{K}_{\tn{pre}}  := \left\{ (\bar \rho, \vartheta ) \in V_{n,h}^0\timesV_{e,h}^0 \, \colon \,
0 \leq \vartheta(t_i,x,y) \leq \theta(\bar \rho(t_i,x), \bar \rho(t_iy))
\, \forall i=0,\ldots,N-1, \ \forall x,y \in \mathcal{X} \right\}
\end{align}
and the edge-based action
\begin{align}
\label{def:ActionU}
\widehat{\mathcal{A}}(\vartheta,m) & := \frac12 \int_0^1 \sum_{x,y \in \mathcal{X}}
\Phi(\vartheta(t,x,y),m(t,x,y)) \, Q(x,y)\,\pi(x) \d t\\
& \text{with} \quad
\Phi(\vartheta,m) := \begin{cases}
\frac{m^2}{\vartheta} & \text{if } \vartheta>0, \\
0 & \text{if } (m,\vartheta) = (0,0), \nonumber \\
+\infty & \text{else.}
\end{cases}
\end{align}
\end{Lemma}
Note that $\Phi$ is the integrand of the Benamou--Brenier action functional and that $\alpha(s,t,m) = \Phi(\theta(s,t),m)$.
\begin{proof}
The first equality is merely the definition of $\mathcal{A}_h$ and using the fact that $\avg_h m = m$ for $m \in V_{e,h}^0$.
For the second equality note that for any $\vartheta \in V_{e,h}^0$ with $(\bar \rho,\vartheta) \in \mathcal{K}_{\tn{pre}}$ one has $\vartheta(t_i,x,y) \leq \theta(\bar \rho(t_i,x),\bar \rho(t_i,y))$. By monotonicity of $\Phi$ in its first argument this implies 
$\Phi(\vartheta(t_i,x,y),m(t_i,x,y)) \geq \alpha(\bar \rho(t_i,x),\bar \rho(t_i,y),m(t_i,x,y))$ and hence
\begin{align}\label{eq:Aest}
\mathcal{A}(\bar \rho,m) \leq \inf \left\{\widehat{\mathcal{A}}(\vartheta,m) + {\mathcal{I}}_{\mathcal{K}_{\tn{pre}}}(\bar \rho,\vartheta) : \vartheta \in \mathbb{R}^{\mathcal{X} \times \mathcal{X}} \right\}.
\end{align}
Further, we obviously have that
$\bar\vartheta(t_i,x,y):=\theta(\bar \rho(t_i,x),\bar \rho(t_i,y))$
satisfies $(\bar \rho,\vartheta)\in\mathcal{K}_{\tn{pre}}$ and
$\widehat{\mathcal{A}}(\bar\vartheta,m)=\mathcal{A}(\bar \rho,m)$. Hence,
we have equality in \eqref{eq:Aest}.
\end{proof}

The proximal mapping of the function $\widehat{\mathcal{A}}$ can be computed separately for each time interval and graph edge. 
However, the set $\mathcal{K}_{\tn{pre}}$ still couples the variables $\avg_h \rho$ and $\vartheta$ according to the graph structure and the averaging operator $\avg_h$ couples the variables of $\rho$ in time.
To resolve this, we introduce a second set of auxiliary variables.
\begin{Lemma}
For $\bar \rho \in \mathbb{R}^\mathcal{X}$, $\vartheta \in \mathbb{R}^{\mathcal{X} \times \mathcal{X}}$ one finds
\begin{align} \nonumber
{\mathcal{I}}_{\mathcal{K}_{\tn{pre}}}(\avg_h \rho,\vartheta)  = \inf \Big\{ &
{\mathcal{I}}_{{\mathcal{J}_{avg}}}(\rho,\bar \rho)
+ {\mathcal{I}}_{{\mathcal{J}_{=}}}(\bar \rho,q)
+ {\mathcal{I}}_{{\mathcal{J}_{\pm}}}(q,\massNodes^{-},\massNodes^{+})
+ {\mathcal{I}}_{\mathcal{K}}(\massNodes^{-},\massNodes^{+},\vartheta)
\, : \,\\
& (\bar \rho,q,\massNodes^{-}, \massNodes^{+}) \in
(V_{n,h}^0)^2 \times (V_{e,h}^0)^2
\Big\} \label{eq:actionReform2}
\end{align}
\begin{align} \intertext{where}
{\mathcal{J}_{avg}} & := \left\{
(\rho, \bar \rho) \in V_{n,h}^{1} \times V_{n,h}^0
\, \colon \,
\bar \rho = \avg_h \rho
     \right	\}, \\
{\mathcal{J}_{=}} & := \left\{(\bar \rho,q) \in (V_{n,h}^0)^2
\,\colon\, \bar \rho=q \right\}, \\
\label{def:SetJ}
{\mathcal{J}_{\pm}} & := \left\{ (q, \massNodes^{-}, \massNodes^{+})
\in V_{n,h}^0 \times (V_{e,h}^0)^2
\,:\,
q(t_i,x) = \massNodes^{-}(t_i,x,y), \,
q(t_i,y) = \massNodes^{+}(t_i,x,y)
\right\},\\
\label{def:SetK}
\mathcal{K} & := \left\{ (\massNodes^{-}, \massNodes^{+}, \vartheta) \in (V_{e,h}^0)^3
\,:\,
(\massNodes^{-}(t_i,x,y), \massNodes^{+}(t_i,x,y), \vartheta(t_i,x,y)) \in K \right\}, \\
\intertext{with}
K & := \{ (\massNodes^{-}, \massNodes^{+}, \vartheta) \in \mathbb{R}^3 \,:\,
0 \leq \vartheta \leq \theta(\massNodes^{-},\massNodes^{+}) \}.
\end{align}
\end{Lemma}
\begin{proof}
For fixed $\rho \in V_{n,h}^{1}$ there is precisely one tuple $(\bar \rho,q,\massNodes^{-},\massNodes^{+})$ such that 
\[(\rho,\bar \rho) \in {\mathcal{J}_{avg}}, \quad (\bar \rho,q) \in {\mathcal{J}_{=}},\; \text{and} \quad (q,\massNodes^{-},\massNodes^{+}) \in {\mathcal{J}_{\pm}}\,,\] 
given by $\bar \rho=\avg_h \rho$, $q=\bar \rho$, $\massNodes^{-}(t_i,x,y)=q(t_i,x)$, $\massNodes^{+}(t_i,x,y) = q(t_i,y)$. For this $(\massNodes^{-},\massNodes^{+})$ one finds $(\massNodes^{-}, \massNodes^{+}, \vartheta) \in \mathcal{K}$ if and only if $(\avg_h \rho,\vartheta) \in \mathcal{K}_{\tn{pre}}$.
\end{proof}

The function ${\mathcal{I}}_{{\mathcal{J}_{avg}}}$ relates the values of $\rho$ on time nodes to the average values on the adjacent time intervals, ${\mathcal{I}}_{{\mathcal{J}_{\pm}}}$ communicates the values of $q$ on graph nodes to the adjacent graph edges and ${\mathcal{I}}_\mathcal{K}$ ensures the mass averaging via the function $\theta$. The additional splitting via ${\mathcal{I}}_{{\mathcal{J}_{=}}}$ will later simplify partition of the final optimization problem into primal and dual component.
The sets ${\mathcal{J}_{avg}}$, ${\mathcal{J}_{=}}$, ${\mathcal{J}_{\pm}}$ and $\mathcal{K}$ are all products of simpler low-dimensional sets, implying simpler computation of the relevant proximal mappings and projections.

This gives us an equivalent formulation for the discrete minimization problem \eqref{def:metricDiscrete}:
\begin{align} \nonumber
\mathcal{W}_h(\bar \rho_A, \bar \rho_B)^2
=  \inf
\Big\{ & (\mathcal{F} + \mathcal{G})(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h)
\,:\,\\
&(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h) \in V_{n,h}^{1} \times (V_{e,h}^0)^4 \times (V_{n,h}^0)^2
\Big\} 
\label{def:metricWithUAndV}
\end{align}
with
\begin{align*}
\mathcal{F}(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h)
& :=
\widehat{\mathcal{A}}(\vartheta_h, m_h)
+ {\mathcal{I}}_{{\mathcal{J}_{\pm}}}(q_h,\massNodes^{-}_h, \massNodes^{+}_h) 
+ {\mathcal{I}}_{{\mathcal{J}_{avg}}}(\rho_h,\bar \rho_h), \\
\mathcal{G}(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h)
& :=
{\mathcal{I}}_{\mathcal{CE}_h(\rho_A, \rho_B)}(\rho_h, m_h) 
+ {\mathcal{I}}_{\mathcal{K}}(\massNodes^{-}_h,\massNodes^{+}_h,\vartheta_h)
+ {\mathcal{I}}_{{\mathcal{J}_{=}}}(\bar \rho_h,q_h).
\end{align*}

The structure of this optimization problem is well suited for the first order primal-dual algorithm presented in \cite{ChPo11}. 
We consider the Hilbert space $H = V_{n,h}^{1} \times (V_{e,h}^0)^4 \times (V_{n,h}^0)^2$ composed of tuples of functions in space and time
with the scalar product
\begin{align}
& \InProdH{\left( \rho_{h,1},m_{h,1},\vartheta_{h,1},\massNodes^{-}_{h,1},\massNodes^{+}_{h,1},\bar \rho_{h,1}, q_{h,1}  \right)}
          {\left( \rho_{h,2},m_{h,2},\vartheta_{h,2},\massNodes^{-}_{h,2},\massNodes^{+}_{h,2},\bar \rho_{h,2}, q_{h,2} \right)} \nonumber \\
& \quad := h \sum_{i=0}^N \InProdV{\rho_{h,1}(t_i,\cdot)}{\rho_{h,2}(t_i,\cdot)}
   + h \sum_{i=0}^{N-1}         
         \InProdV{\bar \rho_{h,1}(t_i,\cdot)}{\bar \rho_{h,2}(t_i,\cdot)}
         +\InProdV{q_{h,1}(t_i,\cdot)}{q_{h,2}(t_i,\cdot)} \nonumber \\
& \qquad + h \sum_{i=0}^{N-1} 
          \InProdE{m_{h,1}(t_i,\cdot)}{m_{h,2}(t_i,\cdot)} 
        +\InProdE{\vartheta_{h,1}(t_i,\cdot)}{\vartheta_{h,2}(t_i,\cdot)} \nonumber \\
  & \qquad + h \sum_{i=0}^{N-1}        \InProdE{\massNodes^{-}_{h,1}(t_i,\cdot)}{\massNodes^{-}_{h,2}(t_i,\cdot)}
         +\InProdE{\massNodes^{+}_{h,1}(t_i,\cdot)}{\massNodes^{+}_{h,2}(t_i,\cdot)} \,.
\label{eq:ScalarProductDiscrete}
\end{align}
and the induced norm denoted by $\NormH{\cdot}$.
Then applying \cite[Algorithm 1]{ChPo11} to solve problem \eqref{def:metricWithUAndV} with $\mathcal{F}, \,\mathcal{G}: \,H \to \R \cup \{\infty\}$
amounts to iteratively compute for initial data $(\iterz{a},\iterz{b}) \in H^2$ and $\iterz{\widebar{a}}=\iterz{a}$  
\begin{align}
\iterll{b} & = \mathrm{prox}_{\sigma\,\mathcal{F}^\ast}(\iterl{b}+\sigma\,\iterl{\widebar{a}}),\nonumber  \\
\iterll{a} & = \mathrm{prox}_{\tau\,\mathcal{G}}(\iterl{a}-\tau\,\iterll{b}), \label{eq:ChambollePock} \\
\iterll{\widebar{a}} & = \iterll{a} + \lambda \cdot (\iterll{a}-\iterl{a}).\nonumber
\end{align}
where $\tau$, $\sigma>0$, $\lambda \in [0,1]$ . 
As demonstrated in \cite{ChPo11}  the iterates converge to a minimizer in \eqref{def:metricWithUAndV} if $\tau \cdot \sigma < 1$.
For some $(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h) \in H$ one finds
\begin{align*}
 \mathcal{F}^\ast(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h) 
  = \widehat{\mathcal{A}}^\ast(\vartheta_h,m_h)
  + {\mathcal{I}}_{{\mathcal{J}_{\pm}}}^\ast(q_h,\massNodes^{-}_h, \massNodes^{+}_h)
  + {\mathcal{I}}_{{\mathcal{J}_{avg}}}^\ast(\rho_h, \bar \rho_h)
\end{align*}
and the proximal mapping 
$(\rho^{\tn{pr}}_h,m^{\tn{pr}}_h,\vartheta^{\tn{pr}}_h,{\massNodes^{-}_h}^{\tn{pr}},{\massNodes^{+}_h}^{\tn{pr}}, \bar \rho^{\tn{pr}}_h, q^{\tn{pr}}_h) 
= \mathrm{prox}_{\sigma\,\mathcal{F}^\ast}(\rho_h, m_h, \vartheta_h, \massNodes^{-}_h, \massNodes^{+}_h, \bar \rho_h, q_h)$ decomposes as follows:
\begin{align*}
 (\vartheta^{\tn{pr}}_h,m^{\tn{pr}}_h)
 & = \mathrm{prox}_{\sigma\,\widehat{\mathcal{A}}^\ast}(\vartheta_h,m_h)\,, \\
 (q^{\tn{pr}}_h,{\massNodes^{-}_h}^{\tn{pr}},{\massNodes^{-}_h}^{\tn{pr}}) 
 & = \mathrm{prox}_{\sigma\,{\mathcal{I}}_{{\mathcal{J}_{\pm}}}^\ast}(q_h,\massNodes^{-}_h,\massNodes^{+}_h) \, , \\
 (\rho^{\tn{pr}}_h, \bar \rho^{\tn{pr}}_h) 
 & = \mathrm{prox}_{\sigma \, {\mathcal{I}}_{{\mathcal{J}_{avg}}}^\ast}(\rho_h,\bar \rho_h) \, .
\end{align*}
Likewise, for 
$(\rho^{\tn{pr}}_h,m^{\tn{pr}}_h,\vartheta^{\tn{pr}}_h,{\massNodes^{-}_h}^{\tn{pr}},{\massNodes^{+}_h}^{\tn{pr}},\bar \rho^{\tn{pr}}_h, q^{\tn{pr}}_h) 
= \mathrm{prox}_{\tau\,\mathcal{G}}(\rho_h,\allowbreak m_h,\allowbreak \vartheta_h,\allowbreak \massNodes^{-}_h,\allowbreak \massNodes^{+}_h, \bar \rho_h, q_h)$ one finds
\begin{align*}
 (\rho^{\tn{pr}}_h,m^{\tn{pr}}_h) & =
 \mathrm{proj}_{\mathcal{CE}_h(\rho_A,\rho_B)}(\rho_h,m_h)\,, \\
 ({\massNodes^{-}}^{\tn{pr}}_h,{\massNodes^{+}}^{\tn{pr}}_h,\vartheta_h) 
 & = \mathrm{proj}_{\mathcal{K}}(\massNodes^{-}_h,\massNodes^{+}_h,\vartheta_h)\, , \\
 (\bar \rho^{\tn{pr}}_h, q^{\tn{pr}}_h) 
 &= \mathrm{proj}_{{\mathcal{J}_{=}}}(\bar \rho_h, q_h) \, .
\end{align*}
Each of the proximal maps is performed with respect to the norm $\NormH{\cdot}$ restricted to the relevant variables.
\bigskip

In what follows, we will study these maps in detail.
In fact, we will observe that 
$\mathrm{prox}_{\widehat{\mathcal{A}}^\ast}$ and $\mathrm{proj}_{\mathcal{K}}$ can be separated into low-dimensional problems over each time-step and edge $(x,y) \in \mathcal{X} \times \mathcal{X}$,
$\mathrm{prox}_{{\mathcal{I}}_{{\mathcal{J}_{\pm}}}^\ast}$ splits into low-dimensional problems for each time-step and node $x \in \mathcal{X}$,
$\mathrm{prox}_{{\mathcal{I}}_{{\mathcal{J}_{=}}}}$ is a simple pointwise update,
$\mathrm{prox}_{{\mathcal{I}}_{{\mathcal{J}_{avg}}}^\ast}$ decouples for each node $x \in \mathcal{X}$ to a sparse linear system in time,
and $\mathrm{proj}_{\mathcal{CE}_h(\rho_A,\rho_B)}$ can be computed solving a linear system, which is sparse if $Q$ is sparse.
Consequently, $\mathrm{prox}_{\mathcal{F}^\ast}$ and $\mathrm{prox}_\mathcal{G}$ can be computed efficiently and ensure that the above scheme is well-suited to solve \eqref{def:metricWithUAndV}.




\subsection[Projection onto CE]{Projection onto $\mathcal{CE}_h(\rho_A,\rho_B)$}
For given $(\rho_h,m_h) \in V_{n,h}^{1} \times V_{e,h}^0$ we need to solve the following problem:
\begin{align} \label{eq:ProjCE}
 \mathrm{proj}_{\mathcal{CE}_h(\rho_A,\rho_B)}(\rho_h,m_h) 
 = \argmin_{(\rho^{\tn{pr}}_h,m^{\tn{pr}}_h) \in  \mathcal{CE}_h(\rho_A,\rho_B)}
 \frac{h}{2} \sum_{i=0}^N \NormV{\rho^{\tn{pr}}_h(t_i,\cdot)-\rho_h(t_i,\cdot)}^2 + \frac{h}{2} \sum_{i=0}^{N-1} \NormE{m^{\tn{pr}}_h(t_i,\cdot)-m_h(t_i,\cdot)}^2
\end{align}
To this end we take into account the following dual formulation.
\begin{Proposition}\label{prop:ProjCEDual}
The solution $(\rho^{\tn{pr}}_h,m^{\tn{pr}}_h)$ to \eqref{eq:ProjCE} is given by
\begin{subequations}
\label{eq:ProjCEPDRelationGF}
\begin{align}
 \label{eq:ProjCEPDRelationGF:A}
 \rho^{\tn{pr}}_h(t_i,x) & = \rho_h(t_i,x) + \frac{\varphi_h(t_i,x) - \varphi_h(t_{i-1},x)}{h} \,, \quad \forall i=1,\ldots,N-1 \, , \\
 \label{eq:ProjCEPDRelationGF:B}
 \rho^{\tn{pr}}_h(t_0,x) & = \rho_A(x) \, , \; \rho^{\tn{pr}}_h(t_N,x) = \rho_B(x) \, , \\
 \label{eq:ProjCEPDRelationGF:C}
 m^{\tn{pr}}_h(t_i,x,y) &= m_h(t_i,x,y) + \nabla_\mathcal{X} \varphi_h(t_i,x,y)  \,, \quad \forall i=1,\ldots,N-1 \, .
\end{align}
\end{subequations}
where $\varphi_h$  solves the space time elliptic equation

\begin{align}
& \pi(x) \frac{ \varphi_h(t_1,x) - \varphi_h(t_{0},x) }{h^2} + \pi(x)\,\triangle_\mathcal{X} \varphi_h(t_0,x)
  = - \pi(x) \left( \frac{\rho_h(t_{1},x) - \rho_A(x)}{h} + \div m_h(t_0,x) \right) \; , \nonumber \\
& \pi(x) \frac{ -\varphi_h(t_{N-1},x) + \varphi_h(t_{N-2},x)}{h^2} + \pi(x) \triangle_\mathcal{X} \varphi_h(t_{N-1},x)  \nonumber \\
& \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad 
  = - \pi(x) \left( \frac{\rho_B(x) - \rho_h(t_{N-1},x)}{h} + \div m_h(t_{N-1},x) \right) \nonumber \\
 &\pi(x) \frac{\varphi_h(t_{i+1},x) - 2 \varphi_h(t_i,x) + \varphi_h(t_{i-1},x) }{h^2} + 
 \pi(x)\,\triangle_\mathcal{X} \varphi_h(t_i,x)   \nonumber\\
& \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad = - \pi(x) \left( \frac{\rho_h(t_{i+1},x) - \rho_h(t_i,x)}{h} + \div m_h(t_i,x) \right) \label{eq:ProjCELM}
\end{align}
for $i=1,\ldots,N-2$ and $x \in \mathcal{X}$.
\end{Proposition}
The factors $\pi(x)$ in \eqref{eq:ProjCELM} could be canceled but they will simplify further analysis.
\begin{proof}
We define the Lagrangian corresponding to \eqref{eq:ProjCE} as
\begin{align*}
 L[\rho_h^{\tn{pr}}, m_h^{\tn{pr}},\varphi_h,\lambda_A,\lambda_B]
 =&  \frac{h}{2} \sum_{i=0}^N \NormV{\rho^{\tn{pr}}_h(t_i,\cdot)-\rho_h(t_i,\cdot)}^2 + \frac{h}{2} \sum_{i=0}^{N-1} \NormE{m^{\tn{pr}}_h(t_i,\cdot)-m_h(t_i,\cdot)}^2 \\
 & +h \sum_{i=0}^{N-1} \sum_{x \in \mathcal{X}} \varphi_h(t_i,x) \left( \frac{\rho_h^{\tn{pr}}(t_{i+1},x) - \rho_h^{\tn{pr}}(t_i,x)}{h} + \div_\mathcal{X} m_h^{\tn{pr}}(t_i,x) \right) \pi(x) \\
 &+  \sum_{x \in \mathcal{X}} \left( \lambda_B(x) ( \rho_h(t_N,x) - \rho_B(x) ) + \lambda_A(x) ( \rho_h(t_0,x) - \rho_A(x) ) \right) \pi(x)
\end{align*}
where $\lambda_A, \lambda_B$ are the Lagrange multipliers for the boundary conditions $\rho_h(t_0,\cdot) = \rho_A$, $\rho_h(t_N,\cdot) = \rho_B$.
The optimality condition in $\rho_h$ and $m_h$ directly imply \eqref{eq:ProjCEPDRelationGF:A} and \eqref{eq:ProjCEPDRelationGF:C}. \eqref{eq:ProjCEPDRelationGF:B} reflects the boundary conditions, which are to be ensured in $\mathcal{CE}_h(\rho_A,\rho_B)$.
Inserting these relations into the continuity equation $\partial_t \rho_h^{\tn{pr}} + \div m_h^{\tn{pr}}=0$ leads to the system of equations \eqref{eq:ProjCELM}.
\end{proof}

The Lagrange multiplier $\varphi_h$ in Proposition \ref{prop:ProjCEDual} lives in $V_{n,h}^0$ which can be identified with $\mathbb{R}^{N  \mathrm{card} \mathcal{X}}$. We equip this space with the canonical basis
\begin{align*}
(\varphi_h^{i,x})_{i=0,\ldots, N-1,\,x\in \mathcal{X}}
\quad \tn{where} \quad
(\varphi_h^{i,x})(t_j,y) = \delta_{i,j} \cdot \delta_{x,y}
\end{align*}
and the standard Euclidean inner product with respect to this basis.
Then the elliptic equation \eqref{eq:ProjCELM} can be written as a linear system $S Z = F$ 
for a coordinate vector $Z = (\varphi_h(t_i,x))_{i=0,\ldots N-1,\, x\in X}$,
a matrix $S \in \mathbb{R}^{(N\mathrm{card} \mathcal{X}) \times (N \mathrm{card} \mathcal{X})}$
and a vector $F \in \mathbb{R}^{N  \mathrm{card} \mathcal{X}}$.
The matrix $S$ is symmetric since $\pi(x)Q(x,y) = \pi(y) Q(y,x)$ and the matrix representation of $\triangle_\mathcal{X}$ is $Q - \text{diag}(\sum_yQ(\cdot,y))$.
Furthermore, $S$ is sparse if $Q$ is sparse.
However, the matrix $S$ is not invertible, its kernel is spanned by functions that are constant in space and time. 
To see this, assume that a non constant $Z$  
is in the kernel of $S$ and denote by $\phi_h$ the associated function in $V_{n,h}^0$. Now, 
let  $I_+(\mu) := \{(i,x) \in \{0,\ldots, N-1\} \times \mathcal{X}\,:\,\phi_h(i,x)>\mu\}$ for $\mu = \min \phi_h(i,x)$
and define $\psi_h\in V_{n,h}^0$ via $\psi(t_i, x) = 1$ if $(i,x)\in I_+(\mu)$ and $\psi_h(t_i, x) = 0$ else. 
Let $W$ be the associated nodal vector to $\psi$. 
By assumption on $Z$ the set $I_+(\mu)$ is non empty and thus it is easy to see that $W^\top S Z <0$ and thus $Z$ 
can not be in the kernel of $S$, which proves the claim.

We impose the additional constraint $\sum_{i=0}^{N-1} \sum_{x \in \mathcal{X}} \varphi_h(t_i,x) = 0$ to remove this ambiguity. 
This can be written as $w^\top \varphi_h=0$ where $w \in \mathbb{R}^{N  \mathrm{card} \mathcal{X}}$ is the vector with entries $w^{i,x} = 1$ leading to the linear system 
\begin{align*}
 \begin{pmatrix}
  S & w \\
  w^T & 0 
 \end{pmatrix}
 \begin{pmatrix}
  Z\\
  \lambda
 \end{pmatrix}
 =
 \begin{pmatrix}
  F \\
  0
 \end{pmatrix}
 \, .
\end{align*}
This system is uniquely solvable and the solution implies $\lambda=0$ if $F \perp w$ (in the Euclidean sense), which is true because $\rho_A$ and $\rho_B$ are assumed to be of equal mass.

\subsection[{Proximal Mapping of A*}]{Proximal Mapping of $\widehat{\mathcal{A}}^\ast$}
\newcommand{\Bc}{\mathcal{B}}
The function $\widehat{\mathcal{A}}$ is convex and 1-homogeneous, hence its Fenchel conjugate is the indicator function of a convex set and the proximal mapping of $\widehat{\mathcal{A}}^\ast$ is a projection.
For $(\vartheta,m) \in (V_{e,h}^0)^2$ one has
\begin{align*}
\widehat{\mathcal{A}}(\vartheta,m) & = h \sum_{i=0}^{N-1} \sum_{x,y \in \mathcal{X}}
\Phi(\vartheta(t_i,x,y),m(t_i,x,y)) \, Q(x,y)\,\pi(x)\,.
\end{align*}
Following \cite{BeBr00} a direct calculation for $(p,q) \in (V_{e,h}^0)^2$ yields
\begin{align*}
\widehat{\mathcal{A}}^\ast(p,q) & = \sup_{(\vartheta,m) \in (V_{e,h}^0)^2}
h \sum_{i=0}^{N-1}
\Big[
\InProdE{p(t_i,\cdot,\cdot)}{\vartheta(t_i,\cdot,\cdot)}
+ \InProdE{q(t_i,\cdot,\cdot)}{m(t_i,\cdot,\cdot)} \\[-1ex]
& \qquad \qquad \qquad \qquad - \frac12 \sum_{(x,y) \in \mathcal{X} \times \mathcal{X}} \Phi(\vartheta(t_i,x,y),m(t_i,x,y))\,Q(x,y)\,\pi(x) \Big] \\
& = \frac{h}{2} \sum_{ \genfrac{}{}{0pt}{}{i=0,\ldots, N-1 }{(x,y) \in \mathcal{X} \times \mathcal{X}} }
\Phi^\ast(p(t_i,x,y),q(t_i,x,y))\,Q(x,y)\,\pi(x)
=\sum_{ \genfrac{}{}{0pt}{}{i=0,\ldots, N-1 }{(x,y) \in \mathcal{X} \times \mathcal{X}} }
{\mathcal{I}}_\Bc(p(t_i,x,y),q(t_i,x,y))
\end{align*}
with $\Phi^\ast={\mathcal{I}}_\Bc$ for $\Bc = \{ (p,q) \in \mathbb{R}^2\,:\, p+\frac{q^2}{4} \leq 0 \}\,$.
Thus the proximal mapping separates into two-dimensional problems for each time interval and graph edge and $(p^{{\tn{pr}}},q^{{\tn{pr}}}) = \mathrm{prox}_{\sigma\,\widehat{\mathcal{A}}^\ast}(p,q)$ precisely if
\begin{align*}
(p^{{\tn{pr}}}(t_i,x,y),q^{{\tn{pr}}}(t_i,x,y)) = \mathrm{proj}_{\Bc}(p(t_i,x,y),q(t_i,x,y))\,,
\end{align*}
where $\mathrm{proj}_{\Bc}$ is the projection with respect to the standard Euclidean distance on $\mathbb{R}^2$ and  a Newton scheme in $\mathbb{R}$ can be used 
to solve for this projection.
Since this proximal mapping is a projection, it is in particular independent of the step size $\sigma$.


\subsection[Projection onto K]{Projection onto $\mathcal{K}$}
For given $(\massNodes^{-},\massNodes^{+},\vartheta) \in (V_{e,h}^0)^3$ we need to solve
\begin{align*}
\mathrm{proj}_{\mathcal{K}}(\massNodes^{-},\massNodes^{+},\vartheta) =\!\! \argmin_{
({\massNodes^{-}}^{\tn{pr}},{\massNodes^{+}}^{\tn{pr}},\vartheta^{\tn{pr}}) \in \mathcal{K}}
\frac{h}{2} \sum_{i=0}^{N-1} \Big( &
\NormE{{\massNodes^{-}}^{\tn{pr}}(t_i,\cdot,\cdot)-\massNodes^{-}(t_i,\cdot,\cdot)}^2
+ \NormE{{\massNodes^{+}}^{\tn{pr}}(t_i,\cdot,\cdot)-\massNodes^{+}(t_i,\cdot,\cdot)}^2 \\[-1ex]
& + \NormE{\vartheta^{\tn{pr}}(t_i,\cdot,\cdot)-\vartheta(t_i,\cdot,\cdot)}^2
\Big)\,.
\end{align*}
Recall that $\mathcal{K}$ is a product of the tree-dimensional closed convex set $K$, as indicated in \eqref{def:SetK}. Therefore 
$({\massNodes^{-}}^{\tn{pr}},{\massNodes^{+}}^{\tn{pr}},\vartheta^{\tn{pr}}) = \mathrm{proj}_{\mathcal{K}}(\massNodes^{-},\massNodes^{+},\vartheta)$
decouples into the edgewise projection in each time step, i.e.
\begin{align*}
({\massNodes^{-}}^{\tn{pr}}(t_i,x,y),{\massNodes^{+}}^{\tn{pr}}(t_i,x,y),\vartheta^{\tn{pr}}(t_i,x,y)) = \mathrm{proj}_{K}(\massNodes^{-}(t_i,x,y),\massNodes^{+}(t_i,x,y),\vartheta(t_i,x,y))
\end{align*}
where this projection is with respect to the standard Euclidean distance on $\mathbb{R}^3$.
Let us denote by $\partial^+ \theta(x)$ the super-differential of $\theta$ at $x \in \mathbb{R}^2$, which is the analogue of the sub-differential for concave functions.
More precisely, $\partial^+ \theta(x) = -\partial (-\theta)(x)$, where $\partial (-\theta)(x)$ is the sub-differential of the convex function $x \mapsto -\theta(x)$ at $x$.
Then the projection $p^{\pr} = \mathrm{proj}_{K}(p)$ of  $p \in \mathbb{R}^3$ is characterized by \cite[Prop.~ 6.46]{ConvexFunctionalAnalysis-11}
\begin{align}
\label{def:nCone}
p - p^{\pr} & \in N_K(p^{\pr}) :=
\{ z \in \mathbb{R}^3 \,:\, \langle z, q-p^{\pr} \rangle \leq 0 \, \forall \, q \in K \}\,,
\end{align}
where $N_K(p^{\pr})$ is the normal cone of $K$ at $p^{\pr}$.
To solve this inclusion we distinguish the following cases:
\begin{Lemma}
\label{lem:ProjKNCones}
For an averaging function $\theta:\mathbb{R}^2 \to \mathbb{R}$ fulfilling the assumptions listed in Section \ref{sec:intro} 
and for $K := \{ p \in \mathbb{R}^3\, : \, 0 \leq p_3 \leq \theta(p_1,p_2) \}$
the normal cone $N_K(p^{\pr})$ for $p^{\pr} \in K$ is given by:
\begin{enumerate}[label=(\roman*)]
\item 
\label{item:LemmaNConeTrivial}
Trivial projection: $p = p^{\pr} \in \inter K \!= \! \{ (p_1,p_2,p_3) \in \mathbb{R}^3 : 0 \! <\!  p_3 \! < \! \theta(p_1,p_2)\}$, then $N_K(p^{\pr}) = \{0\}$.
\item 
\label{item:LemmaNConeBottom}
Projection onto `bottom facet' of $K$: $p^{\pr} \in (0,+\infty) \times (0,+\infty) \times \{0\}$, then
$N_K(p^{\pr}) = \{0\} \times \{0\} \times \mathbb{R}^-_0$.
\item 
\label{item:LemmaNConeAxis}
Projection onto coordinate axis: $p^{\pr} = (p^{\pr}_1,0,0)$ for $p^{\pr}_1 \in (0,+\infty)$, then 
\begin{align*}
N_K(p^{\pr})= \{0\} \times \mathbb{R}^-_0 \times \mathbb{R}^-_0
\cup \left\{ (0,q_2,q_3) \in \{0\} \times \mathbb{R}^-_0 \times (0,+\infty) \,:\, 
(0,-q_2/q_3) \in \partial^+ \theta(p^{\pr}_1,0) \right\}.
\end{align*}
Note that $(0,q) \in \partial^+ \theta(p^{\pr}_1,0)$ is equivalent to $q \geq \lim_{z \searrow 0} \partial_2 \theta(p^{\pr}_1,z)$ and that $\partial^+ \theta(p^{\pr}_1,0)$ is empty if $\lim_{z \searrow 0} \partial_2 \theta(p^{\pr}_1,z)=\infty$.
The analogous representation holds for the second axis.
\item \label{item:LemmaNConeOrigin}
Projection onto origin: $p^{\pr} = (0,0,0)$, then
\[
N_K(p^{\pr})= (\mathbb{R}^-_0)^3 \cup \{ (q_1,q_2,q_3) \in  \mathbb{R}^-_0 \times \mathbb{R}^-_0 \times (0,+\infty) \,:\,
(q_1/q_3,q_2/q_3) \in -\partial^+ \theta(0) \}\,.
\]
\item 
\label{item:LemmaNConeTop}
Projection onto `upper surface' of $K$: $p^{\pr} = (p^{\pr}_1,p^{\pr}_2,\theta(p^{\pr}_1,p^{\pr}_2))$ for $(p^{\pr}_1,p^{\pr}_2) \in (0,+\infty)^2$, then
\[
N_K(p^{\pr})=\{ \lambda \cdot (-\partial_1 \theta(p^{\pr}_1,p^{\pr}_2),-\partial_2 \theta(p^{\pr}_1,p^{\pr}_2),1) \,:\, \lambda \in \mathbb{R}^+\}\,.
\]
\end{enumerate}
\end{Lemma}
\begin{proof}
For $p^{\pr} \in \inter K$ one finds $N_K(p^{\pr})=\{0\}$ and thus $p^{\pr}=p$, which implies \ref{item:LemmaNConeTrivial}.

In case \ref{item:LemmaNConeBottom} the set $\mathbb{R} \times \mathbb{R} \times \{0\}$ is obviously the only supporting plane of $K$ that contains $p^{\pr}$. Thus the normal cone is just the ray
in direction $(0,0,-1)$.

Assume $p^{\pr}=(p^{\pr}_1,0,0)$, $p^{\pr}_1>0$. Then there is some $\varepsilon>0$ such that $\{(p^{\pr}_1+\varepsilon,0,0),(p^{\pr}_1-\varepsilon,0,0),(p^{\pr}_1,\varepsilon,0)\} \subset K$.
Therefore $N_K(p^{\pr}) \subset \{0\} \times \mathbb{R}^-_0 \times \mathbb{R}$.
Since $\mathbb{R} \times \{0\} \times \mathbb{R}$ and $\mathbb{R} \times \mathbb{R} \times \{0\}$ are supporting planes of $K$ that contain $p^{\pr}$, one must have $\{0\} \times \mathbb{R}^-_0 \times \mathbb{R}^-_0 \subset N_K(p^{\pr})$.
Moreover, for $\lim_{z \searrow 0} \partial_2 \theta(p^{\pr}_1,z)<\infty$ 
let $z=(z_1,z_2) \in \partial^+ \theta(p^{\pr}_1,0)$. One must have $z_1=0$ and $z_2 \in \partial^+ f(0)$ with auxiliary function $f : t \mapsto \theta(p^{\pr}_1,t)$.
Then $\{ q \in \mathbb{R}^3 \,:\, \langle q-p^{\pr},(0,-z_2,1)\rangle = 0\}$ is a supporting plane of $K$ and consequently $(0,-z_2,1) \in N_K(p^{\pr})$. Conversely, from $z_2 \notin \partial^+ f(0)$ follows $(0,-z_2,1) \notin N_K(p^{\pr})$.
So \[N_K(p^{\pr}) = \{0\} \times \mathbb{R}^-_0 \times \mathbb{R}^-_0 \cup \{(0,-\lambda \cdot z, \lambda)\,:\, z \in \partial^+ f(0), \lambda \in (0,+\infty) \}.\]
The auxiliary function $f$ is concave and by monotonicity of the super-differential we find $\partial^+ f(0) = [\lim_{z \searrow 0} \partial_2 \theta(p^{\pr}_1,z), +\infty)$. 
With this characterization we arrive at the expression for $N_K(p^{\pr})$ as given in \ref{item:LemmaNConeAxis}. 
The proof for the second axis is analogous.

For $p^{\pr} = (0,0,0)$ we find $(\mathbb{R}^-_0)^3 \subset N_K(0) \subset \mathbb{R}^-_0 \times \mathbb{R}^-_0 \times \mathbb{R}$ 
with arguments analogous to those in case \ref{item:LemmaNConeAxis}. 
For every $z =(z_1,z_2) \in \partial^+ \theta(0)$ a supporting plane through $0$ is given by $\{q \in \mathbb{R}^3 : \langle q, (-z_1,-z_2,1) \rangle = 0 \}$ and hence $(-z_1,-z_2,1) \in N_K(0)$. Conversely, $z=(z_1,z_2) \notin \partial^+ \theta(0)$ implies $(-z_1,-z_2,1) \notin N_K(0)$. With this, one obtains the expression for $N_K(0)$ given in \ref{item:LemmaNConeOrigin}.

Finally, we consider $p^{\pr} = (p^{\pr}_1,p^{\pr}_2,\theta(p^{\pr}_1,p^{\pr}_2))$ with $(p^{\pr}_1,p^{\pr}_2) \in (0,+\infty)^2$. 
In a neighbourhood of $p^{\pr}$, $K$ is the subgraph of a concave, differentiable function. 
The unique supporting plane of $K$ through $p^{\pr}$ is given by 
$\{ q \in \mathbb{R}^3 : \langle q-p^{\pr}, (-\partial_1 \theta(p^{\pr}_1,p^{\pr}_2),-\partial_2 \theta(p^{\pr}_1,p^{\pr}_2),1) \rangle = 0\}$
and $(-\partial_1 \theta(p^{\pr}_1,p^{\pr}_2),-\partial_2 \theta(p^{\pr}_1,p^{\pr}_2),1)$ is the unique associated outer normal as stated in \ref{item:LemmaNConeTop}.
\end{proof}

Using Lemma \ref{lem:ProjKNCones} one can devise an algorithm for the projection onto $K$.
For $p=(p_1,p_2,p_3) \in \mathbb{R}^3$ the projection $p^{\pr}=\mathrm{proj}_{K}(p)$ can be determined as follows:
\bigskip

\hspace{10ex} \begin{minipage}{0.8\linewidth}
\begin{algorithmic}
\Function{Project$K$}{$p_1$,$p_2$,$p_3$}
\LineIf{$0 \leq p_3 \leq \theta(p_1,p_2)$} \Return $(p_1,p_2,p_3)$
\LineIf{$p_3 \leq 0$} \Return $(\max\{p_1,0\},\max\{p_2,0\},0)$
\If{$(p_1 > 0) \wedge (p_2 \leq 0)$} 
\LineIf{$-p_2/p_3 \geq \lim_{z \searrow 0} \partial_2 \theta(p_1,z)$}
\Return $(p_1,0,0)$
\EndIf
\If{$(p_1 \leq 0) \wedge (p_2 > 0)$}
\LineIf{$-p_1/p_3 \geq \lim_{z \searrow 0} \partial_1 \theta(z,p_2)$}
\Return $(0,p_2,0)$
\EndIf 
\If{$(p_1 \leq 0) \wedge (p_2 \leq 0)$}
\LineIf{$(-p_1/p_3,-p_2/p_3) \in \partial^+ \theta(0)$}
\Return $(0,0,0)$
\EndIf
\State \Return \Call{Project$K$Top}{$p_1$,$p_2$,$p_3$}
\EndFunction
\end{algorithmic}
\end{minipage}
\bigskip

The function \Call{Project$K$Top}{$p_1$,$p_2$,$p_3$}  in the above algorithm corresponds to case \ref{item:LemmaNConeTop} of Lemma \ref{lem:ProjKNCones}, where $p^{\pr}$ lies on the `upper surface' of $K$, defined by the graph surface of  $\theta$. It will be described in more detail below.
In the following we will occasionally use the curve $c: (0,\infty) \to \mathbb{R}^2; \; q \mapsto (q^{-1/2},q^{1/2})$ to parametrize orientations in $(0,\infty)^2$. Due to the 1-homogeneity of $\theta$, often it suffices to look at its values at $\theta(c(q))$. Alternative choices for $c$ are feasible as well.

\begin{Lemma}[Projection onto `upper surface' of $K$]\label{lem:SetKSmall}
Let $p \in \mathbb{R}^3$ with projection on $K$ given by $p^{\pr} = (p^{\pr}_1,p^{\pr}_2,\theta(p^{\pr}_1,p^{\pr}_2))$ with $(p^{\pr}_1,p^{\pr}_2) \in (0,+\infty)^2$.
Further, let $w(q) = (q^{1/2},q^{-1/2},\theta(q^{1/2},q^{-1/2}))$ be a parametrized curve on the `upper surface' 
and $n(q) = (-\partial_1 \theta(q^{1/2},q^{-1/2}), -\partial_2 \theta(q^{1/2},q^{-1/2}), 1)$ be the corresponding normal.
Then there exists a unique $(q,\tau) \in (0,\infty)^2$ s.t. $p^{\pr} = \tau w(q)$. We have that $q$ is the unique root of $q \mapsto \langle p  , w(q)  \times  n(q)\rangle$
and $\tau = \langle p , \frac{w(q)}{\|w(q)\|^2}\rangle$.
\end{Lemma}
\begin{proof}
Since $\theta$ is 1-homogeneous, any $p^{\pr}$ of the form $(p^{\pr}_1,p^{\pr}_2,\theta(p^{\pr}_1,p^{\pr}_2))$, $(p^{\pr}_1,p^{\pr}_2) \in (0,+\infty)^2$, can be written as $p^{\pr} = \tau \cdot w(q)$ for unique $q \in (0,+\infty)$ and $\tau \in (0,+\infty)$. In explicit, $q= p^{\pr}_1/p^{\pr}_2$ and $\tau=(p^{\pr}_1 \cdot p^{\pr}_2)^\frac12$. 
Now, $n(q)$ is orthogonal on the graph of $\theta$ and outward pointing.
Hence, $p$ lies in the plane spanned by $w(q)$ and $n(q)$.
This is equivalent to $\langle p, w(q) \times n(q) \rangle = 0$.
Since $p^{\pr}$ is unique, this must be the unique root of $q \mapsto \langle p, w(q) \times n(q) \rangle$.
Once $q$ is determined, we know the ray on which $p^{\pr}$ lies. To find $\tau$, one must solve the remaining one-dimensional projection onto the ray. 
Consequently, $\tau$ is the unique minimizer of $\tau \mapsto \frac{1}{2} \|p - \tau \cdot w(q)\|^2$, which concludes the proof.
\end{proof}
For case \ref{item:LemmaNConeOrigin} of Lemma \ref{lem:ProjKNCones}  we need to characterize the super-differential of $\theta$ at the origin.
\begin{Lemma} \label{lem:ProjKSupDiff}
The super-differential of $\theta$ at the origin is given by
\[\partial^+ \theta(0) = \ol{\{ \nabla \theta(q^{-1/2},q^{1/2}) \,:\, q \in (0,\infty) \} } +  (\mathbb{R}^+_0)^2\,.\]
\end{Lemma}
\begin{proof}
Due to the 1-homogeneity  of $\theta$
\begin{align*}
\langle \nabla \theta(\lambda p), \lambda r \rangle = \lim_{\epsilon \to 0} \frac{\theta(\lambda (p+\epsilon r))-\theta(\lambda p)}{\epsilon} 
= \lambda \lim_{\epsilon \to 0} \frac{\theta(p+\epsilon r)- \theta(p)}{\epsilon}  = \lambda\langle \nabla \theta( p), r \rangle
\end{align*}
for $p \in (0,+\infty)^2$, $\lambda >0$, and all $r\in \mathbb{R}^2$, which leads to $\nabla \theta(\lambda p) = \nabla \theta(p)$ for $p \in (0,+\infty)^2$ and $\lambda >0$.
Thus, for the curve $c: (0,\infty) \to \mathbb{R}^2; \; q \mapsto (q^{-1/2},q^{1/2})$ the set of tangent planes at $(c(q),\theta(c(q)))$ spanned by $(\nabla \theta(c(q)),1)$ and $(c(s),\theta(c(q)))$ for $q\in (0,\infty)$
is already the complete set of affine tangent planes to the graph of $\theta$ over $(0,\infty)^2$. Thus, by continuity of $\theta$ on $[0,\infty)^2$ 
we get $\theta(0) + \langle r,p\rangle \geq \theta(p)$ for $r \in \{ \nabla \theta(c(q)) \,:\, q \in (0,\infty) \}$. From this we deduce that
$\partial^+ \theta(0) \supset  \{ \nabla \theta(c(q)) \,:\, q \in (0,\infty) \}  +  (\mathbb{R}^+_0)^2\,$.
Since $\partial^+ \theta(0)$ is a closed set \cite[Prop.~16.3]{ConvexFunctionalAnalysis-11}, this implies
\[\partial^+ \theta(0) \supset \ol{ \{ \nabla \theta(c(s)) \,:\, q \in (0,\infty) \} } + (\mathbb{R}^+_0)^2\,.\]
Furthermore, for any $w\in \mathbb{R}^2\setminus \{0\}$ with $w_1,\, w_2 \leq 0$ there exists a $p'$ with  
$\theta(0) + \langle (r+w),p\rangle < \theta(p')$
Since $\theta(z)=0$ for $z \in (\{0\} \times \mathbb{R}^+_0) \cup (\mathbb{R}^+_0 \times \{0\})$ and $\theta(z)=-\infty$ outside $[0,\infty)^2$ we finally obtain
that $\theta(0) + \langle r,p\rangle \geq \theta(p)$ if and only if $r \in \ol{ \{ \nabla \theta(c(q)) \,:\, q \in (0,\infty) \} } + (\mathbb{R}^+_0)^2$, which proves the claim.
\end{proof}

\paragraph*{Logarithmic Mean.}
Now, we turn to the specific case when $\theta=\theta_{\tn{log}}$ is the logarithmic mean \eqref{def:means}.
For $s>0$ $\lim_{t \searrow 0} \partial_1 \theta(t,s)= \lim_{t \searrow 0} \partial_2 \theta(s,t) = +\infty$. That is, $N_K(s,0,0)=\{0\} \times \mathbb{R}^-_0 \times \mathbb{R}^-_0$ and analogous  
$N_K(0,s,0)=  \mathbb{R}^-_0 \times \{0\} \times \mathbb{R}^-_0$. Consequently, the algorithm simplifies as follows:
\bigskip

\hspace{10ex} \begin{minipage}{0.8\linewidth}
\begin{algorithmic}
\Function{Project$K$}{$p_1$,$p_2$,$p_3$}
\LineIf{$0 \leq p_3 \leq \theta(p_1,p_2)$} \Return $(p_1,p_2,p_3)$
\LineIf{$p_3 \leq 0$} \Return $(\max\{p_1,0\},\max\{p_2,0\},0)$
\LineIf{$(p_1 \leq 0) \wedge (p_2 \leq 0) \wedge (-p_1/p_3,-p_2/p_3) \in \partial^+ \theta(0)$}
\Return $(0,0,0)$
\State \Return \Call{Project$K$Top}{$p_1$,$p_2$,$p_3$}
\EndFunction
\end{algorithmic}
\end{minipage}
\bigskip

The inclusion in $\partial^+ \theta(0)$  can be tested as follows.
\begin{Lemma}
\label{lem:ProjKLogSupDiff}
Let $z=(z_1,z_2) \in \mathbb{R}^2$. If $\min\{z_1,z_2\} \leq 0$ then $z \notin \partial^+ \theta(0)$. Otherwise, there is a unique $q_1 \in (0,+\infty)$ such that $\partial_1 \theta(q_1^{-1/2},q_1^{1/2})=z_1$ and then
$z \in \partial^+ \theta(0)$ if and only if  $z_2 \geq \partial_2 \theta(q^{-1/2},q^{1/2})$.
\end{Lemma}
\begin{proof}
Note that for the logarithmic mean
$\partial^+ \theta(0) \subset (0,+\infty)^2$
and therefore $z \notin \partial^+ \theta(0)$ if $\min\{z_1,z_2)\leq 0$.
One finds that \[\partial_1 \theta(q^{-1/2},q^{1/2})=\frac{q-1-\log(q)}{\log^2(q)}\] is monotone increasing with $\partial_1 \theta(q^{-1/2},q^{1/2}) \to 0$ as $q \to 0$ and $\partial_1 \theta(q^{-1/2},q^{1/2}) \to +\infty$ as $q \to +\infty$.
Indeed,  for $\beta(q)=\partial_1 \theta(q^{-1/2},q^{1/2})$ with $\beta(1):=\frac12$ we obtain a continuous extension on $(0,\infty)$. Furthermore, we consider 
$\beta'(q) = \frac{2(1-q) + \log(q)(1+q)}{q \log^3(q)}$ with continuous extension $\frac16$ for $q=1$ and verify that $2(1-q) + \log(q)(1+q)$ is negative for $q<1$ and 
positive for  and $q>1$. This implies that $\beta'(q) >0$.
Furthermore, by symmetry we obtain that $\partial_2 \theta(q^{-1/2},q^{1/2})$ is monotone decreasing with $\partial_2 \theta(q^{-1/2},q^{1/2}) \to +\infty$ as $q \to 0$ and $\partial_2 \theta(q^{-1/2},q^{1/2}) \to 0$ as $q \to +\infty$.  
By Lemma \ref{lem:ProjKSupDiff} 
\begin{align*}
\partial^+ \theta(0) = \{ \nabla \theta(q^{-1/2},q^{1/2}) \,:\, q \in (0,+\infty) \} +  (\mathbb{R}_+)^2\,.
\end{align*}	
Thus, for every $z\in (0,+\infty)^2$ there is a unique $q_1\in (0,+\infty)$
such that $\partial_1 \theta(q_1^{-1/2},q_1^{1/2})=z_1$ and $z_1 \geq \partial_1 \theta(q^{-1/2},q^{1/2})$ if and only if $q\leq q_1$.
Furthermore, there is a unique $q_2\in (0,+\infty)$
such that $\partial_2 \theta(q_2^{-1/2},q_2^{1/2})=z_2$ and $z_2 \geq \partial_2 \theta(q^{-1/2},q^{1/2})$ if and only if $q\geq q_2$.
Hence, $z\in \partial^+ \theta(0)$ if and only if  $q_2 \leq q_1$, which is equivalent to $z_2 \geq \partial_2 \theta(q_1^{-1/2},q_1^{1/2})$.
\end{proof}
\begin{Remark}[Comments on Numerical Implementation]
The sought-after $q$ in Lemma \ref{lem:ProjKLogSupDiff} can be determined with a one-dimensional Newton iteration.
The function $q \mapsto \partial_1 \theta(q^{-1/2},q^{1/2})$ becomes increasingly steep as $q \to 0$ which leads to increasingly unstable Newton iterations as $z_1$ approaches $0$. On $q \in [1,+\infty)$ the function is rather flat and easy to invert numerically.
To avoid these numerical problems, note that the roles of $z_1$ and $z_2$ in Lemma \ref{lem:ProjKLogSupDiff} can easily be swapped which corresponds to the transformation $q \leftrightarrow q^{-1}$.
Moreover, for $\max\{z_1,z_2\} < \frac{1}{2}$ one has $z \notin \partial^+ \theta(0)$.
With this rule and by swapping the values of $z_1$ and $z_2$ if $z_1 < z_2$ one can always remain in the regime $q \in [1,+\infty)$.
Additionally, we recommend to replace the function $\theta(s,t)$ and its derivatives by a local Taylor expansion near the numerically unstable diagonal $s=t$.
\end{Remark}

\paragraph*{Geometric Mean.}
Furthermre, let us consider the case where $\theta=\theta_{\tn{geo}}$ is the geometric mean \eqref{def:means}.
For $s>0$ we again find $\lim_{t \searrow 0} \partial_1 \theta(t,s)= \lim_{t \searrow 0} \partial_2 \theta(s,t) = +\infty$ and consequently the same simplification of the algorithm applies as in the case of the logarithmic mean.
For the test of the inclusion $z=(z_1,z_2) \in \partial^+ \theta(0)$, we argue as in the proof of Lemma \ref{lem:ProjKLogSupDiff}.
The functions $\partial_1 \theta(q^{-1/2},q^{1/2}) = \frac{1}{2} q^\frac12$  and  
$\partial_2 \theta(q^{-1/2},q^{1/2}) = \frac{1}{2} q^{-\frac12}$ have the same monotonicity properties as for the logarithmic mean.
Therefore, if $\min\{z_1,z_2\} \leq 0$ then $z \notin \partial^+ \theta(0)$.
Otherwise, $q_1=4\,z_1^2$ and thus the condition $\partial_2 \theta(q_1^{-1/2},q_1^{1/2}) \leq z_2$ is equivalent to $z_1 \cdot z_2 \geq \frac14$.
To summarize, we have obtained
\begin{align*} 
\partial^+ \theta(0) = \left\{ z\in \mathbb{R}^2 \,\colon\, 
z_1 \cdot z_2 \geq \tfrac{1}{4}  \wedge \min\{z_1,z_2\}>0 \right\}\,.
\end{align*}

\subsection[Proximal Mapping of I*J+/-]{Proximal Mapping of ${{\mathcal{I}}^{\ast}_{{\mathcal{J}_{\pm}}}}$}
Note that ${\mathcal{I}}_{\mathcal{J}_{\pm}}$ is a 1-homogeneous function. Hence, ${{\mathcal{I}}^\ast_{{\mathcal{J}_{\pm}}}}$ will once again be an indicator function and $\mathrm{prox}_{{\mathcal{I}}_{\mathcal{J}_{\pm}}^\ast}$ a projection. 
Consequently, the proximal mapping is independent of the step size $\sigma$, i.e.~$\mathrm{prox}_{\sigma\,{\mathcal{I}}_{\mathcal{J}_{\pm}}^\ast}=\mathrm{prox}_{{\mathcal{I}}_{\mathcal{J}_{\pm}}^\ast}$.
To compute $\mathrm{prox}_{{\mathcal{I}}_{\mathcal{J}_{\pm}}^\ast}$ we use Moreau's decomposition \cite[Thm.~14.3]{ConvexFunctionalAnalysis-11} that implies
\begin{align}
\label{eq:ProjJDualMoreau}
\mathrm{prox}_{{\mathcal{I}}_{\mathcal{J}_{\pm}}^\ast}=\id-\mathrm{prox}_{{\mathcal{I}}_{\mathcal{J}_{\pm}}}=\id -\mathrm{proj}_{\mathcal{J}_{\pm}}
\end{align}
where $\id$ is the identity map on $V_{n,h}^0 \times (V_{e,h}^0)^2$. To compute 
$\mathrm{proj}_{\mathcal{J}_{\pm}}(\rho, \massNodes^{-}, \massNodes^{+} )$ for a point 
$(\rho, \massNodes^{-}, \massNodes^{+} ) \in V_{n,h}^0 \times (V_{e,h}^0)^2$ one has to find the minimizer 
$(\rho^{\tn{pr}}, {\massNodes^{-}}^{\tn{pr}}, {\massNodes^{+}}^{\tn{pr}}) \in {\mathcal{J}_{\pm}}$ of
\begin{align*}
\sum_{i=0}^{N-1} \NormV{\rho^{\tn{pr}}(t_i,\cdot)-\rho(t_i,\cdot)}^2 + \NormE{{\massNodes^{-}}^{\tn{pr}}(t_i,\cdot,\cdot)-\massNodes^{-}(t_i,\cdot,\cdot)}^2
+ \NormE{{\massNodes^{+}}^{\tn{pr}}(t_i,\cdot,\cdot)-\massNodes^{+}(t_i,\cdot,\cdot)}^2\,.
\end{align*}
Recall that for any $\rho^{\tn{pr}} \in V_{n,h}^0$ there is precisely one pair $({\massNodes^{-}}^{\tn{pr}},{\massNodes^{+}}^{\tn{pr}}) \in (V_{e,h}^0)^2$ such that $(\rho^{\tn{pr}},{\massNodes^{-}}^{\tn{pr}},{\massNodes^{+}}^{\tn{pr}}) \in {\mathcal{J}_{\pm}}$, see \eqref{def:SetJ}. Therefore, one has to find $\rho^{\tn{pr}} \in V_{n,h}^0$ which minimizes
\begin{align*}
&\sum_{i=0}^{N-1} \sum_{x \in \mathcal{X}} |\rho^{\tn{pr}}(t_i,x) - \rho(t_i,x)|^2 \pi(x)
 + \frac{1}{2}\!\! \sum_{(x,y) \in \mathcal{X}^2} |\rho^{\tn{pr}}(t_i,x) - \massNodes^{-}(t_i,x,y)|^2 Q(x,y) \pi(x) \\
& + \frac{1}{2}\!\! \sum_{(x,y) \in \mathcal{X}^2} |\rho^{\tn{pr}}(t_i,y) - \massNodes^{+}(t_i,x,y)|^2 Q(x,y) \pi(x)\,.
\end{align*}
The optimality condition in $\rho^{\tn{pr}}$ in combination with the reversibility $Q(x,y) \pi(x) = Q(y,x) \pi(y)$ yields for $i=0, \ldots, N-1$, $x \in \mathcal{X}$
\begin{align*}
\rho^{\tn{pr}}(t_i,x) = \frac{1}{2} \left( \rho(t_i,x) + \frac{1}{2} \sum_{y \in \mathcal{X}} (\massNodes^{-}(t_i,x,y) + \massNodes^{+}(t_i,y,x) ) Q(x,y) \right)
\end{align*}
and subsequently ${\massNodes^{-}}^{\tn{pr}}(t_i,x,y)=\rho^{\tn{pr}}(t_i,x)$, ${\massNodes^{+}}^{\tn{pr}}(t_i,x,y)=\rho^{\tn{pr}}(t_i,y)$ for $(x,y) \in \mathcal{X} \times \mathcal{X}$.
Finally, for  $(\rho^{\tn{pr}},{\massNodes^{-}}^{\tn{pr}},{\massNodes^{+}}^{\tn{pr}})=\mathrm{proj}_{{\mathcal{J}_{\pm}}}(\rho,\massNodes^{-},\massNodes^{+})$ using \eqref{eq:ProjJDualMoreau} one gets 
$\mathrm{prox}_{{\mathcal{I}}_{{\mathcal{J}_{\pm}}}^\ast}(\rho,\massNodes^{-},\massNodes^{+})=(\rho,\massNodes^{-},\massNodes^{+})-(\rho^{\tn{pr}},{\massNodes^{-}}^{\tn{pr}},{\massNodes^{+}}^{\tn{pr}})$.


\subsection[Proximal Mapping of I*Javg]{Proximal Mapping of ${{\mathcal{I}}^{\ast}_{{\mathcal{J}_{avg}}}}$}
Once more, we use Moreau's decomposition, \eqref{eq:ProjJDualMoreau}, to compute the proximal mapping of ${\mathcal{I}}^{\ast}_{{\mathcal{J}_{avg}}}$ via the projection onto ${\mathcal{J}_{avg}}$.
Note that the original problem \eqref{def:metricWithUAndV} does not change if we add the constraint $\rho_h(t_0,\cdot)=\rho_A$ and $\rho_h(t_N,\cdot)=\rho_B$ to the set ${\mathcal{J}_{avg}}$. That is, we consider the projection onto the set
\begin{align*}
{\hat{\mathcal{J}}_{avg}} =\left\{ (\rho_h,\bar \rho_h \in {\mathcal{J}_{avg}} \,\colon\,
\rho_h(t_0,\cdot)=\rho_A,\,
\rho_h(t_N,\cdot)=\rho_B \right\}.
\end{align*}
To compute the projection we have to solve
\begin{align*}
  \argmin_{ ( \rho_h^{\tn{pr}}, \bar \rho_h^{\tn{pr}} ) \in {\hat{\mathcal{J}}_{avg}}   } 
  \frac{1}{2} \sum_{i=0}^N  \sum_{x \in \mathcal{X}} |\rho_h^{\tn{pr}}(t_i,x) \!-\! \rho_h(t_i,x)|^2 \pi(x) 
  + \frac{1}{2} \sum_{i=0}^{N-1}  \sum_{x \in \mathcal{X}} |\bar \rho_h^{\tn{pr}}(t_i,x) \!-\! \bar \rho_h(t_i,x)|^2 \pi(x) \, . 
\end{align*}
Thus, we introduce a Lagrange multiplier $\lambda_h \in V_{n,h}^0$ and define the corresponding Lagrangian
\begin{align*}
 L( \rho_h^{\tn{pr}}, \bar \rho_h^{\tn{pr}}, \lambda_h ) 
 & =  \frac{1}{2} \sum_{i=0}^N  \sum_{x \in \mathcal{X}} |\rho_h^{\tn{pr}}(t,x) \!-\! \rho_h(t,x)|^2 \pi(x) 
  + \frac{1}{2} \sum_{i=0}^{N-1}  \sum_{x \in \mathcal{X}} |\bar \rho_h^{\tn{pr}}(t,x) \!-\! \bar \rho_h(t,x)|^2 \pi(x) \\
 & - \sum_{i=0}^{N-1} \sum_{x \in \mathcal{X}} \lambda_h(t_i,x) \left( \avg_h \rho_h^{\tn{pr}}(t_i,x) - \bar \rho_h^{\tn{pr}}(t_i,x) \right) \pi(x)
 \, .
\end{align*}
We know directly from the added boundary constraints that
\begin{align*}
\rho_h^{\tn{pr}}(t_0,x) & = \rho_A, &
\rho_h^{\tn{pr}}(t_N,x) & = \rho_B.
\end{align*}
The optimality condition in $\rho_h^{\tn{pr}}$ for all $x \in \mathcal{X}$ and for interior time steps $i=1,\ldots,N-1$ reads as
\begin{align}\label{eq:AverageMassProjMassNodes}
 \rho_h^{\tn{pr}}(t_i,x) = \rho_h(t_i,x) + \tfrac{1}{2} (\lambda_h(t_{i-1},x) + \lambda_h(t_i,x) ) \, .
\end{align}
Further, the optimality condition in $\bar \rho_h^{\tn{pr}}$ implies that on each interval 
\begin{align}\label{eq:AverageMassProjMassElement}
 \bar \rho_h^{\tn{pr}}(t_i,x) = \bar \rho_h(t_i,x) - \lambda_h(t_i,x) \, .
\end{align}
Combining both with the constraint $\avg_h \rho_h^{\tn{pr}}(t_i,x) = \bar \rho_h^{\tn{pr}}(t_i,x)$, we obtain
\begin{align*}
 \bar \rho_h(t_i,x) - \lambda_h(t_i,x) 
 &= \bar \rho_h^{\tn{pr}}(t_i,x)
 = \avg_h \rho_h^{\tn{pr}}(t_i,x)\\
 &= \avg \rho_h(t_i,x) + \tfrac{1}{4} (\lambda_h(t_{i-1},x) + 2 \lambda_h(t_i,x) + \lambda_h(t_{i+1},x) ) 
\end{align*}
for all interior elements $\interval{i}$ with $i=1,\ldots,N-2$ and for all $x \in \mathcal{X}$.
Analogously, using the boundary conditions we get
\begin{align*}
 & \bar \rho_h(t_0,x) - \lambda_h(t_0,x) = \tfrac{1}{2} (\rho_A(x) + \rho_h(t_1,x) ) + \tfrac{1}{4} ( \lambda_h(t_0,x) + \lambda_h(t_1,x) ) \\
 & \bar \rho_h(t_{N-1},x) - \lambda_h(t_{N-1},x) = \tfrac{1}{2} (\rho_B(x) + \rho_h(t_{N-1},x) ) + \tfrac{1}{4} ( \lambda_h(t_{N-2},x) + \lambda_h(t_{N-1},x) ) \, .
\end{align*}
Thus, for each $x \in \mathcal{X}$ the Lagrange multiplier $\lambda_h$ satisfies the linear system of equations
\begin{align*}
 & \tfrac{1}{4} ( 5 \lambda_h(t_0,x) + \lambda_h(t_1,x) ) = \bar \rho_h(t_0,x) - \tfrac{1}{2} ( \rho_A(x) + \rho_h(t_1,x) ) \\
 & \tfrac{1}{4} ( \lambda_h(t_{i-1},x) + 6 \lambda_h(t_i,x) + \lambda_h(t_{i+1},x) ) = \bar \rho_h(t_i,x) - \tfrac{1}{2} ( \rho_h(t_{i+1},x) + \rho_h(t_i,x) ) \quad \forall i=1,\ldots,N-2 \\
 & \tfrac{1}{4} ( \lambda_h(t_{N-2},x) + 5 \lambda_h(t_{N-1},x) ) = \bar \rho_h(t_{N-1},x) - \tfrac{1}{2} ( \rho_B(x) + \rho_h(t_{N-1},x) )
\end{align*}
This system is solvable, since the corresponding matrix with diagonal $(5,6,\ldots,6,5)$ and off-diagonal $1$ is strictly diagonal dominant. 
Then, given the Lagrange multiplier $\lambda_h$, the solution of the projection problem is given by \eqref{eq:AverageMassProjMassNodes} and \eqref{eq:AverageMassProjMassElement}.
Finally, the proximal map of ${\mathcal{I}}_{{\mathcal{J}_{avg}}}^{\ast}$ can be computed by Moreau's identity, \eqref{eq:ProjJDualMoreau}.
Thus, to compute the proximal mapping of ${{\mathcal{I}}^{\ast}_{{\hat{\mathcal{J}}_{avg}}}}$ one must solve a sparse system in time for each graph node separately. Since the involved matrix is constant, it can be pre-factored.

\subsection[Proximal Mapping of I*J=]{Proximal Mapping of ${\mathcal{I}}_{{\mathcal{J}_{=}}}$}
The proximal map of ${\mathcal{I}}_{{\mathcal{J}_{=}}}$ is  given by the projection
\begin{align*}
 \mathrm{proj}_{{\mathcal{J}_{=}}}( \bar \rho_h, q_h )
& = \argmin_{ ( \bar \rho_h^{{\tn{pr}}}, q_h^{{\tn{pr}}}) \in V_{e,h}^0 \times V_{e,h}^0 \; : \; \bar \rho_h^{\tn{pr}} = q_h^{\tn{pr}} } 
 \quad \frac{1}{2} h \sum_{i=0}^{N-1} \sum_{x,y \in \mathcal{X}}  \left( |\bar \rho_h - \bar \rho_h^{\tn{pr}}|^2 + |q_h - q_h^{\tn{pr}}|^2 \right) Q(x,y) \pi(x)\\
 &= \frac{1}{2} ( \bar \rho_h + q_h, \bar \rho_h + q_h ) \, .
\end{align*}

\section{Numerical Results}\label{sec:numerics}
In what follows we compare the numerical solution based on our
discretization with the explicitly known solution for a simple model
with just two nodes.  Furthermore, we apply our method to a set of
characteristic test cases to study the qualitative and quantitative
behaviour of the discrete transportation distance.
\paragraph*{Comparison with the exact solution for the 2-node case.}
Consider a two point graph $\mathcal{X} = \{a,b\}$ with Markov chain and stationary distribution
\begin{align*}
 Q = \begin{pmatrix}
      0 & p \\
      q & 0
     \end{pmatrix} \, ,
  \quad 
  \pi = \begin{pmatrix}
         \frac{q}{p+q} \\
         \frac{p}{p+q}
        \end{pmatrix} \, ,
\end{align*}
where $p,q \in (0,1]$.
For this case, Maas \cite{Ma11} constructed an explicit solution for the geodesic from 
$\rho_A = \left( \frac{p+q}{q}, 0 \right)$ to $\rho_B = \left( 0, \frac{p+q}{p} \right)$. 
Note that every probability measures on $\mathcal{X}$ can be described by a single parameter $r \in [-1,1]$ via
\begin{align*}
  \rho( r ) = ( \rho_a(r), \rho_b(r) ) = \left( \frac{p+q}{q} \frac{1-r}{2} \, , \frac{p+q}{p} \frac{1+r}{2} \right) .
\end{align*}
Especially, we have $\rho_A = \rho(-1)$ and $\rho_B = \rho(1)$.
Using this representation, Maas showed that for $-1 \leq \alpha \leq \beta \leq 1$
the optimal transport distance is given by
\begin{align}\label{eq:Distance2Point}
 \mathcal{W}( \rho(\alpha), \rho(\beta) ) = \frac{1}{2} \sqrt{\frac{1}{p} + \frac{1}{q}} \int_\alpha^\beta \frac{1}{\sqrt{\theta(\rho_a(r), \rho_b(r) )}} \d r
\end{align}
and the optimal transport geodesic from $\rho(\alpha)$ to $\rho(\beta)$ is given by $\rho(\gamma(t))$ for $t \in [0,1]$, 
where $\gamma$ satisfies the differential equation
\begin{align}\label{eq:ODE2Point}
 \gamma'(t) = 2 (\beta - \alpha) \mathcal{W}( \rho(\alpha), \rho(\beta) ) \sqrt{\frac{pq}{p+q} \theta( \rho(\gamma_a(t)), \rho_b(\gamma(t)) } \, .
\end{align}
For the special case, where $\theta$ is the logarithmic mean
$\theta_{\text{log}}$ and $p=q$, one obtains that
$\theta_{\text{log}} ( \rho_a(r), \rho_b(r) ) =
\frac{r}{\arctanh(r)} \, .$
and consequently the discrete transport distance is given by
$ \mathcal{W}( \rho(\alpha), \rho(\beta) ) = \frac{1}{\sqrt{2p}}
\int_\alpha^\beta \sqrt{ \frac{\arctanh(r)}{r} } \d r\,.  $
Furthermore, the optimal transport geodesic from $\rho(\alpha)$
to $\rho(\beta)$ is given by $\rho(\gamma(t))$ for
$t \in [0,1]$, where $\gamma$ satisfies the differential equation
$ \gamma'(t) = \sqrt{2p} (\beta - \alpha) \mathcal{W}( \rho(\alpha),
\rho(\beta) ) \sqrt{ \frac{\gamma(t)}{\arctanh(\gamma(t))} } \,
.  $
For this two point graph we numerically compute the optimal transport
geodesic. This allows us to evaluate directly the distance $\mathcal{W}$,
which we can compare with a numerical quadrature of
\eqref{eq:Distance2Point}.  Using the approximation of $\mathcal{W}$, we use
an explicit Euler scheme to compute the solution
$\rho_h^{\text{ODE}}$ of the ODE \eqref{eq:ODE2Point}.  For the
case $p=q=1$ we compare our numerical solution to the Euler
approximation for the ODE for $N=2000$ in Fig.~\ref{fig:testone}.
\begin{figure}

\centering
\begin{minipage}{0.8\textwidth}
\resizebox{1.0\textwidth}{!}{
\begin{tikzpicture}[x=\textwidth,y=\textwidth]
 \node at (0.0, 0){\includegraphics[scale=1.0]{fig_Line2_twoPointGeodesic.pdf}};
 \node at (1.0, 0){\includegraphics[scale=1.0]{fig_Line2_twoPointGeodesicDiff.pdf}};
\end{tikzpicture}
}
\end{minipage}
\caption{The mass distribution at $b$ is plotted over time $t\in [0,1]$.
Left: Numerical solution for a 2-point graph $\mathcal{X} = \{a,b\}$ for the logarithmic (red) and geometric (green). The black line represents the diagonal, which 
is the solution in the case of the  (non admissible) arithmetic averaging. 
Right: Difference of the numerical solution for the logarithmic (red) and geometric (green) mean with the Euler scheme solution $\rho_h^{\text{ODE}}$ for the logarithmic mean. 
}
\label{fig:testone}
\end{figure}
\paragraph*{Geodesics on some selected graphs.}
Let us consider four different graphs whose nodes and edges form a triangle, the $3\times 3$ lattice, a cube, and a hypercube, respectively.
Figure~\ref{fig:labels} depicts these graphs with labeled nodes and edges. 
In all cases, we set for each node $x$ with $m$ outgoing edges $\pi(x) = \frac{m}{|E|}$ and $Q(x,y) = \frac{1}{\pi(x)|E|}$.  
Figure~\ref{fig:simple} shows numerically computed geodesic paths. The underlying time step size is $h=\frac1{100}$.
\begin{figure}
\centering
\resizebox{0.8\textwidth}{!}{
\begin{tikzpicture}[x=\textwidth,y=\textwidth]
 \node at (0, 0.45){\includegraphics[page=3, height=0.225\textwidth]{fig_All_LabelingNodesAndEdgesTriangle.pdf}};
 \node at (0, 0.24){\includegraphics[page=3, height=0.225\textwidth]{fig_All_LabelingNodesAndEdgesGrid.pdf}};
 \node at (0, 0){\includegraphics[page=3, height=0.3\textwidth]{fig_All_LabelingNodesAndEdgesCube.pdf}};
  \node at (0.45, 0.2){\includegraphics[page=3, height=0.6\textwidth]{fig_HyperCube_LabelingNodesAndEdges.pdf}};
\end{tikzpicture}
}

\caption{Labeling of nodes and edges for four different graphs: a triangle, the $3x3$ lattice, a cube, and a hypercube.}
\label{fig:labels}
\end{figure}
The solution $(\rho, m)$ is displayed at intermediate
time steps indicated on the arrow in the first row.  For each of these
time steps, blue discs and red arrows superimposed over the graph
display mass and momentum at nodes and on edges, respectively.  The
area of a disc is proportional to the mass $\rho(x) \pi(x)$.  A
red arrow connecting nodes $x$ and $y$ renders the momentum
$m(x,y)$.  The direction of the arrow indicates the direction
of the flow, i.e.~it points from $x$ to $y$ if
$m(x,y) = - m(y,x) > 0$
(cf. Lemma~\ref{Lem:AntisymmetryM}).  The thickness of an arrow is
proportional to $|m(x,y)|Q(x,y)\pi(x)$.  Underneath these
graph drawings both, mass and the momentum on nodes and edges, are
plotted in histograms.  The numbering of the columns in these plots
refers to the numbering of nodes and edges in Figure~\ref{fig:labels}.
The plots associated with $t=0$ and $t=1$ show the prescribed boundary
conditions in time.  As the stopping criteria for the iterative
algorithm in \eqref{eq:ChambollePock} we choose
$\int_0^1\NormV{\rho^{k+1} - \rho^k}^2\d t$ with threshold $10^{-10}$,
where $k$ denotes the iteration step.  Figure~\ref{fig:hypercube}
visualizes in the same fashion an optimal transport path on the graph
of the hypercube.  Note that for the cube, the hypercube, and the
$3\times 3$ lattice the computed solutions are symmetric. In explicit,
mass and momentum at time $t$ equal the mass and the momentum at time
$1-t$ on point reflected nodes and edges, respectively.  Furthermore,
for the cube and the hypercube the distribution of mass is constant on
all nodes at time $t=\frac12$.  Finally, in
Figure~\ref{fig:MChangesSign} we depict an example of graph with four
nodes, which shows that the sign of the momentum variable on a fixed
edge may change along a geodesic path.

\begin{figure}
\resizebox{0.9\textwidth}{!}{
\begin{tikzpicture}[x=\paperwidth,y=\paperwidth]
 \node at (0, 0){\includegraphics[page=1, scale=1]{fig_All_SolutionsGeodesicTriangleGridCube.pdf}};
\end{tikzpicture}
}
\caption{Numerically computed geodesics on a triangle, a square
  lattice and a cube for prescribed boundary conditions at time $0$
  and $1$. Note in particular the symmetry under time reversal and the spreading of mass at intermediate
  times (equidistribution at $t=\frac12$ for the cube).}
\label{fig:simple}
\end{figure}

\begin{figure}
\centering
\resizebox{1.0\textwidth}{!}{
\begin{tikzpicture}[x=\paperwidth,y=\paperwidth]
 \node at (0, 0){\includegraphics[page=1, scale=0.8]{fig_HyperCube_SolutionGeodesicTotal.pdf}};
\end{tikzpicture}
}
\caption{Top: Numerically computed geodesic on a hypercube. Bottom: Distribution of mass and momentum, note again the symmetry under time reversal and the spreading of mass, with equidistribution at time $t=\tfrac12$.}
\label{fig:hypercube}
\end{figure}

\begin{figure}
\centering
\resizebox{1.0\textwidth}{!}{
\begin{tikzpicture}[x=\paperwidth,y=\paperwidth]
 \node at (0, 0){\includegraphics[page=1, scale=0.8]{fig_MChangesSign_SolutionGeodesicTotal.pdf}};
\end{tikzpicture}
}
\caption{Numerically computed geodesic on a graph with four nodes. Note that the sign of $m$ for edge $2$ changes (cf.~$t=\frac{1}{5}$ and $t=\frac{4}{5}$).}
\label{fig:MChangesSign}
\end{figure}

\paragraph*{Experimental results related to the Gromov-Hausdorff convergence for simple graphs.} 
In \cite{GiMa13} it was shown that for the $d$-dimensional torus
$\mathbf{T}^d$ the discrete transportation distance $\mathcal{W}$ on a
discretized torus $\mathbf{T}^d_M$ with uniform mesh size
$\frac{1}{M}$ converges in the Gromov-Hausdorff metric to the
classical $L^2$-Wasserstein distance on $\mathbf{T}^d$.  In fact, the
optimal transport with respect to the classical $L^2$-Wasserstein
distance between two point masses is a point mass travelling along the
connecting straight line.  Concerning the expected concentration of
the transport along this line we perform the following numerical
experiments for $d=1,\,2$.  We first consider for $d=1$ the unit
interval $I = [0,1]$ and a sequence of space discretizations
$\mathcal{X}_M = \{ x_0, \ldots, x_{M} \}$ with uniform mesh size
$\frac1M$ with $M\in \mathbb{N}$.  The corresponding Markov kernel $Q_M$ for
$\mathcal{X}_M$ is defined by
$Q_M(x_i, x_{i+1}) = Q_M(x_i, x_{i-1}) = \frac{1}{2}$ for
$i=1,\ldots, x_{M-1}$ and $Q_M(x_0,x_1) = 1 = Q_M(x_{M},x_{M-1})$.
The continuous $L^2$-Wasserstein geodesic connecting
$\rho_A = \delta_0$ and $\rho_B = \delta_1$ is given by
the transport of the Dirac measure with constant speed:
\begin{align*}
 \rho(t,x) = \delta_t(x) \, .
\end{align*}
In Figure~\ref{fig:GromovHausdorff} we plot the density distribution of the discrete optimal transport geodesic at time $t=\frac{1}{2}$ for different grid sizes 
$\frac1M$. One observes the onset of mass concentration in space at that time at the location $x=\frac12$ for increasing $M$.
\begin{figure}
\centering
\resizebox{0.5\textwidth}{!}{
\begin{tikzpicture}[x=\textwidth,y=\textwidth]
 \node at (0, 0){\includegraphics[scale=3]{fig_convergenceTest_convergenceTest.pdf}};
\end{tikzpicture}
}
\caption{Linearly interpolated densities for the $\mathcal{W}$ geodesic on a one dimensional chain graph between a Dirac mass at the beginning and the end, at $t=0.5$ with $M=2$ (blue), $4$ (red), $8$ (green), $16$ (orange), $32$ (yellow), and $64$ (black).}
\label{fig:GromovHausdorff}
\end{figure}
For $d=2$ we consider a square lattice of uniform grid size $\frac1M$ with $M\in \mathbb{N}$ and nodes $\mathcal{X}_M = \{ (i/M,j/M) \; : \; i,j \in (0,\ldots,M) \}$. 
The weights of the Markov kernel $Q$ are proportional to the number of adjacent edges.
Now, we investigate a discrete geodesic connecting the Dirac masses $\delta_{(0,0)}$ and $\delta_{(1,1)}$. One expects that for increasing $M$ mass on bands parallel to the space diagonal will decrease.  
In Figure \ref{fig:GromovHausdorffSquare} we plot for decreasing mesh size $\frac{1}{M}$ the in time accumulated density values along the diagonal and the off-diagonals bands of nodes.
More precisely, we define the bands of nodes $l^i_M = \{ (x_1,x_2) \in \mathcal{X}_M \times \mathcal{X}_M \; : \; x_2 = x_1 + \frac{i}{M} \}$ ($i=0$ being the diagonal) and compare the values $\int_0^1 \sum_{x \in l^i_M} \rho(t,x) \pi(x) \d t$.
\begin{figure}
\centering
\resizebox{0.9\textwidth}{!}{
\begin{tikzpicture}[x=\textwidth,y=\textwidth]
 \node at (0, 0){\includegraphics[scale=1]{fig_convergenceSquare_SolRhoAllSlices.pdf}};
\end{tikzpicture}
}
\caption{
 Geodesics in the distance $\mathcal{W}$ on a two dimensional grid graph between Dirac masses at diagonally opposite ends. We show accumulated densities along the diagonal and the off-diagonals (see text for details).
From left to right: $M=4$, $8$, $16$, $32$. 
The width of the bars is scaled with the number of lines.
}
\label{fig:GromovHausdorffSquare}
\end{figure}
\paragraph*{Discrete geodesics on an internet network of Europe.}
In Figure~\ref{fig:Internet} we apply the investigated optimal
transport model to a coarse scale internet network of Europe and show
experimental results with masses (data packages) transported from
Dublin, Lisbon, and Madrid to Athens, Stockholm, and Kiev.  Also here,
we set for each node $x$ with $m$ outgoing edges
$\pi(x) = \frac{m}{|E|}$ and $Q(x,y) = \frac{1}{\pi(x)|E|}$, with
$|E|$ the total number of (directed) edges.

\begin{figure}
\resizebox{0.95\textwidth}{!}{
\begin{tikzpicture}[x=\paperwidth,y=\paperwidth]
 \node at (0, 0){\includegraphics[scale=0.6]{fig_Bics.jpg}};
 \node at (0.3, 0){\includegraphics[page=1, scale=1]{fig_InternetGraph_SolutionGeodesic.pdf}};
\end{tikzpicture}
}
\caption{Extraction of a discrete optimal transport geodesic.}
\label{fig:Internet}
\end{figure}


\section{Simulation of the gradient flow of the entropy}\label{sec:gradflow}

The entropy functional on $\mathcal{P}(\mathcal{X})$ is given by
\begin{align*}
 \mathcal{H}(\rho) = \sum_{x \in \mathcal{X}} \rho(x) \log(\rho(x)) \pi(x) \, .
\end{align*}
with the usual convention `$0\log0 = 0$'.  Maas \cite{Ma11} proved
that for the logarithmic mean $\theta_{\tn{log}}(\cdot,\cdot)$ and
$\rho \in \mathcal{P}(\mathcal{X})$ the heat flow
$t \mapsto e^{t \triangle_\mathcal{X}} \rho$ is a gradient flow
trajectory for the entropy $\mathcal{H}(\rho)$ with respect to the
discrete transportation distance $\mathcal{W}$.  In \cite{MaasDiscretePME2014}
it was shown that a similar result holds true for the Renyi entropy
\begin{align*}
 \mathcal{H}_m(\rho) = \frac{1}{m-1} \sum_{x \in \mathcal{X}} \rho(x)^m \pi(x) \, . 
\end{align*}
In fact, for $m=\tfrac12$ and the gradient flow of $\mathcal{H}_m$ with
respect to the metric $\mathcal{W}$ constructed with $\theta$ being the
geometric mean $\theta_{\tn{geom}}(\cdot,\cdot)$ is given by the
Fokker-Planck equation
$\partial_t \rho = \triangle_\mathcal{X} \rho^m$.

To verify this property numerically, we consider a line of five points with stationary distribution $\pi = \frac{1}{5}(1,2,2,2,1)$, Markov kernel $Q(x,y) = \frac{1}{10 \pi(x)}$ for $x,y$ adjacent, 
and initial mass $\rho = \frac{1}{10}(1,1,5,1,1)$.
Following \cite{JKO1998,AmGiSa08}, for an initial density $\rho_0 \in \mathcal{P}(\mathcal{X})$ and a time step size $\tau > 0$ an implicit time-discrete gradient flow scheme for $\mathcal{H}$ can be defined by
\begin{align}\label{eq:GradStep}
 \rho_{k+1} = \argmin_{\rho_B} \frac{1}{2} \mathcal{W}_h(\rho_{k},\rho_B)^2 + \tau \cdot  \mathcal{H}(\rho_B)
\end{align}
with an inner time step size $h$ appearing in the discretization $\mathcal{W}_h$ of $\mathcal{W}$.
To minimize this functional numerically, we simultaneously carry out the external optimization over $\rho$ and the internal optimization within $\mathcal{W}_h$ . 
To this end, we define a discrete continuity equation with one free endpoint. For initial datum $\rho_A \in \mathcal{P}(\mathcal{X})$ let
\begin{align*}
 \mathcal{CE}_h(\rho_A) = \left\{ (\rho_h, m_h, \rho_B) \in  V_{n,h}^{1} \times V_{e,h}^0 \times \mathbb{R}^{\mathcal{X}} \,:\, (\rho_h,m_h) \in \mathcal{CE}_h(\rho_A,\rho_B) \right\}\,.
\end{align*}
Analogous to \eqref{def:metricWithUAndV}, problem \eqref{eq:GradStep} can be written as
\begin{multline*}
 \min
 \; \{ \mathcal{F}(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h,\rho_B) + 
 	\mathcal{G}(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h,\rho_B)
 \,:\,  \\
 (\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h,\rho_B) \in V_{n,h}^{1} \times (V_{e,h}^0)^4 \times (V_{n,h}^0)^2 \times \mathbb{R}^{\mathcal{X}}
 \} 
\end{multline*}
with
\begin{align*}
 \mathcal{F}(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h,\rho_B)  := &
  \widehat{\mathcal{A}}(\vartheta_h, m_h) 
  + {\mathcal{I}}_{{\mathcal{J}_{\pm}}}(q_h,\massNodes^{-}_h, \massNodes^{+}_h)
  + {\mathcal{I}}_{{\mathcal{J}_{avg}}}(\rho_h,\bar \rho_h)
 + 2\,\tau \cdot \mathcal{H}(\rho_B)\,, \\
  \mathcal{G}(\rho_h,m_h,\vartheta_h,\massNodes^{-}_h,\massNodes^{+}_h,\bar \rho_h, q_h,\rho_B)   := &
  {\mathcal{I}}_{\mathcal{CE}_h(\rho_k)}(\rho_h, m_h, \rho_B)
  + {\mathcal{I}}_{\mathcal{K}}(\massNodes^{-}_h,\massNodes^{+}_h,\vartheta_h )
  + {\mathcal{I}}_{{\mathcal{J}_{=}}}(\bar \rho_h,q_h)\,.
\end{align*}
Again, this is amenable for algorithm \eqref{eq:ChambollePock}. 
We extend the space $H$ by a factor $\mathbb{R}^{\mathcal{X}}$ and adapt the scalar product on $H$ \eqref{eq:ScalarProductDiscrete} adding the term $h\,\InProdV{\rho_{B,1}(\cdot)}{\rho_{B,2}(\cdot)}$ with respect or the additional variable $\rho_B$. 
The proximal step of $\mathcal{F}^\ast$ then entails an additional proximal step of $(2\,\tau \cdot \mathcal{H})^\ast$ with respect to $h \NormV{\cdot}$
and in the proximal step of $G$ the projection onto $\mathcal{CE}_h(\rho_A,\rho_B)$ is replaced by a projection onto $\mathcal{CE}_h(\rho_k)$.
Next, we detail these modifications.

\newcommand{\gamma}{\gamma}
Let us recall that the proximal mapping of $(\gamma \cdot \mathcal{H})^\ast$ and $\gamma \cdot \mathcal{H}$ are linked by Moreau's decomposition, cf.~\eqref{eq:ProjJDualMoreau}.
The computation of the the proximal mapping for $\gamma \cdot  \mathcal{H}$ decouples in space and the resulting one dimensional problem can be solved via Newton's method.
This decoupling is possible since we do not enforce the constraint $\rho_B \in \mathcal{P}(\mathcal{X})$ in the formulation of $\mathcal{H}$ but enforce it via the discrete continuity equation constraint.


To implement the projection 
\begin{multline}
\mathrm{proj}_{\mathcal{CE}_h(\rho_A)}(\rho,m,\rho_B) = 
  \argmin_{(\rho^{\tn{pr}},m^{\tn{pr}},\rho^{\tn{pr}}_B) \in \mathcal{CE}_h(\rho_A)}
\frac{h}{2} \sum_{i=0}^N \NormV{\rho^{\tn{pr}}_h(t_i,\cdot)-\rho_h(t_i,\cdot)}^2 \\
+ \frac{h}{2} \sum_{i=0}^{N-1} \NormE{m^{\tn{pr}}_h(t_i,\cdot)-m_h(t_i,\cdot)}^2
+ \frac{h}{2} \NormV{\rho^{\tn{pr}}_B-\rho_B}^2
\label{eq:ProjCEFree}
\end{multline}
onto the set $\mathcal{CE}_h(\rho_A)$ of solutions of the discrete continuity equation with initial data $\rho_A$
the following modifications apply.
Analogous to Proposition \ref{prop:ProjCEDual}, a space time discrete elliptic equation
\begin{align*}\label{eq:ProjCELMFree}
 &\frac{ \varphi_h(t_1,x) - \varphi_h(t_{0},x) }{h^2} + \triangle_\mathcal{X} \varphi_h(t_0,x)
  =  - \left( \frac{\rho_h(t_{1},x) - \rho_A(x)}{h} + \div m_h(t_0,x) \right) \; , \\
 &\frac{ - \tfrac{3}{2} \varphi_h(t_{N-1},x) - \varphi_h(t_{N-2},x)}{h^2} + \triangle_\mathcal{X} \varphi_h(t_{N-1},x)  \\
& \; \qquad\qquad\qquad\qquad \qquad\qquad\qquad = - \left( ( \frac{ \tfrac{1}{2}(\rho_B(x) + \rho_h(t_N,x)) - \rho_h(t_{N-1},x)}{h} + \div m_h(t_{N-1},x) \right) \; ,\\
& \frac{\varphi_h(t_{i+1},x) - 2 \varphi_h(t_i,x) + \varphi_h(t_{i-1},x) }{h^2} + \triangle_\mathcal{X} \varphi_h(t_i,x)\\
&  \; \qquad\qquad\qquad\qquad \qquad\qquad\qquad  = - \left( \frac{\rho_h(t_{i+1},x) - \rho_h(t_i,x)}{h} + \div m_h(t_i,x) \right)  \\
\end{align*}
with $i=1,\ldots,N-2$ and $x \in \mathcal{X}$
has to be solved for the Lagrange multiplier $\varphi_h \in V_{n,h}^0$.
Note that this system is no longer degenerate due to the additional freedom of $\rho_B$ and thus no regularization as before is required.
Then the solution $(\rho^{\tn{pr}},m^{\tn{pr}},\rho^{\tn{pr}}_B)$ to \eqref{eq:ProjCEFree} is given by
\begin{align*}
 \rho_B^{\tn{pr}}(x) &= \frac{1}{2} \left( \rho_h(t_N,x) + \rho_B(x) - \frac{\varphi_h(t_{N-1},x)}{h}  \right) \, , \\
 \rho^{\tn{pr}}_h(t_i,x)  & = \rho_h(t_i,x) + \frac{\varphi_h(t_i,x) - \varphi_h(t_{i-1},x)}{h}   \, , \\
 \rho^{\tn{pr}}_h(t_0,x)  & = \rho_A(x) \, , \; \rho^{\tn{pr}}_h(t_N,x) = \rho_B^{\tn{pr}}(x) \, , \\
 m^{\tn{pr}}_h(t_i,x,y) &= m_h(t_i,x,y) + \nabla_\mathcal{X} \varphi_h(t_i,x,y)    
\end{align*}
for all  $i=1,\ldots,N-2$ and $x,y \in \mathcal{X}$.
\begin{figure}
\resizebox{1.0\textwidth}{!}{
\begin{tikzpicture}[x=\paperwidth,y=0.3\paperwidth]
\node at (0, 0){\includegraphics[page=1, scale=1]{fig_GradientFlow_GFLineTotal.pdf}};
\end{tikzpicture}
}
\caption{Numerical solution of the heat flow (top) and  the Fokker-Planck equation (bottom) based on an explicit Euler scheme (blue) with time step size $10^{-3}$
and for the gradient flow of the associated entropy using the logarithmic mean (red) and the geometric mean (green), respectively, with $\tau=10^{-3}$ and $h=100$.
Panels on the left show the mass distributions on the graph at different times, panels on the right show the values of the entropies over time. }
\label{fig:heatLine}
\end{figure}
In Figure~\ref{fig:heatLine} we compare the numerical results for this natural discretization of the gradient flow of the entropy to 
the flow computed numerically with a simple explicit Euler discretization applied to the heat equation and the Fokker-Planck equation, respectively, with respect to the underlying Markov kernel.


\bibliographystyle{alpha
