# An ordinary differential equation for entropic optimal transport and its
linearly constrained variants

Joshua Zoen-Git Hiew  Luca Nenna and Brendan Pass Department of Mathematical
and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta,
Canada, T6G 2G1. email<EMAIL_ADDRESS>Paris-Saclay, CNRS,
Laboratoire de mathématiques d’Orsay, ParMA, Inria Saclay, 91405, Orsay,
France. email<EMAIL_ADDRESS>of Mathematical
and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta,
Canada, T6G 2G1. email<EMAIL_ADDRESS>
###### Abstract

We characterize the solution to the entropically regularized optimal transport
problem by a well-posed ordinary differential equation (ODE). Our approach
works for discrete marginals and general cost functions, and in addition to
two marginal problems, applies to multi-marginal problems and those with
additional linear constraints. Solving the ODE gives a new numerical method to
solve the optimal transport problem, which has the advantage of yielding the
solution for all intermediate values of the ODE parameter (which is equivalent
to the usual regularization parameter). We illustrate this method with several
numerical simulations. The formulation of the ODE also allows one to compute
derivatives of the optimal cost when the ODE parameter is $0$, corresponding
to the fully regularized limit problem in which only the entropy is minimized.

Keywords. optimal transport, multi-marginal optimal transport, entropic
regularization, ODE, convex analysis.
2020 Mathematics Subject Classification. Primary: 49Q22 ; Secondary: 49N15,
94A17, 49K40.

##  1 Introduction

Given probability measures $\mu^{1}$ and $\mu^{2}$ on domains
$X^{1},X^{2}\subseteq\mathbb{R}^{d}$, respectively, and a cost function
$c:X^{1}\times X^{2}\rightarrow\mathbb{R}$, the Monge-Kantorovich optimal
transport problem is to minimize

$\int_{X^{1}\times X^{2}}c(x^{1},x^{2})d\gamma(x^{1},x^{2})$ (1.1)

over the set $\Pi(\mu^{1},\mu^{2})$ of joint measures on $X^{1}\times X^{2}$
having the $\mu^{i}$ as marginals. This simply stated problem has grown
exponentially in recent years. It has an extremely wide variety of
applications and has spawned many extensions and variants; see the books of
Villani [34, 35] and Santambrogio [33] for detailed surveys.

A particularly popular modification of optimal transport is _entropic
regularization_ , which arises when one penalizes the _entropy_ ,

$H_{\mu^{1}\otimes\mu^{2}}(\gamma)=\int_{X^{1}\times
X^{2}}\frac{d\gamma}{d(\mu^{1}\otimes\mu^{2})}\log\Big{(}\frac{d\gamma}{d(\mu^{1}\otimes\mu^{2})}\Big{)}d(\mu^{1}\otimes\mu^{2}),$

($H$ is taken to be $+\infty$ when $\gamma$ is not absolutely continuous with
respect to $\mu^{1}\otimes\mu^{2}$) in the transport functional, and so
minimizes $\int_{X^{1}\times X^{2}}c(x^{1},x^{2})d\gamma(x^{1},x^{2})+\eta
H_{\mu^{1}\otimes\mu^{2}}(\gamma)$ for some small $\eta>0$. As the
regularization parameter $\eta$ goes to $0$, solutions are well known to tend
to solutions of the unregularized problem. Much of the early motivation for
entropic regularization was computational; in contrast to the linear
unregularized problem, the optimization problem is now strictly convex, and
the unique solution can be efficiently computed via the celebrated Sinkhorn
algorithm [16, 7]. The theoretical properties of the curve $\gamma_{\eta}$ of
solutions, as a function of the regularization parameter $\eta$, including its
rate of convergence and Taylor expansion about $\eta=0$, are also topics of
current research interest [23, 30, 8, 14, 11, 26].

In this paper, we show that, for discrete marginals111General marginals will
be dealt with in a separate, forthcoming paper., solutions to the regularized
problem can be characterized via a well-posed _ordinary differential equation_
(ODE); see Theorem 3.2. In fact, we do this for problems that are considerably
more general than (1.1) in two ways: we allow _multiple marginals_ , as
reviewed in [31, 18] and _additional linear constraints_ as formulated in
[36]. This general framework includes many applications in addition to the
classical, unconstrained two marginal problem, including multi-marginal
problems arising in the evaluation of Wasserstein barycenters [2] and in
density functional theory [9, 15], for example, as well as constrained
problems such as martingale optimal transport, vectorial martingale optimal
transport, adapted optimal transport and relaxed weak optimal transport, among
others [1, 3, 5, 6, 21, 20]. The precise general problem we address, after
regularizing as above, is formulated in (2.3) below.

The development of our ODE results in two distinct contributions. First, it
yields a new numerical method to compute solutions to entropically regularized
optimal transport and its linearly constrained variants, by solving the
corresponding ODE. This is similar in spirit to the method proposed by two of
the present authors in [25], but applies to much more general cost functions.
Note that the work in [25] applies only to multi-marginal problems with
pairwise cost functions, and that the initial condition for the ODE derived
there requires solving $n-1$ two marginal optimal transport problems, whereas
the initial condition in our formulation here has a simple closed form
solution (at least for problems without extra constraints). In particular, the
approach in [25] does not apply to the two marginal problem, and additional
linear constraints were not addressed there. In fact, we would like to
emphasize that the ODE approach seems very promising for computing other
curves of measures. To illustrate this point, we show in Section 5.5 below
that a slight modification of our general framework here can be used to
compute (and provide a new perspective on) Wasserstein geodesics (also known
as displacement interpolants [24]) and barycenters [2].

We verify, through several numerical simulations, the feasability of the
method presented here to compute solutions to a variety of problems. A key
computational advantage of our approach is that it yields the _entire curve_
$\gamma(\varepsilon)$ of solutions, which interpolates between product measure
and optimal transport as the ODE variable $\varepsilon$ varies ($\varepsilon$
is essentially equivalent to $1/\eta$ – see (2.5) below for the precise
definition). This curve is of potential independent interest, due to current
research on the behaviour of $\gamma$ as a function of $\eta$; constructing it
using the Sinkhorn method would require a separate Sinkhorn calculation for
each $\varepsilon$ and would therefore be much less efficient.

A second contribution arising from the formulation of the ODE is that it
allows one to calculate derivatives of the optimal cost around $\varepsilon=0$
(the fully regularized limit, corresponding to $\eta\rightarrow\infty$).
Taylor expansions of the cost around the fully unregularized, optimal
transport limit ($\varepsilon\rightarrow\infty$, or $\eta\rightarrow 0$) have
been an active topic of recent research, whereas expansions around the fully
regularized limit have received relatively limited attention (see [14] for a
first order expansion). The ODE formulation easily allows one to calculate
derivatives (albeit with discrete marginals) in closed form. We illustrate
this by obtaining a formula for the second derivative $C^{\prime\prime}(0)$ of
the optimal cost for unconstrained, two marginal, optimal transport.

Our work here has another, modest and somewhat incidental, consequence. Aside
from developing a general ODE method to solve linearly constrained optimal
transport problems, we also provide a general framework for applying the
Sinkhorn algorithm to this class of problems, unifying and extending previous
work on particular cases [17, 19]. The key to this are two general lemmas on
the strong convexity of the objective functional after an appropriate
reduction in the variables (Lemmas 2.3 and 3.1 below), enabling the
implementation of the general Sinkhorn algorithm (also known as the block
coordinate descent method), found, for instance, in [29] – see Proposition 3.4
below. This fact also justifies the use of the Sinkhorn algorithm in Section 5
to generate reference solutions (especially for the multi-period martingale
problem in Section 5.4, as to the best of our knowledge, a version of the
Sinkhorn algorithm has not been developed for this type of problem before).

This paper is structured as follows. In Section 2, we introduce optimal
transport with extra linear constraints and the entropic regularization of the
corresponding linear programming problem. In this setting, we derive the ODE
and prove its well-posedness in Section 3. In Section 4, we derive formulas
for higher-order derivatives of the optimal value of unconstrained optimal
transport as a function of $\varepsilon$. Finally, in Section 5, we
demonstrate several numerical examples based on our ODE methodology and
compare them with the traditional Sinkhorn algorithm.

##  2 Optimal Transport with Extra Linear constraints

###  2.1 The basic problem

Let $\mu^{1},...,\mu^{n}$ be probability measures222From here on, the
superscript index will denote various marginals or spaces, or different linear
constraints, while the subscript index will represent the elements within
certain marginals or spaces. on respective bounded domains $X^{1},...,X^{n}$
in ${\mathbb{R}}^{d}$. Set $X\coloneqq X^{1}\times...\times X^{n}$. Let
$\mathcal{P}(X)$ be the set of probability measures on $X$ and
$\Pi(\mu^{1},...,\mu^{n})$ be the set of couplings of $\mu^{1},...,\mu^{n}$;
that is, measures $\gamma\in\mathcal{P}(X)$ whose marginals are the $\mu^{i}$.
Let $Q\subset C_{b}(X)$ be a subspace of bounded continuous functions.

If $c:X\to{\mathbb{R}}$ is the cost function, then the optimal transport
problem with extra linear constraints is the following optimization problem:

$\inf_{\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})}\int_{X}cd\gamma,$ (2.1)

where $\Pi^{Q}(\mu^{1},...,\mu^{n})\subset\Pi(\mu^{1},...,\mu^{n})$ is the set
of couplings such that $\int qd\gamma=0$ for all $q\in Q$.

######  Remark 2.1.

The marginal constraints themselves can be represented within the linear
constraint context; letting $Q^{i}=\\{f^{i}-\int f^{i}d\mu^{i}|f^{i}\in
C_{b}(X^{i})\\}$, the subspace $Q\coloneqq\oplus Q^{i}$ of functions enforces
the marginal constraints. We choose here to distinguish these constraints,
denoting them by $\Pi(\mu^{1},...,\mu^{n})$, to emphasize our focus on optimal
transportation. This is contrast to the notational convention in [4], where
they include the marginal constraints within the subspace $Q$.

######  Remark 2.2.

Note that if $Q=\\{0\\}$, then the problem reduces to unconstrained, multi-
marginal optimal transport.

It is easy to see that the existence of an optimizer is equivalent to the non-
emptiness of the feasible set, provided that the function $c$ satisfies
certain mild conditions[36].

The corresponding dual problem is:

$\sup_{(\psi,q)\in\Psi(Q,c)}\sum_{i=1}^{n}\int_{X_{i}}\psi^{i}d\mu^{i},$ (2.2)

where $\Psi(Q,c)\coloneqq\\{(\psi,q)|\sum_{i=1}^{n}\psi^{i}+q\leq
c,\psi=(\psi^{1},...,\psi^{n}),\psi^{i}\in L^{1}(\mu^{i}),q\in Q\\}$

As proven in [36], no duality gap exists between the primal (2.1) and dual
(2.2) problems; that is, the optimal values in the two problems coincide.

For a given $\eta>0$, we are interested here in the following entropically
regularized version of (2.1):

$\boxed{\inf_{\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})}\int_{X}cd\gamma+\eta
H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma),}$ (2.3)

where $H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma)$ denotes the relative entropy of
$\gamma$ with respect to the product measure $\otimes_{i=1}^{n}\mu^{i}$ of the
corresponding marginals, defined as:

$H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma)=\int_{X}\frac{d\gamma}{d(\otimes_{i=1}^{n}\mu^{i})}\log\left(\frac{d\gamma}{d(\otimes_{i=1}^{n}\mu^{i})}\right)d(\otimes_{i=1}^{n}\mu^{i}),$

when $\gamma$ is absolutely continuous with respect to the product measure,
and $+\infty$ otherwise. This regularization has been studied extensively in
the classical optimal transport problem and many variants. The dual of (2.3)
is the unconstrained optimization problem:

$\boxed{\sup_{(\psi,q)\in\Psi(Q,c)}\sum_{i=1}^{n}\int_{X^{i}}\psi^{i}d\mu^{i}-\eta\int_{X}\text{exp}\left(\frac{\sum_{i=1}^{n}\psi^{i}+q-c}{\eta}\right)d(\otimes_{i=1}^{n}\mu^{i}).}$
(2.4)

As the objective function in (2.3) is strictly convex with respect to $\gamma$
for $\eta>0$, the optimization problem has a unique optimal solution, denoted
by $\gamma_{\eta}$. Similarly, as we show below (in the discrete case – see
Lemmas 2.3 and 3.1), if we restrict to an appropriate subspace of $\Psi(Q,c)$
then the objective function in (2.4) is strictly concave with respect to
$(\psi,q)$, also resulting in a unique solution $(\psi_{\eta},q_{\eta})$ when
$\eta>0$. It is well known in the unconstrained setting $Q=\\{0\\}$ that
cluster points of solutions of $\gamma_{\eta}$ as $\eta\rightarrow 0$ solve
the unregularized problem (2.1); furthermore, under strong additional
assumptions, $\gamma_{\eta}$ converges to the unique measure $\gamma_{0}$
having minimal entropy _among all minimizers of ( 2.1)_ (see Theorem 18 in
[22] for a fairly general result in this direction). These results carry over
easily to the more general problem with nontrival $Q$; we provide a proof in
Appendix A for the sake of completeness. We also show there that cluster
points of $(\psi_{\eta},q_{\eta})$ solve (2.2), at least in the discrete
setting that we specialize to below (for convergence in the unconstrained
case, see, for example, [28]). We are interested here in the curve
$\eta\mapsto\gamma_{\eta}$, which interpolates between the optimal transport
$\gamma_{0}$ at $\eta=0$ and the entropy minimizing element of
$\Pi^{Q}(\mu^{1},...,\mu^{n})$ as $\eta\rightarrow\infty$. In fact, it will be
convenient for us to view this interpolation in a slightly different (but
equivalent) way: we fix $\eta$ and introduce a second parameter
$\varepsilon\geq 0$ to obtain the problem

$\inf_{\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})}\int_{X}\varepsilon cd\gamma+\eta
H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma).$ (2.5)

For a fixed $\eta$, as $\varepsilon\to 0$, this becomes equivalent to finding
the couplings that satisfy the extra linear constraints while minimizing the
relative entropy:

$\inf_{\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})}H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma).$

Note that in the unconstrained problem $Q=\\{0\\}$, a trivial solution
arises—the product measure $\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})$ itself.

Conversely, if we let $\varepsilon\to\infty$, the problem reverts to one
equivalent to the original unregularized optimal transport problem with extra
linear constraints (2.1), while at $\varepsilon=1$ we recover (2.3). We are
interested here in the minimmizer $\gamma(\varepsilon)$ in (2.5) as a function
of $\varepsilon$, and will show in Section 3 below that for discrete marginals
it can be characterized by an ODE in $\varepsilon$.

The dual of (2.5) is simply (2.4) with the cost function $\varepsilon c$
replacing $c$:

$\sup_{(\psi,q)\in\Psi(Q,c)}\sum_{i=1}^{n}\int_{X_{i}}\psi^{i}d\mu^{i}-\eta\int_{X}\text{exp}\left(\frac{\sum_{i=1}^{n}\psi^{i}+q-\varepsilon
c}{\eta}\right)d(\otimes_{i=1}^{n}\mu^{i}).$ (2.6)

###  2.2  The discrete regularized problem

We now specialize to the discrete marginal setting. Assume that each
$\mu^{i}=\sum_{x_{i}\in X^{i}}\mu^{i}_{x^{i}}\delta_{x^{i}}$ is supported on a
finite subset $X^{i}\subseteq\mathbb{R}^{d}$; in this case, $C_{b}(X)$ is
finite dimensional. The constraint subspace $Q$ becomes finite dimensional
with basis $\\{q^{j}\\}_{j=1}^{K}$ for some $K\geq 0$ ($K=0$ is unconstrained
optimal transport). The cost function $c=(c_{x})_{x\in X}$, constraint vectors
$q^{j}=(q_{x}^{j})_{x\in X}$ and product measure
$\boldsymbol{\mu}=(\boldsymbol{\mu}_{x})_{x\in
X}=(\otimes_{i=1}^{n}\mu^{i}_{x^{i}})_{x\in X}$ become vectors in
$\mathbb{R}^{m}$ for $m=\Pi_{i=1}^{n}N^{i}=|X|$, where $N^{i}=|X^{i}|$.
Problem (2.5) then becomes

$\begin{array}[]{lrll}\min&\sum_{x\in X}\varepsilon
c_{x}\gamma_{x}&+\eta\sum_{x\in
X}\log\left(\frac{\gamma_{x}}{\boldsymbol{\mu}_{x}}\right)\gamma_{x}&\\\
\text{subject to}&\sum_{x\in X|_{x_{i}}}\gamma_{x}&=\mu_{x_{i}}^{i}&\forall
x_{i}\in X_{i},\quad i=1,...,n\\\ &\sum_{x\in
X}q_{x}^{j}\gamma_{x}&=0&j=1,...,K\\\ &\gamma_{x}&\geq 0&\forall x\in
X,\end{array}$ (2.7)

while its corresponding dual problem becomes an unconstrained optimization
problem:

$\max_{\psi,p}\sum_{i}^{n}\sum_{x_{i}\in
X_{i}}\psi_{x_{i}}^{i}\mu_{x_{i}}^{i}-\eta\sum_{x\in
X}\exp\left(\frac{\sum_{i}^{n}\psi_{x_{i}}^{i}+\sum_{j=1}^{K}q_{x}^{j}p^{j}-\varepsilon
c_{x}}{\eta}\right)\boldsymbol{\mu}_{x}.$ (2.8)

Note that in this setting, (2.7) is known as an exponentially penalized linear
program, and duality is well known (see for instance [13]). It will be
convenient to express the problem more compactly; we enumerate
$X=\\{x_{\ell}\\}_{\ell=1}^{m}$, concatenate $\varphi:=(\psi,p)$ and denote
product measure by $\boldsymbol{\mu}:=\otimes_{i=1}^{n}\mu^{i}$, so that
$\boldsymbol{\mu_{\ell}}=\otimes_{i=1}^{n}\mu^{i}_{x_{l}^{i}}$ for all
$x_{l}=(x^{1}_{l},x^{2}_{l},...,x^{n}_{l})\in X$. We then note that
$\sum_{i=1}^{n}\psi_{x^{i}_{\ell}}^{i}+\sum_{j=1}^{K}q_{x_{\ell}}^{j}p^{j}=B_{\ell}\cdot\psi+C_{\ell}\cdot
p=A_{\ell}\cdot\varphi$ for a vector $A_{\ell}=[B_{\ell},C_{\ell}]$ and define
the matrix $A=[B,C]$ whose rows are the $A_{\ell}$:

$A=\begin{bmatrix}A_{1}\\\ \vdots\\\ A_{m}\end{bmatrix}.$ (2.9)

It is easy to see that the above primal problem (2.7) can be expressed as:

$\begin{array}[]{lrll}\min&\lx@intercol\hfil\varepsilon
c^{T}\gamma+\eta\sum_{\ell=1}^{m}\log(\frac{\gamma_{\ell}}{\boldsymbol{\mu}_{\ell}})\gamma_{\ell}\hfil\lx@intercol&\\\
\text{subject to}&A^{T}\gamma&=b\\\ &\gamma&\geq 0&,\end{array}$ (2.10)

and the corresponding dual (2.8) as:

$\max_{\varphi}b^{T}\varphi-\eta\sum_{\ell=1}^{m}\exp\left(\frac{A_{\ell}\varphi-\varepsilon
c_{\ell}}{\eta}\right)\boldsymbol{\mu}_{\ell},$ (2.11)

where $b=[\mu,0]^{T}=[\mu^{1}\mu^{2}...\mu^{n},0,...,0]^{T}$ is the vector
formed by concatenating the vectors $\mu^{i}\in\mathbb{R}^{N^{i}}$ with a $0$
vector of dimension $K$ (the dimension of $Q$).

The dual program (2.11) is not strictly concave in general. A common strategy
is to restrict the variable $\varphi$ to a subspace where strict concavity
holds. As we will see in the subsequent section, the following lemma implies
that this can always be done without affecting the maximum in (2.11) by fixing
some components of $\varphi$ to be $0$ or, equivalently, removing some columns
from $A$ and considering the reduced problem on the remaining corresponding
components of $\varphi$. In the argument below, it will be useful to note that
the matrix $B$ arises from marginal constraints and $C$ from additional linear
constraints. Therefore, for any $\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})$ we
have $\mu=\gamma^{T}B$ and $0=\gamma^{T}C$.

######  Lemma 2.3.

Assume that $\Pi^{Q}(\mu^{1},\mu^{2},...,\mu^{n}$) is non-empty. Then there is
a matrix $\hat{A}=[\hat{B},\hat{C}]$ and vector $\hat{b}=[\hat{\mu},0]^{T}$ of
the corresponding dimensions such that the function

$\hat{\varphi}\mapsto\hat{b}^{T}\hat{\varphi}-\eta\sum_{\ell=1}^{m}\exp\left(\frac{\hat{A}_{\ell}\hat{\varphi}-\varepsilon
c_{\ell}}{\eta}\right)\boldsymbol{\mu}_{\ell}$ (2.12)

has the same range as the objective function in (2.11), and $\hat{A}$ has full
column rank. Furthermore, $\hat{B}$ and $\hat{C}$ are obtained by removing
some columns of $B$ and $C$ respectively and $\hat{b}$ by removing the
corresponding entries of $b$.

###### Proof.

It is easy to see that $\eqref{discreteOTExtraLinConsEntDualShort}$ and (2.12)
have the same range provided the matrices

$P=\begin{bmatrix}b^{T}\\\ A\\\ \end{bmatrix}=\begin{bmatrix}\mathbf{\mu}&0\\\
B&C\\\ \end{bmatrix}\text{\qquad and
\qquad}\hat{P}=\begin{bmatrix}\hat{b}^{T}\\\ \hat{A}\\\ \end{bmatrix}$

have the same column space.

We first claim that if we drop columns of $B$ until the resulting matrix
$\hat{B}$ has full column rank and the same range as $B$, and drop the
corresponding entries of $\mu$ to obtain $\hat{\mu}$, then
$\begin{bmatrix}\hat{\mu}\\\ \hat{B}\\\ \end{bmatrix}$ has the same range as
$\begin{bmatrix}\mu\\\ B\\\ \end{bmatrix}$. To see this, it suffices to show
that if $\psi$ is in the null space of $B$, it is also orthogonal to $\mu$.
Noting that as $B$ captures the marginal constraints, we have
$B_{\ell}\psi=\sum_{i=1}^{n}\psi_{x_{\ell}^{i}}$, where
$x_{\ell}=(x_{\ell}^{1},x_{\ell}^{2},...,x_{\ell}^{n})$, so if $B\psi=0$,
$\sum_{i=1}^{n}\psi_{x_{\ell}^{i}}=0$ for each $l$.

Multiplying by $\boldsymbol{\mu_{\ell}}$ and summing over $l$, we have

$0=\sum_{l=1}^{m}\boldsymbol{\mu_{\ell}}\sum_{i=1}^{n}\psi_{x_{\ell}^{i}}=\sum_{i=1}^{n}\sum_{l:x_{\ell}^{i}=x^{i}}\boldsymbol{\mu_{\ell}}\psi_{x^{i}}=\sum_{i=1}^{n}\mu^{i}_{x^{i}}\psi_{x^{i}}=\mu\psi,$

establishing the claim.

We now claim that we can drop a subset of the columns of $C$ to obtain a
matrix $\hat{C}$ such that $\begin{bmatrix}\hat{\mu}&0\\\ \hat{B}&C\\\
\end{bmatrix}$ and $\begin{bmatrix}\hat{\mu}&0\\\ \hat{B}&\hat{C}\\\
\end{bmatrix}$ have the same range while $[\hat{B},\hat{C}]$ has full column
rank. Take $q^{1}$ as the first column of $C$. If $q^{1}$ does not belongs to
the column space of $\hat{B}$, we set $\hat{C}^{1}=q^{1}$. Otherwise, we claim
that $\begin{bmatrix}\hat{\mu}\\\ \hat{B}\\\ \end{bmatrix}$ and
$\begin{bmatrix}\hat{\mu}&0\\\ \hat{B}&q^{1}\\\ \end{bmatrix}$ have the same
range. If $\hat{B}\varphi=q^{1}$ and $\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})$,
then $\hat{B}^{T}\gamma=\hat{\mu}^{T}$ and ${q^{1}}^{T}\gamma=0$ by the
marginal constraints and extra linear constraints, respectively. Hence we have

$\hat{\mu}\varphi=\gamma^{T}\hat{B}\varphi=\gamma^{T}q^{1}=0.$

Then for any $\begin{bmatrix}u\\\ v\\\ \end{bmatrix}$ of the appropriate
dimension, we have

$\begin{bmatrix}\hat{\mu}&0\\\ \hat{B}&q^{1}\\\
\end{bmatrix}\begin{bmatrix}u\\\ v\\\
\end{bmatrix}=\begin{bmatrix}\hat{\mu}u\\\ \hat{B}u+q^{1}v\\\
\end{bmatrix}=\begin{bmatrix}\hat{\mu}(u+\varphi v)\\\ \hat{B}(u+\varphi v)\\\
\end{bmatrix}=\begin{bmatrix}\hat{\mu}\\\ \hat{B}\\\
\end{bmatrix}\begin{bmatrix}u+\varphi v\\\ \end{bmatrix},$ (2.13)

where we have used the fact that $\hat{\mu}\varphi=0$ and
$\hat{B}\varphi=q^{1}$ in the second last equation. Therefore, we can conclude
that it is safe to drop this column without changing the range. Continuing in
this way, assume $\hat{C}^{j}$ is a matrix formed by the columns we keep after
considering the first $j$ columns of $C$. If the $j+1$th column $q^{j+1}$ of
$C$ does not belong to the column space of $[\hat{B},\hat{C}^{j}]$, then we
keep it and form a new matrix $\hat{C}^{j+1}$. Otherwise, we can write
$\hat{B}\varphi+\hat{C}^{j}\psi=q^{j+1}$ for some $\varphi$ and $\psi$. Then,
given that $\hat{B}^{T}\gamma=\hat{\mu}^{T}$ and $\hat{C}^{jT}\gamma=0$

$\hat{\mu}\varphi+0\psi=\gamma^{T}\hat{B}\varphi+\gamma^{T}\hat{C}^{j}\psi=\gamma^{T}q^{j+1}=0.$
(2.14)

Then with the similar argument as above, $\begin{bmatrix}\hat{\mu}&0\\\
\hat{B}&\hat{C}^{j}\\\ \end{bmatrix}$ and $\begin{bmatrix}\hat{\mu}&0&0\\\
\hat{B}&\hat{C}^{j}&q^{j+1}\\\ \end{bmatrix}$ have the same range. Then we can
set $\hat{C}^{j+1}=\hat{C}^{j}$ and in any case, $[\hat{B},\hat{C}^{j+1}]$ has
full column rank. We continuous this process until we check all columns of $C$
to obtain the final matrix $\hat{C}=\hat{C}^{K}$. By construction,
$[\hat{B},\hat{C}]$ has full column rank and

$P=\begin{bmatrix}\mu&0\\\ B&C\\\ \end{bmatrix}\text{\qquad and
\qquad}\hat{P}=\begin{bmatrix}\hat{\mu}&0\\\ \hat{B}&\hat{C}\\\ \end{bmatrix}$

have the same range. Now, we set $\hat{b}^{T}=[\hat{\mu},0]$ and
$\hat{A}=[\hat{B},\hat{C}]$, and the result follows. ∎

##  3 ODE for entropic Optimal Transport

###  3.1 Formulation of the ODE

We now formulate a system of ODEs characterizing the optimal dual Kantorovich
potential for problem (2.11). For a fixed $\eta>0$, we define the function:

$\Phi(\varphi,\varepsilon)=-b^{T}\varphi+\eta\sum_{\ell=1}^{m}\exp\left(\frac{A_{\ell}\varphi-\varepsilon
c_{\ell}}{\eta}\right)\boldsymbol{\mu}_{\ell}.$ (3.1)

By Lemma 2.3, we can from now on always assume that $A$ has full column rank.
We then consider the optimization problem:

$\min_{\varphi}\Phi(\varphi,\varepsilon).$ (3.2)

Then, (2.11) is equivalent to (3.2) with the sign inverted. The first order
optimality condition necessitates that at the minimizing $\varphi$:

$\nabla_{\varphi}\Phi(\varphi,\varepsilon)=0.$ (3.3)

It is evident that $\Phi(\varphi,\varepsilon)$ is smooth in $\varphi$ and
$\varepsilon$. The lemma below implies its strict convexity in $\varphi$, and
therefore, for each $\varepsilon$, uniqueness of the solution
$\varphi(\varepsilon)$ to (3.3). Invertibility of the Hessian (also implied by
the following lemma) then ensures $\varphi(\varepsilon)$ is smooth, by the
implicit function theorem.

######  Lemma 3.1.

The function
$\Phi(\varphi,\varepsilon)=-b^{T}\varphi+\eta\sum_{\ell}^{m}\exp(\frac{1}{\eta}(A_{\ell}\varphi-\varepsilon
c_{\ell}))\boldsymbol{\mu}_{\ell}$ is uniformly convex in $\varphi$ over any
compact domain $V\times T\subset{\mathbb{R}}^{e}\times{\mathbb{R}}$ where $e$
is the dimension of $\varphi$.

###### Proof.

A straightforward calculation shows that the gradient and Hessian of
$\Phi(\varphi,\varepsilon)$ in $\varphi$ are given by:

$\displaystyle\nabla_{\varphi}\Phi(\varphi,\varepsilon)$
$\displaystyle=-b^{T}+\sum_{\ell=1}^{m}A_{\ell}^{T}\exp(\frac{1}{\eta}(A_{\ell}\varphi-\varepsilon
c_{\ell}))\boldsymbol{\mu}_{\ell}$ $\displaystyle
D_{\varphi\varphi}^{2}\Phi(\varphi,\varepsilon)$
$\displaystyle=\frac{1}{\eta}\sum_{\ell=1}^{m}A_{\ell}^{T}\exp(\frac{1}{\eta}(A_{\ell}\varphi-\varepsilon
c_{\ell}))A_{\ell}\boldsymbol{\mu}_{\ell}.$

The Hessian can be rewritten as:

$D_{\varphi\varphi}^{2}\Phi(\varphi,\varepsilon)=\frac{1}{\eta}A^{T}DA,$ (3.4)

where $D$ is a diagonal matrix with the $(\ell,\ell)$ th element being
$\exp(\frac{1}{\eta}(A_{\ell}\varphi-\varepsilon
c_{\ell}))\boldsymbol{\mu}_{\ell}$. This implies that $D$ and therefore
$A^{T}DA$ is positive definite, and as the entries are continuous, uniformly
bounded below $A^{T}DA\geq d>0$ on the compact domain $V\times T$. ∎

Differentiating (3.3) with respect to $\varepsilon$ we obtain:

$\frac{{\partial}}{{\partial}\varepsilon}\nabla_{\varphi}\Phi(\varphi(\varepsilon),\varepsilon)+D_{\varphi\varphi}^{2}\Phi(\varphi(\varepsilon),\varepsilon)\frac{d\varphi}{d\varepsilon}(\varepsilon)=0.$

Using the lemma, we obtain the Cauchy problem that governs the behavior of the
optimal potential under the change of $\varepsilon$:

$\begin{cases}\frac{d\varphi}{d\varepsilon}(\varepsilon)&=-\left[D_{\varphi\varphi}^{2}\Phi(\varphi(\varepsilon),\varepsilon)\right]^{-1}\frac{{\partial}}{{\partial}\varepsilon}\nabla_{\varphi}\Phi(\varphi(\varepsilon),\varepsilon)\\\
\varphi(0)&=\phi_{0}.\end{cases}$ (3.5)

Here, $\phi_{0}$ represents the initial condition of the ODE corresponding to
the optimal $\varphi$ solving (3.2) when $\varepsilon=0$. This corresponds to
finding the measure $\gamma$ which minimizes the relative entropy among the
admissable class $\Pi^{Q}(\mu^{1},\mu^{2},...,\mu^{n})$, and is often easier
to find than the solution to (3.2) for $\varepsilon>0$. In particular, note
that when $Q=\\{0\\}$, representing the unconstrained (multi-marginal) optimal
transport problem, the minimizer $\gamma_{0}$ is product measure, making it
straightforward to verify that $\phi_{0}=0$. It is straightforward to
calculate:

$\frac{{\partial}}{{\partial}\varepsilon}\nabla_{\varphi}\Phi(\varphi,\varepsilon)=-\frac{1}{\eta}\sum_{\ell=1}^{m}c_{\ell}A_{\ell}^{T}\exp(\frac{1}{\eta}(A_{\ell}\varphi-\varepsilon
c_{\ell}))\boldsymbol{\mu}_{\ell}.$

Combined with equation (3.4) for the Hessian, the ODE (3.5) can be written
explicitly.

###  3.2 Well-posedness of the ODE

We now establish that the ODE (3.5) has a unique solution. This is a
straightforward application of the Cauchy-Lipschitz Theorem, after the
verification of a few preliminary properties.

As $\varphi(\varepsilon)$ is clearly continuous, there exist an $M>0$ such
that $||\varphi(\varepsilon)||\leq M,\varepsilon\in[0,1]$. Therefore, we aim
to prove the well-posedness of (3.2) over $V\times[0,1]$ where:

$V\coloneqq\left\\{\varphi\in{\mathbb{R}}^{e}|\,||\varphi||\leq M\right\\}.$

######  Theorem 3.2.

Let $\varphi(\varepsilon)$ be the solution of (3.2) for $\varepsilon\in[0,1]$.
Then the trajectory $\varepsilon\mapsto\varphi(\varepsilon)$ is smooth and is
characterized as the unique solution to the initial value problem:

$\begin{cases}\frac{d\varphi}{d\varepsilon}(\varepsilon)&=-\left[D_{\varphi\varphi}^{2}\Phi(\varphi(\varepsilon),\varepsilon)\right]^{-1}\frac{{\partial}}{{\partial}\varepsilon}\nabla_{\varphi}\Phi(\varphi(\varepsilon),\varepsilon)\\\
\varphi(0)&=\phi_{0}\end{cases}.$ (3.6)

Here, $\phi_{0}$ is the solution to the problem:

$\min_{\varphi}\Phi(\varphi,0),$ (3.7)

corresponding to the dual solution to the entropic minimization with respect
to the product measure satisfying the extra linear constraints.

###### Proof.

The smoothness of $\varphi(\varepsilon)$ and the fact that it satisfies the
ODE is demonstrated in section 3.1.

For existence and uniqueness, to apply the Cauchy-Lipschitz Theorem, we only
need that the function

$(\varphi,\varepsilon)\mapsto-\left[D_{\varphi\varphi}^{2}\Phi(\varphi(\varepsilon),\varepsilon)\right]^{-1}\frac{{\partial}}{{\partial}\varepsilon}\nabla_{\varphi}\Phi(\varphi(\varepsilon),\varepsilon)$

is Lipschitz with respect to $\varphi$ and continuous with respect to
$\varepsilon$ on the domain $V\times[0,1]$. This follows immediately from the
smoothness and uniform convexity of $\Phi$ (noting that the latter property
ensures an upper bound on
$[D_{\varphi\varphi}^{2}\Phi(\varphi(\varepsilon),\varepsilon)]^{-1}$). ∎

######  Remark 3.3.

Due to the algebraic properties of the entropy, we can recover the curve of
the solution to the primal $\gamma(\varepsilon)$ by using the optimal dual
variables, that is

$\gamma(\varepsilon)=\exp\left(\frac{\sum_{i}^{n}\psi^{i}(\varepsilon)+\sum_{j=1}^{K}q^{j}p^{j}(\varepsilon)-\varepsilon
c}{\eta}\right)d(\otimes_{i=1}^{n}\mu^{i}).$

The convergence of the Sinkhorn algorithm for unconstrained two or multi-
marginal problems is well known [10]. On the other hand, while the Sinkhorn
algorithm has been applied to various linearly constrained variants [17],
[19], there does not seem to be a general convergence result in this context.
The generalization of the Sinkhorn algorithm to this setting is a block
coordinate descent method, and our Lemma’s 2.3 and 3.1 imply that the entropic
dual problem (2.11) can always be formulated in a way that this scheme
converges, as the following result confirms.

######  Proposition 3.4.

For any initial $\phi_{0}$, the block coordinate descent iteration (Sinkhorn
algorithm) will converge to the unique optimal solution of (3.2) with fixed
$\varepsilon\geq 0$.

###### Proof.

Recall that we are assuming that $A$ has full column rank (and that this
assumption is justified by Lemma 2.3). By Lemma 3.1,
$\Phi(\varphi,\varepsilon)$ is strongly convex in $\varphi$ for fixed
$\varepsilon\geq 0$. Together with the differentiability of $\Phi$ in
$\varphi$, by Theorem 14.6.7 of [29], the result follows. ∎

In particular, note that this result will be crucial in Section 5.4 below, as,
to the best of our knowledge, a version of the Sinkhorn algorithm has not been
implemented on multi-period martingale optimal transport problems before.

##  4 Derivatives of the optimal cost

In this section, we illustrate how the ODE (3.5) can be used to compute
derivatives of the optimal cost at $\varepsilon=0$. For simplicity, we
restrict our attention to the unconstrained two marginal optimal transport
problem (that is, $n=2$, $Q=\\{0\\}$), although the techniques apply more
generally. In this setting, adopting slightly different notation than before,
the primal problem is

$P(\varepsilon)\coloneqq\inf_{\gamma\in\Pi(\mu,\nu)}\int_{X\times
Y}\varepsilon c(x,y)d\gamma(x,y)+\eta H_{\mu\otimes\nu}(\gamma),$ (4.1)

while the dual is:

$D(\varepsilon)\coloneqq\sup_{u,v}\int_{X}ud\mu+\int_{Y}vd\nu-\eta\int_{X\otimes
Y}e^{\frac{u+v-\varepsilon c}{\eta}}d\mu d\nu.$ (4.2)

For the discrete marginals $\\{\mu_{r}\\},\\{\nu_{s}\\}$, if we denote
$c_{rs}=c(x_{r},y_{s})$, then the dual can be written as an equivalent
problem:

$C(\varepsilon)\coloneqq\inf_{u,v}\Phi(u,v,\varepsilon)=-D(\varepsilon)=-P(\varepsilon),$
(4.3)

where

$\Phi(u,v,\varepsilon)=-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\nu_{s}+\eta\sum_{r,s}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}-\varepsilon
c_{rs})\right]\mu_{r}\nu_{s}.$ (4.4)

We will see that one can use the ODE approach to evaluate all derivatives of
$C(\varepsilon)$ at $\varepsilon=0$ in closed form (we provide the explicit
calculations for derivatives up to order $2$, but it will be clear that higher
order derivatives can be calculated similarly). A potential application is a
Taylor expansion for $C(\varepsilon)$ around $\varepsilon=0$ to arbitrary
order.

In this setting, it is actually simpler to eliminate $v$ and express the
problem in terms of $u$ only. Letting $(u(\varepsilon),v(\varepsilon))$ be the
minimizer in (4.4), or equivalently, solutions to our ODE, we have, at
$(u,v)=(u(\varepsilon),v(\varepsilon))$, by the first order optimality
condition, $\forall j$,

$\frac{{\partial}}{{\partial}v_{j}}\Phi=-\nu_{j}+\sum_{r}\exp\left[\frac{1}{\eta}(u_{r}+v_{j}-\varepsilon
c_{rj})\right]\mu_{r}\nu_{j}=0.$

Hence we have

$v_{j}=-\eta\log\left[\sum_{r}\exp\left[\frac{1}{\eta}(u_{r}-\varepsilon
c_{rj})\right]\mu_{r}\right].$ (4.5)

After we substitute (4.5) in (4.4), we get:

$\Phi(u,\varepsilon)\coloneqq-\sum_{r}u_{r}\mu_{r}+\eta\sum_{s}\log\left[\sum_{r}\exp\left[\frac{1}{\eta}(u_{r}-\varepsilon
c_{rs})\right]\mu_{r}\right]\nu_{s}-\eta,$ (4.6)

where we have abused the notation for $\Phi$; note that we can now write
$C(\varepsilon)=\inf_{u}\Phi(u,\varepsilon)$.

For simplicity, we denote $e_{ij}=\exp[\frac{1}{\eta}(u_{i}-\varepsilon
c_{ij})]\mu_{i}$. Note that,

$\displaystyle-P^{\prime}(\varepsilon)=-D^{\prime}(\varepsilon)=C^{\prime}(\varepsilon)$
$\displaystyle=\frac{d}{d\varepsilon}\Phi(u(\varepsilon),\varepsilon)$
$\displaystyle=\nabla_{u}\Phi(u(\varepsilon),\varepsilon)\frac{du}{d\varepsilon}+\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)$
$\displaystyle=\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)$
$\displaystyle=-\sum_{s}\left[\frac{\sum_{r}e_{rs}c_{rs}}{\sum_{r}e_{rs}}\right]\nu_{s},$
(4.7)

where we have used the first order optimality condition
$\nabla_{u}\Phi(u(\varepsilon),\varepsilon)=0,\forall\varepsilon\geq 0$. Since
it is easy to verify that $\forall i,u_{i}(0)=0$, we have $e_{ij}=\mu_{i}$
when $\varepsilon=0$, therefore:

$C^{\prime}(0)=-\sum_{s}\left[\frac{\sum_{r}c_{rs}\mu_{r}}{\sum_{r}\mu_{r}}\right]\nu_{s}=-{\mathbb{E}}[c(X,Y)].$
(4.8)

The expectation in (4.8) is taken with respect to the product measure
$\mu\otimes\nu$. 444Note that (4.8) could more easily be seen by applying the
envelope theorem to the primal problem (4.1) and noting the unique optimality
of product measure at $\varepsilon=0$ in that equation. However, we will
require (4.7) to compute higher order derivatives.

For the second order derivative,

$\displaystyle C^{\prime\prime}(\varepsilon)$
$\displaystyle=\frac{d}{d\varepsilon}\left(\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right)$
$\displaystyle=\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]^{T}\frac{du}{d\varepsilon}+\frac{{\partial}^{2}}{{\partial}\varepsilon^{2}}\Phi(u(\varepsilon),\varepsilon)$
$\displaystyle=\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]^{T}\left[-D^{2}_{uu}\Phi(u(\varepsilon),\varepsilon)\right]^{-1}\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]+\frac{{\partial}^{2}}{{\partial}\varepsilon^{2}}\Phi(u(\varepsilon),\varepsilon).$
(4.9)

The last equality is due to the ODE (3.5). We note that evaluating this
quantity at the point $\varepsilon=0$ (where $u(0)=0$ is known) involves only
derivatives of $\Phi$. We relegate the detailed calculation to Appendix B and
provide only the result here:

$C^{\prime\prime}(0)=\frac{1}{\eta}\biggl{(}\left({\mathbb{E}}[c(X,Y)]\right)^{2}+{\mathbb{E}}[c^{2}(X,Y)]-{\mathbb{E}}\left[({\mathbb{E}}[c(X,Y)|X])^{2}\right]-{\mathbb{E}}\left[({\mathbb{E}}[c(X,Y)|Y])^{2}\right]\biggr{)},$
(4.10)

where the expectations are with respect to product measure $\mu\otimes\nu$,
or, equivalently,

$\displaystyle C^{\prime\prime}(0)=$
$\displaystyle\frac{1}{\eta}\Biggl{(}\left(\int_{X\times
Y}c(x,y)d\mu(x)d\nu(y)\right)^{2}+\int_{X\times Y}c^{2}(x,y)d\mu(x)d\nu(y)$
$\displaystyle-$
$\displaystyle\int_{Y}\left(\int_{X}c(x,y)d\mu(x)\right)^{2}d\nu(y)-\int_{X}\left(\int_{Y}c(x,y)d\nu(y)\right)^{2}d\mu(x)\Biggr{)}.$

We note that higher order derivatives can be evaluated similarly; for
instance, one can find $C^{\prime\prime\prime}(\varepsilon)$ by
differentiating (4.9) using the chain rule, and then use the ODE (3.5) to
eliminate the $\frac{du}{d\varepsilon}$ terms and substitute $u(0)=0$ to
determine $C^{\prime\prime\prime}(0)$. The result is a closed form expression
involving only derivatives of the function $\Phi$ at $(u,\varepsilon)=(0,0)$,
which can readily be computed.

######  Remark 4.1.

$P(\varepsilon)$ as defined by (4.1) is clearly concave, as an infimum of
affine functions. Therefore, we must have
$C^{\prime\prime}(0)=-P^{\prime\prime}(0)\geq 0$.

We do not see another way to prove that the expression in (4.10) is non-
negative. However, the fact that expectations are taken with respect to
product measure is crucial, as the expression can in fact be negative for more
general couplings. To see this, note that

$\displaystyle\text{Var}[{\mathbb{E}}[c(X,Y)|X]]$
$\displaystyle={\mathbb{E}}[({\mathbb{E}}[c(X,Y)|X])^{2}]-(E[c])^{2}$
$\displaystyle{\mathbb{E}}[\text{Var}[c(X,Y)|X]]$
$\displaystyle={\mathbb{E}}[c^{2}]-{\mathbb{E}}[({\mathbb{E}}[c(X,Y)|X])^{2}],$

and consider the case $\mu=\nu$ such that $\mu$ is supported on more than one
point, $c(x,y)=x+y$ and the coupling is given by $X=Y$. In this case,
${\mathbb{E}}[({\mathbb{E}}[c(X,Y)|Y])^{2}]={\mathbb{E}}[({\mathbb{E}}[c(X,Y)|X])^{2}]$
and the non-negativity of the quantity in (4.10) becomes

${\mathbb{E}}[\text{Var}[c(X,Y)|X]]\geq\text{Var}[{\mathbb{E}}[c(X,Y)|X]].$

However, it is easy to calculate that ${\mathbb{E}}[\text{Var}[c(X,Y)|X]]=0$
but $\text{Var}[{\mathbb{E}}[c(X,Y)|X]]=4\text{Var}[X]\neq 0$, which violates
the above inequality.

##  5 Numerical examples

In this section we exploit the ODE (3.5) to compute numerical approximations
of solutions to the optimal transport problem and several variants. In all the
examples below, we discretize the ODE in $\varepsilon$ and solve using a 4th-
order Runge-Kutta method.

To verify the accuracy of our new numerical method, we also compute solutions
via the Sinkhorn algorithm in each example in the first four subsections
below; in all cases, the values of the optimal costs obtained from the
Sinkhorn are very close to the values obtained from the ODE method555 A
similar comparison is not provided for the example in the last subsection,
since its purpose is somewhat different; in particular, one of the marginals
of the curve $\gamma(\varepsilon)$ is in fact the main object of interest
there, rather than the final value $\gamma(1)$, and this curve of marginals is
not as readily obtained by the Sinkhorn method.. Furthermore, we have chosen
examples for which explicit solutions $\gamma_{0}$ of the unregularized
problem (2.1) are known. This then yields a lower bound on the optimal value
in (2.3). Since $\gamma_{0}$ is a competitor in (2.3), it can also be used to
compute an upper bound for this problem; thus, we obtain an interval in which
the optimal value must fall:

$\int_{X}cd\gamma_{0}\leq\inf_{\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})}\int_{X}cd\gamma+\eta
H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma)\leq\int_{X}cd\gamma_{0}+\eta
H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma_{0}).$ (5.1)

In the examples in the first four subsections below, the values obtained via
both the Sinkhorn and ODE method do indeed fall with this interval.

Throughout this section, we will consider problems with at most three
marginals which will be represented by the symbols $\mu$, $\nu$ and $\theta$.
We will denote the corresponding spaces $X$, $Y$ and $Z$. We perform all
numerical calculations in Python on 13th Gen Intel(R) Core(TM) i7-13620H 2.40
GHz Notebook.

###  5.1 Two marginal optimal transport

We begin with the two marginal optimal transport problem. Let
$\mu=\sum_{r}\mu_{r}\delta_{x_{r}},\nu=\sum_{s}\nu_{s}\delta_{y_{s}}$ be
discrete measures and let the cost matrix be $c_{rs}=c(x_{r},y_{s})$. For
fixed $\eta>0$, we aim to solve:

$\inf_{u,v}\Phi(u,v,\varepsilon),$ (5.2)

where

$\Phi(u,v,\varepsilon)=-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\nu_{s}+\eta\sum_{r,s}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}-\varepsilon
c_{rs})\right]\mu_{r}\nu_{s}.$ (5.3)

It well known that if we fix $u_{0}=v_{0}=0$, then there will be a unique
minimizer of (5.3). This corresponds to removing the corresponding column to
the $A$ in (3.1) so that the reduced version has full column rank, as
discussed more generally in Lemma 2.3.

As in the previous section, it will be convenient to further reduce the number
of variables, using the first order optimality condition:

$\displaystyle v_{j}$
$\displaystyle=-\eta\log\left[\sum_{r}\exp\left[\frac{1}{\eta}(u_{r}-\varepsilon
c_{rj})\right]\mu_{r}\right]=:v_{j}(u).$ (5.4)

We simplify by substituting (5.4) into (5.3) to yield the new dual objective
function:

$\bar{\Phi}(u,\varepsilon)=\Phi(u,v(u),\varepsilon)=-\sum_{r}u_{r}\mu_{r}+\eta\sum_{s}\log\left[\sum_{r}\exp\left[\frac{1}{\eta}(u_{r}-\varepsilon
c_{rs})\right]\mu_{r}\right]\nu_{s}+\eta.$ (5.5)

The unique minimizer $u(\varepsilon)$ of $u\mapsto\bar{\Phi}(u,\varepsilon)$
will clearly still satisfify the ODE (3.5) (with $\bar{\Phi}$ in place of
$\Phi$).

######  Remark 5.1.

The ODE induced by the function (5.5) inherits the well-posedness related to
the original function (5.3). The gradient and the Hessian of the function
(5.5) are:

$\displaystyle\nabla\bar{\Phi}(u,\varepsilon)$
$\displaystyle=\nabla_{u}\Phi(u,v(u),\varepsilon)+\nabla_{v}\Phi(u,v(u),\varepsilon)Dv(u)$
$\displaystyle D^{2}\bar{\Phi}(u,\varepsilon)$
$\displaystyle=D_{uu}^{2}\Phi(u,v(u),\varepsilon)+D_{uv}^{2}\Phi(u,v(u),\varepsilon)Dv(u)+Dv(u)^{T}D_{vu}^{2}\Phi(u,v(u),\varepsilon)$
$\displaystyle\qquad+Dv(u)^{T}D_{vv}^{2}\Phi(u,v(u),\varepsilon)Dv(u)+\nabla_{v}\Phi(u,v(u),\varepsilon)D^{2}v(u).$

The last term vanishes since $\nabla_{v}\Phi=0$. Therefore, we have:

$D^{2}\bar{\Phi}=[I,Dv^{T}]D^{2}\Phi\begin{bmatrix}I\\\ Dv\\\ \end{bmatrix}.$
(5.6)

$D^{2}\bar{\Phi}$ then inherits lower bounds from $D^{2}\Phi$, and so
$\bar{\Phi}$ is uniformly convex on compact sets. The well posedness of the
ODE satisfied by $u$ then follows exactly as in Theorem 3.2. We will abuse
notation below by referring to the function of $u$ and $\varepsilon$ only as
$\Phi$ (rather than $\bar{\Phi}$ as above), and, in subsequent subsections, we
will make similar reductions in the number of variables; the well-posedness of
the corresponding ODEs will follow similarly and we will not comment on it in
detail.

Let

$e_{ij}=\exp\left[\frac{1}{\eta}(u_{i}-\varepsilon
c_{ij})\right]\mu_{i},\qquad
e_{ij}^{c}=\exp\left[\frac{1}{\eta}(u_{i}-\varepsilon
c_{ij})\right]c_{ij}\mu_{i}.$

Then the Hessian matrix for $\Phi$ will be the symmetric matrix:

$\biggl{(}D_{uu}^{2}\Phi(u(\varepsilon),\varepsilon)\biggr{)}_{i_{1}i_{2}}=\frac{1}{\eta}\sum_{s}\left[\frac{e_{i_{1}s}}{\sum_{r}e_{rs}}\delta_{i_{1}i_{2}}-\frac{e_{i_{1}s}\cdot
e_{i_{2}s}}{(\sum_{r}e_{rs})^{2}}\right]\nu_{s},$ (5.7)

where $\delta_{i_{1}i_{2}}$ is the Kronecker delta and the mixed derivative
is:

$\biggl{(}\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\biggr{)}_{i}=-\frac{1}{\eta}\sum_{s}\left[\frac{e_{is}^{c}}{\sum_{r}e_{rs}}-\frac{e_{is}\cdot
e_{s}^{c}}{(\sum_{r}e_{rs})^{2}}\right]\nu_{s}.$ (5.8)

It is easy to verify that the initial condition for the system of equations is
$u_{i}=0,\forall i$.

We take $\eta=0.002$ and discretize $[0,1]$ with 100 time steps for the
variable $\varepsilon$. Both the marginals $\mu$ and $\nu$ are uniformly
supported on 100 evenly spaced grid points on $[0,1]$. We solved the ODE with
two costs: the attractive cost $c(x,y)=(y-x)^{2}$ and the repulsive cost
$c(x,y)=-\log(0.1+|x-y|)$. It is well known that the optimal measure for the
attractive cost is induced by the identity map in our setting. A closed form
for the optimal measure for the repulsive cost is also known [12]. We provide
the evolution of the primal coupling at 16 evenly distributed moments for
$\varepsilon$ between 0 and 1 for both attractive and repulsive type costs in
Figure 1, and 2, respectively. Note that the solutions at $\varepsilon=1$
agree qualitatively with the known solutions for both costs. Tables 1 and 2
summarize the corresponding numerical calculations for the attractive and
repulsive type costs, respectively. In both cases, the optimal values are
similar to those found by the Sinkhorn algorithm, and lie within the range
given by (5.1).

Figure 1: Evolution of the primal coupling at 16 evenly distributed values of $\varepsilon$ between 0 and 1 calculated by the ODE method for two marginal optimal transport with attractive cost $c(x,y)=|y-x|^{2}$. Figure 2: Evolution of the primal coupling at 16 evenly distributed values of $\varepsilon$ between 0 and 1 calculated by the ODE method for two marginal optimal transport with repulsive cost $c(x,y)=-\log(0.1+|x-y|)$. | ODE | Sinkhorn
---|---|---
Computed regularized primal value | 0.0050 | 0.0052
Optimal unregularized primal value | 0 | 0
Optimal unregularized primal value + Entropy | 0.0092 | 0.0092
Iterations | 100 | 324
CPU time (sec) | 1.40 | 0.04
Table 1: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for two marginal optimal transport with attractive cost $c(x,y)=(y-x)^{2}$. | ODE | Sinkhorn
---|---|---
Computed regularized primal value | 0.5033 | 0.5080
Optimal unregularized primal value | 0.5024 | 0.5024
Optimal unregularized primal value + Entropy | 0.5117 | 0.5117
Iterations | 100 | 684
CPU time (sec) | 1.49 | 0.15
Table 2: Comparison of the performance between 4-th order Runge-Kutta ODE
method and Sinkhorn algorithm for two marginal optimal transport with
repulsive cost $c(x,y)=-\log(0.1+|x-y|)$.

###  5.2 Multi-marginal optimal transport

We turn now to a three marginal optimal transport problem; in this setting our
approach here is closely related to the work in [25]. However, the present
framework can accommodate completely general cost functions (whereas the work
in [25] applies only to pairwise costs). Even for pairwise costs, the solution
$\gamma(\varepsilon)$ as a function of the regularization parameter here is
likely of more direct interest than the interpolation between two and multi-
marginal solutions produced by the method in [25].

Let
$\mu=\sum_{r}\mu_{r}\delta_{x_{r}},\theta=\sum_{s}\theta_{s}\delta_{y_{s}},\nu=\sum_{t}\nu_{t}\delta_{z_{t}}$
and $c_{rst}=c(x_{r},y_{s},z_{t})$ be the cost tensor. For fixed $\eta>0$, we
minimize the following objective function:

$\Phi(u,v,w,\varepsilon)=-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\theta_{s}-\sum_{t}w_{t}\nu_{t}+\eta\sum_{r,s,t}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}+w_{t}-\varepsilon
c_{rst})\right]\mu_{r}\theta_{s}\nu_{t}.$ (5.9)

As in the previous section, we substitute the first order condition with
respect to one variable,

$\displaystyle w_{k}$
$\displaystyle=-\eta\log\left[\sum_{rs}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}-\varepsilon
c_{rsk})\right]\mu_{r}\theta_{s}\right],$ (5.10)

into (5.9) to reduce the number of unknown variables:

$\Phi(u,v,\varepsilon)=-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\theta_{s}+\eta\sum_{t}\log\left[\sum_{r,s}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}-\varepsilon
c_{rst})\right]\mu_{r}\theta_{s}\right]\nu_{t}+\eta.$ (5.11)

As in the two marginal case (see Remark 5.1) , the minimizers
$u(\varepsilon,v(\varepsilon))$ will satisfy the well-posed ODE (3.5).

We refer the reader to appendix B.1 for the Hessian of $\Phi$ with respect to
the dual potential variables $\varphi=(u,v)$ and the mixed derivative with
respect to $(u,v)$ and $\varepsilon$.

We set the initial values to be $u_{i}=v_{j}=0,\forall i,j$ which can be
verified to be the unique solution for $\varepsilon=0$ satisfying our
normalization $u_{0}=v_{0}=0$. For the numerical simulation, each marginal to
be uniform on 99 evenly spaced grid points (for simplicity of calculating the
true optimal measure [12]) on $[0,1]$, $\eta=0.006$ and $\varepsilon\in[0,1]$
with 100 time steps. The cost function is $c(x,y,z)=d(x,y)+d(y,z)+d(x,z)$
where $d(x,y)=-\log(0.1+|x-y|)$, which is the same as the first example in
[25]. Figure 3 shows the 2D projection to all the pairs of marginals at
$\varepsilon=1$, which agrees well with the known explicit solution to the
unregularized problem in [12]. From Table 4, the values obtained from the ODE
and Sinkhorn method agree quite well and both fall within the interval given
by (5.1). Note in addition that in this case, the ODE method actually takes
less running time and fewer iterations than the Sinkhorn.

Figure 3: The optimal measure for three marginal optimal transport with cost $c(x,y,z)=d(x,y)+d(y,z)+d(x,z)$ where $d(x,y)=-\log(0.1+|x-y|)$ generated by ODE | ODE | Sinkhorn
---|---|---
Computed regularized primal value | 1.9163 | 1.9193
Optimal unregularized primal value | 1.9137 | 1.9137
Optimal unregularized primal value + Entropy | 1.9647 | 1.9647
Iterations | 100 | 1006
CPU time (sec) | 12.40 | 47.46
Table 3: Comparison of the performance between 4-th order Runge-Kutta ODE
method and Sinkhorn algorithm for three marginal optimal transport problem.

###  5.3 Martingale optimal transport

Given probability measures $\mu$ and $\nu$, with $\mu$ dominated by $\nu$ in
convex order (denoted by $\mu\preceq_{c}\nu)$, the (one-period) martingale
optimal transport falls into the class of optimal transport problems with
extra linear constraints (2.1) when we take the subspace $Q$ to be the set $M$
of martingale test functions:

$M=\\{g(x)(y-x)|g(x)\in C_{b}(X)\\}.$ (5.12)

Note that $\int md\pi=0,\forall m\in M$ is equivalent to the classical
martingale condition:

$\int_{Y}yd\pi_{x}(y)=x,$ (5.13)

where $\pi=\pi_{x}\otimes\mu$ is the regular disintegration of $\pi$ with
respect to $\mu$, or in probabilistic notation:

${\mathbb{E}}_{\pi}[Y|X]=X.$ (5.14)

The primal problem of martingale optimal transport problem is:

$\inf_{\pi\in\Pi^{M}(\mu,\nu)}\int_{X\times Y}c(x,y)d\pi(x,y),$ (5.15)

and the corresponding dual problem:

$\sup\int_{X}u(x)d\mu(x)+\int_{Y}v(y)d\nu(y),$ (5.16)

subject to $u(x)+v(y)+g(x)(y-x)\leq c(x,y)$, among all $u\in L(\mu),v\in
L(\nu),g\in C(X)$. For the numerical example, let
$\mu=\sum_{r}\mu_{r}\delta_{x_{r}},\nu=\sum_{s}\nu_{s}\delta_{y_{s}}$ be
discrete measures such that $\mu$ is less than $\nu$ in convex order. Then the
cost matrix will be $c_{rs}=c(x_{r},y_{s})$. The entropic regularized problem
with fixed $\eta$ will be:

$\inf_{u,v,g}\Phi(u,v,g,\varepsilon),$ (5.17)

where

$\Phi(u,v,\varepsilon)=-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\nu_{s}+\eta\sum_{r,s}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}+g_{r}(y_{s}-x_{r})-\varepsilon
c_{rs})\right]\mu_{r}\nu_{s}.$ (5.18)

The first order optimality condition implies

$\displaystyle u_{i}$
$\displaystyle=-\eta\log\left[\sum_{s}\exp\left[\frac{1}{\eta}(v_{s}+g_{i}(y_{s}-x_{i})-\varepsilon
c_{is})\right]\nu_{s}\right]$ (5.19) $\displaystyle v_{j}$
$\displaystyle=-\eta\log\left[\sum_{r}\exp\left[\frac{1}{\eta}(u_{r}+g_{r}(y_{j}-x_{r})-\varepsilon
c_{rj})\right]\mu_{r}\right]$ (5.20) $\displaystyle 0$
$\displaystyle=\sum_{s}\exp\left[\frac{1}{\eta}(u_{i}+v_{s}+g_{i}(y_{s}-x_{i})-\varepsilon
c_{is})\right](y_{s}-x_{i})\nu_{s}.$ (5.21)

We will use this set of equations in the Sinkhorn algorithm. One important
note should be made towards the fitting of $g_{i}$ from $u_{i}$ and $v_{i}$ is
not as straight forward as in the unconstrained optimal transport. Instead, we
need to solve the root of the equation (5.21) for each $i$. Slightly longer
computation time is expected if we apply Newton’s method for finding the root.

On the other hand, we can further simplify (5.18) by substituting (5.19) into
it. By abuse of notation, the dual cost becomes

$\Phi(v,g,\varepsilon)\coloneqq-\sum_{s}v_{s}\nu_{s}+\eta\sum_{r}\log\left[\sum_{s}\exp\left[\frac{1}{\eta}(v_{s}+g_{r}(y_{s}-x_{r})-\varepsilon
c_{rs})\right]\nu_{s}\right]\mu_{r}+\eta.$ (5.22)

The minimizers $v(\varepsilon),g(\varepsilon)$ satisfy the well-posed ODE
(3.5).

We refer the reader to Appendix B.2 for the Hessian of $\Phi$ with respect to
the dual potential variables $v,g$ and the mixed derivatives with respect to
$v,g$ and $\varepsilon$. Unlike in the unconstrained optimal transport
problem, $v_{j}=g_{i}=0$ is no longer an initial condition for the system of
equations since the product measure is in general not a martingale measure
(nevertheless, setting $v_{0}=g_{0}=0$ will reduce $A$ from (2.9) into a full
column rank matrix). Instead, we can solve the initial system of equations,
(5.19) - (5.21) with $\varepsilon=0$, by the Sinkhorn algorithm. Note that by
our numerical experiment, the running time for the Sinkhorn algorithm is much
faster for $\varepsilon=0$ then $\varepsilon=1$. This keeps the overall
calculation time for the ODE method competitive with the pure Sinkhorn method
to minimize the $\varepsilon=1$ objective function.

The cost function is $c(x,y)=\exp(-x)\cdot y^{2}$ which is of martingale
Spence–Mirrlees type, hence the unregularized optimal coupling is the left-
monotone coupling [6]. The $x$-marginal $\mu$ is uniformly distributed on 100
evenly spaced grid points on $[-0.3,0.3]$ while the $y$-marginal $\nu$ is
uniformly distributed on 200 evenly spaced grid points on $[-1,1]$. We set
$\eta=0.006$ and $\varepsilon=1$ for this numerical example. We employ 25
steps to discretize $\varepsilon$ on $[0,1]$. The plot in Figure 4, the plot
shows that the ODE method does indeed approximate the left-monotone measure.
From Table 4, we see that both the ODE method and Sinkhorn algorithm yield
comparable optimal values, and both fall inside the desired interval (5.1).
Again, in this case, the ODE method takes fewer iterations and less running
time.

Figure 4: The optimal measure for martingale optimal transport with cost $c(x,y)=\exp(-x)\cdot y^{2}$ generated by ODE | ODE | Sinkhorn
---|---|---
Computed regularized primal value | 0.2990 | 0.2990
Optimal unregularized primal value | 0.2964 | 0.2964
Optimal unregularized primal value + Entropy | 0.3211 | 0.3211
Iterations | 25 | 152
CPU time (sec) | 1.44 | 3.38
Table 4: Comparison of the performance between 4-th order Runge-Kutta ODE
method and Sinkhorn algorithm for martingale optimal transport problem.

###  5.4 Multi-period martingale optimal transport

We consider a 3-period martingale optimal transport problem. Let
$\mu=\sum_{r}\mu_{r}\delta_{x_{r}}\in\mathcal{P}(X),\theta=\sum_{s}\theta_{s}\delta_{y_{s}}\in\mathcal{P}(Y),\nu=\sum_{t}\nu_{t}\delta_{z_{t}}\in\mathcal{P}(Z)$
such that $\mu\preceq_{c}\theta\preceq_{c}\nu$. The linear constraints are:

$Q=MM=\\{g(x)(y-x)+h(x,y)(z-y)|g(x)\in C_{b}(X),h(x,y)\in C_{b}(X,Y)\\}.$

Note that

${\mathbb{E}}_{\pi}[Y|X]=X,\qquad{\mathbb{E}}_{\pi}[Z|X,Y]=Y\qquad\iff\pi\in\Pi^{MM}(\mu,\theta,\nu).$

With the cost tensor $c_{ijk}=c(x_{i},y_{j},z_{k})$, the entropic regularized
problem with fixed $\eta$ for the dual problem is

$\inf_{u,v,w,g,h}\Phi(u,v,w,g,h,\varepsilon)$

where

$\displaystyle\Phi(u,v,w,g,h,\varepsilon)=$
$\displaystyle-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\theta_{s}-\sum_{t}w_{t}\nu_{t}$
(5.23)
$\displaystyle\qquad+\eta\sum_{r,s,t}\exp\biggl{[}\frac{1}{\eta}\bigl{(}u_{r}+v_{s}+w_{t}+g_{r}(y_{s}-x_{r})$
$\displaystyle\qquad\qquad+h_{rs}(z_{t}-y_{s})-\varepsilon
c_{rst}\bigr{)}\biggr{]}\mu_{r}\theta_{s}\nu_{t}.$

To implement the Sinkhorn algorithm, we need the first order optimality
condition:

$\displaystyle u_{i}$
$\displaystyle=-\eta\log\left[\sum_{st}\exp\left[\frac{1}{\eta}(v_{s}+w_{t}+g_{i}(y_{s}-x_{i})+h_{is}(z_{t}-y_{s})-\varepsilon
c_{ist})\right]\theta_{s}\nu_{t}\right]$ (5.24) $\displaystyle v_{j}$
$\displaystyle=-\eta\log\left[\sum_{rt}\exp\left[\frac{1}{\eta}(u_{r}+w_{t}+g_{r}(y_{j}-x_{r})+h_{rj}(z_{t}-y_{j}))-\varepsilon
c_{rjt}\right]\mu_{r}\nu_{t}\right]$ $\displaystyle w_{k}$
$\displaystyle=-\eta\log\left[\sum_{rs}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}+g_{r}(y_{s}-x_{r})+h_{rs}(z_{k}-y_{s})-\varepsilon
c_{rsk})\right]\mu_{r}\theta_{s}\right]$ $\displaystyle 0$
$\displaystyle=\sum_{s,t}\exp\left[\frac{1}{\eta}(u_{i}+v_{s}+w_{t}+g_{i}(y_{s}-x_{i})+h_{is}(z_{t}-y_{s})-\varepsilon
c_{ist})\right](y_{s}-x_{i})\theta_{s}\nu_{t}$ $\displaystyle 0$
$\displaystyle=\sum_{t}\exp\left[\frac{1}{\eta}(u_{i}+v_{j}+w_{t}+g_{i}(y_{j}-x_{i})+h_{ij}(z_{t}-y_{j}))-\varepsilon
c_{ijt}\right](z_{t}-y_{j})\nu_{t}.$

The fitting of $u_{i},v_{j},w_{k}$ are straight forward while we use Newton’s
method to solve for $g_{i}$ and $h_{ij}$ in each of the iteration.

On the other hand, by combining (5.23) and (5.24), with a slight abuse of
notation again we set:

$\displaystyle\Phi(v,w,g,h,\varepsilon)=$
$\displaystyle-\sum_{s}v_{s}\theta_{s}-\sum_{t}w_{t}\nu_{t}+\eta\sum_{r}\log\biggl{[}\sum_{s,t}\exp\biggl{[}\frac{1}{\eta}\bigl{(}v_{s}+w_{t}$
(5.25) $\displaystyle\qquad+g_{r}(y_{s}-x_{r})+h_{rs}(z_{t}-y_{s})-\varepsilon
c_{rst}\bigr{)}\biggr{]}\theta_{s}\nu_{t}\biggr{]}\mu_{r}+\eta,$

and note that the minimizing
$v(\varepsilon),w(\varepsilon),g(\varepsilon),h(\varepsilon)$ solve the well-
posed ODE (3.5) after we fix $v_{0},w_{0},g_{0}$ and $h_{00}$ all to be 0. We
refer the reader to Appendix B.3 for the Hessian of $\Phi$ with respect to the
dual potential variables $v,w,g,h$ and the mixed derivatives with respect to
$v,w,g,h$ and $\varepsilon$.

For the numerical experiment, we took the marginals to be uniform on 30, 60
and 90 evenly spaced grid points on $X=[-0.1,0.1]$, $Y=[-0.4,0.4]$ and
$Z=[-1,1]$, respectively. We set $\eta=0.006$ and $\varepsilon\in[0,1]$ with
25 time steps. The cost function is $c(x,y,z)=y^{2}\exp(-x)+z^{2}\exp(-x)$
which is of martingale Spence-Mirrlees type. It is known that the optimal
measure is the multi-period left-monotone coupling [27]. Since $\mu$ is a
discrete measure, the multi-period left-monotone coupling is non-unique, but
we know from Lemma A.2, the optimal measures for regularized problems will
converge to the left-monotone coupling with minimal entropy. Explicitly, the
left-monotone condition uniquely determines the two fold marginals
$\gamma^{xy}=\gamma_{x}^{y}\otimes\mu$,
$\gamma^{xz}=\gamma_{x}^{z}\otimes\mu$, and so the only degrees of freedom are
in the conditional probabilities $\gamma_{x}^{yz}$ of
$\gamma=\gamma_{x}^{yz}\otimes\mu$, which must all be martingale couplings of
$\gamma_{x}^{y}$ and $\gamma_{x}^{z}$. We compute the marginal couplings of
these conditional probabilities withe minimal entropy and use these values in
Table 5.

In Figure 5, we compare the optimal measures calculated by the ODE and
Sinkhorn algorithms. The top row of sub-graphs represents 2D projections of
pairs of marginals obtained using the ODE method. The bottom row shows sub-
graphs generated using the Sinkhorn algorithm. Theoretically, both the $xy$
and $xz$ projections are left monotone couplings for these particular
marginals (at least in the unregularized limit, by the Lemma 6.5 of [27] and
the definition of multi-period left-monotone coupling), and so the conditional
probabilities $\gamma_{x}^{y}$ and $\gamma_{x}^{z}$ are supported on two
disjoint consecutive intervals for each $x$. Although the $xy$ projection
looks a little bit spread for the ODE generated measure, we encounter a
similar situation for the measure generated by Sinkhorn algorithm.
Nevertheless, as shown in the table 5, the values generated by both methods
still fall in the desired range.

Figure 5: 2D projections for pairs of marginals of the optimal measures for multi-period martingale transport problem. Top: generated by ODE method. Bottom: generated by Sinkhorn algorithm | ODE | Sinkhorn
---|---|---
Computed regularized primal value | 0.3807 | 0.3807
Optimal unregularized primal value | 0.3767 | 0.3767
Optimal unregularized primal value + Entropy | 0.4127 | 0.4127
Iterations | 25 | 60
CPU time (sec) | 8.71 | 7.40
Table 5: Comparison of the performance between 4-th order Runge-Kutta ODE
method and Sinkhorn algorithm for multi-period martingale transport problem.

###  5.5 Wasserstein geodesics and barycenters

A slight variant of the theory developed in this paper can be used to compute,
and provide a new perspective on, Wasserstein geodesics (also known as
displacement interpolants [24]) and barycenters [2]. Let us consider, firstly,
the computation of a geodesic (since we are considering a regularized problem,
it would be more appropriate to refer to it as the _entropic interpolant_)
between two measures $\mu^{1}$ and $\mu^{2}$. The geodesic at time
$\varepsilon$ can be obtained as the second marginal of a coupling
$\gamma\in{\mathcal{P}}(X^{1}\times Z\times X^{2})$ whose first and third
marginals are $\mu^{1}$ and $\mu^{2}$, respectively, and such that it
minimizes an optimal transport problem with cost

$c_{\varepsilon}(x^{1},z,x^{2})=(1-\varepsilon)|x^{1}-z|^{2}+\varepsilon|z-x^{2}|^{2}.$

After regularizing, we minimize

$\int_{X^{1}\times Z\times X^{2}}c_{\varepsilon}d\gamma+\eta
H_{\mu^{1}\otimes\mathcal{U}\otimes\mathcal{\mu}^{2}}(\gamma).$ 666Note that
here, since the $z$ marginal of $\gamma$ is not prescribed in the problem, we
take the second marginal of the reference measure to be uniform $\mathcal{U}$.

Note that our formulation of the geodesic problem actually falls slightly
outside the general ODE framework developed in Section 3 in two ways. First,
the problem is not a typical optimal transport problem, since the marginal of
the variable $z$ is not prescribed (this, the displacement interpolant, is in
fact what we are trying to find). Secondly, the cost takes the form
$c_{\varepsilon}=c_{0}+\varepsilon c$, rather than the form
$c_{\varepsilon}=\varepsilon c$ dealt with above. Nonetheless, the development
of the ODE (3.5) and proof of its well posedness follow almost exactly as
above, and so we can solve the ODE to approximate the displacement
interpolant. Notice that now the parameter $\varepsilon$ plays the role of the
time. In Figure 6 we plot the interpolant, between the uniform on $[0,1]$ and
the sum of two gaussians, computed by using the ODE approach.

|  |  |  |
---|---|---|---|---
Figure 6: Entropic interpolant at different time steps $\varepsilon$

In a very similar way one can extend this approach to the computation of
Wasserstein barycenter [2]: in this case we look for a coupling
$\gamma\in{\mathcal{P}}(\times_{i=1}^{n}X^{i}\times Z)$ having $n$ fixed
marginals equal to $\mu^{i}$ such that it minimizes the optimal transport
problem with cost

$c(x^{1},\ldots,x^{n},z)=\sum_{i=1}^{n}\lambda_{i}(\varepsilon)|x^{i}-z|^{2},$

where $\lambda_{i}(\varepsilon)$ are the usual weights such that
$\lambda_{i}(\varepsilon)\geq 0$ and $\sum_{i}\lambda_{i}(\varepsilon)=1$ for
every $\varepsilon$. In this case the ODE returns the barycenter for every
weights $\lambda_{i}(\varepsilon)$ as they vary in $\varepsilon$.

######  Remark 5.2 (Initial condition).

One difference with the theory we have developed above concerns the initial
condition for the Cauchy problem: here, one has to solve a simplified entropic
optimal transport problem where the cost is not zero. For instance, in the
geodesic case, one should solve a three marginal problem with cost
$|x^{1}-z|^{2}$. In particular, by the optimality conditions, the optimal dual
variable $u$ associated to the marginal $\mu^{1}$ has the following explicit
form

$u_{r}=-\eta\log\left[\sum_{s}Z\exp\left[\frac{1}{\eta}(-\varepsilon
c_{rs})\right]\right],$ (5.26)

where $c_{rs}$ is the cost $c_{\varepsilon}$ above evaluated at
$\varepsilon=0$. Notice that since the potential associated to the marginal
$\mu^{2}$ is a constant we can easily set it to be equal to $0$ (or,
equivalently, $Z=1$).

As for the barycenter problem, choose $\lambda_{i}$ such that at
$\varepsilon=0$ there is only one weight equals to $1$, so that the curve
starts at one of the marginal measures; say $\mu^{1}$. The cost function then
becomes $|x^{1}-z|^{2}$, and the initial potentials can be determined
similarly to above; that is, all potentials are equal to $0$ except for the
one corresponding to $x^{1}$, which is determined again by formula (5.26)above
with $Z=1$.

Acknowledgments. L.N. is partially on academic leave at Inria (team
Matherials) for the year 2023-2024 and acknowledges the hospitality if this
institution during this period. His work benefited from the support of the
FMJH Program PGMO, from H-Code, Université Paris-Saclay and from the ANR
project GOTA (ANR-23-CE46-0001). The work of B.P. was partially supported by
National Sciences and Engineering Research Council of Canada Discovery Grant
number 04658-2018. The work of J.H. was completed in partial fulfillment of
the requirements of a doctoral degree in applied mathematics at the University
of Alberta.

##  Appendix A Convergence of entropic optimal transport with extra linear
constraints to the unregularized limit

We aim to show that the optimal measure for the entropic optimal transport
problem with extra linear constraints will converge to an optimal solution of
the unregularized problem that has the minimum entropy relative to the product
measure.

######  Lemma A.1.

$\Pi^{Q}(\mu^{1},...,\mu^{n})$ is compact in the weak topology.

###### Proof.

Since $\Pi(\mu^{1},...,\mu^{n})$ is a weakly compact set, we only need to show
that $\Pi^{Q}$ is a closed subset of $\Pi$. If
$\gamma_{k}\to\gamma_{\infty}\in\Pi$ weakly, then

$\int_{X}qd\gamma_{\infty}=\lim_{k\to\infty}\int_{X}qd\gamma_{k}=0,\qquad\forall
q\in Q$

since $Q\subset C_{b}(X)$. Hence $\gamma_{\infty}\in\Pi^{Q}$. ∎

Given a lower semi-continuous cost $c$. Let $\eta_{k}\downarrow 0$ and
$\gamma_{\eta_{k}}$ solve:

$\inf_{\gamma\in\Pi^{Q}(\mu^{1},...,\mu^{n})}\int_{X}cd\gamma+\eta_{k}H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma).$
(A.1)

Note that for general marginals, it is often the case that all minimizers of
(A.1) have infinite entropy. To make sense of the statement that cluster
points minimize the entropy among solutions to the unregularized problem, we
therefore require an additional assumption (uniform boundedness of the entropy
in the Lemma below). This assumption is certainly satisfied in the discrete
case, which is the main focus of this paper.

######  Lemma A.2.

Suppose that $\eta_{k}$ is a sequence converging to $0$ such that the
minimizers $\gamma_{\eta_{k}}$ in (A.1) converge weakly to some
$\gamma_{0}=\lim\gamma_{\eta_{k}}$, and that the entropy of the
$\gamma_{\eta_{k}}$ is uniformly bounded, that is, $\exists M>0$ such that
$\forall k,H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma_{\eta_{k}})<M$. Then
$\gamma_{0}$ solves the unregularized optimal transport problem with
additional linear constraints (2.1). Moreover, $\gamma_{0}$ has minimal
entropy among all solutions of (2.1):

$\gamma_{0}\in argmin_{\gamma\text{ solves }\eqref{eqn: unregularized
primal}}H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma).$

###### Proof.

By the assumption of $\gamma_{\eta_{k}}$, for any measure $\gamma\in\Pi^{Q}$
such that $H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma)<\infty$, we have

$\int_{X}cd\gamma_{\eta_{k}}+\eta_{k}H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma_{\eta_{k}})\leq\int_{X}cd\gamma+\eta_{k}H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma).$

By taking the limit of both side for $k$, with the lower semi-continuity of
the cost and uniformly boundedness of the relative entropy of
$\gamma_{\eta_{k}}$, we thus have:

$\int_{X}cd\gamma_{0}\leq\int_{X}cd\gamma.$

Note that by Lemma A.1, $\gamma_{0}\in\Pi^{Q}$. Therefore $\gamma_{0}$ solve
the original optimal transport problem with extra linear constraints.

Now, suppose $\bar{\gamma}\in\Pi^{Q}$ is another optimal measure that solves
the unregularized problem with
$H_{\otimes_{i=1}^{n}\mu^{i}}(\bar{\gamma})<\infty$. Again by the optimality
assumption of $\gamma_{\eta_{k}}$,

$\int_{X}cd\gamma_{\eta_{k}}+\eta_{k}H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma_{\eta_{k}})\leq\int_{X}cd\bar{\gamma}+\eta_{k}H_{\otimes_{i=1}^{n}\mu^{i}}(\bar{\gamma}).$
(A.2)

The optimality assumption for $\bar{\gamma}$ to the unregularized problem
implies

$\int_{X}cd\gamma_{\eta_{k}}\geq\int_{X}cd\bar{\gamma}.$ (A.3)

Combining (A.2) and (A.3) we get

$H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma_{\eta_{k}})\leq
H_{\otimes_{i=1}^{n}\mu^{i}}(\bar{\gamma}).$

Taking the limit on both sides again and by the lower semi-continuity of the
relative entropy functional [32], we deduce that

$H_{\otimes_{i=1}^{n}\mu^{i}}(\gamma_{0})\leq
H_{\otimes_{i=1}^{n}\mu^{i}}(\bar{\gamma}).$

This means $\gamma_{0}$ has the minimum relative entropy among all optimal
measures. ∎

In the following, we prove the convergence of the dual potentials on discrete
domain. We adapt to the simplified notation to represent the entropic optimal
transport problem with extra linear constraints introduced in section 2.2:

$\min-b^{T}v+\eta\sum_{\ell=1}^{m}\exp\left(\frac{A_{\ell}v-\varepsilon
c_{\ell}}{\eta}\right)\boldsymbol{\mu}_{\ell},$ (A.4)

where the matrix $A$ has full column rank. We assume the set of solutions to
the dual unregularized optimal transport problem with extra linear constraints
is non-empty and bounded in the sense that if

$D_{0}=\inf-b^{T}v$ (A.5)

subject to $A_{\ell}v\leq c_{\ell}$ and

$V_{0}\coloneqq\left\\{v\in{\mathbb{R}}^{d}|A_{\ell}v\leq
c_{\ell},-b^{T}v=D_{0}\right\\},$

then there exist an $M>0$ such that $||v||<M,\forall v\in V_{0}$. Moreover we
assume that $\forall\ell,c_{\ell}\geq 0$. Let $v_{\eta}$ be the optimal
solution to the (A.4) with parameter $\eta$ and $\eta_{k}\downarrow 0$. We
denote $v_{k}=v_{\eta_{k}}$ for simplicity.

######  Proposition A.3.

The sequence $\\{v_{k}\\}_{k}$ is bounded. If it converges to some $v_{0}=\lim
v_{k}$, then $v_{0}$ is optimal in (A.5).

###### Proof.

Let $\bar{v}\in V_{0}$. Then for all $k$,

$-b^{T}v_{k}+\eta_{k}\sum_{\ell=1}^{m}\exp\left(\frac{A_{\ell}v_{k}-c_{\ell}}{\eta_{k}}\right)\boldsymbol{\mu}_{\ell}\leq-b^{T}\bar{v}+\eta_{k}\sum_{\ell=1}^{m}\exp\left(\frac{A_{\ell}\bar{v}-c_{\ell}}{\eta_{k}}\right)\boldsymbol{\mu}_{\ell}\leq-b^{T}\bar{v}+\eta_{k}\leq\bar{M}$
(A.6)

for some $\bar{M}>0$. The second inequality comes from the fact that
$A_{\ell}\bar{v}\leq c_{\ell}$ and hence the exponential term is less than 1.
Since exponential function and $\boldsymbol{\mu}_{\ell}$ are always positive,
from (A.6) we obtain for all $\ell$:

$A_{\ell}v_{k}\leq\eta_{k}\log\left[\frac{1}{\eta_{k}}(\bar{M}+b^{T}v_{k})\right]+c_{\ell}.$

If $\\{v_{k}\\}_{k}$ is unbounded and $||v_{k}||\to\infty$ as $k\to\infty$ and
without loss of generality if $v_{k}/||v_{k}||\to\hat{v}\neq 0$, then

$A_{\ell}\frac{v_{k}}{||v_{k}||}\leq\frac{\eta_{k}}{||v_{k}||}\log\left[\frac{1}{\eta_{k}}(\bar{M}+b^{T}v_{k})\right]+\frac{c_{\ell}}{||v_{k}||}\to
0.$

We can conclude that $A\hat{v}\leq 0$. Moreover, the positive of exponential
function and $\mu_{i}$ also implies that $-b^{T}\hat{u}\leq 0$ by the same
argument. Therefore, $\bar{u}+\lambda\hat{u}$ is also an optimal solution to
the unregularized problem which violate the boundedness of $V_{0}$.

Assuming $v_{k}\to v_{0}$ as $k\to\infty$. Since $\\{v_{k}\\}_{k}$ is bounded,
we can find another constant $\hat{M}>0$ such that

$\displaystyle\eta_{k}\sum_{\ell}\exp\left(\frac{A_{\ell}v_{k}-\varepsilon
c_{\ell}}{\eta_{k}}\right)\boldsymbol{\mu}_{\ell}$ $\displaystyle\leq\hat{M}$
$\displaystyle A_{\ell}v_{k}$
$\displaystyle\leq\eta_{k}\log\left(\frac{\hat{M}}{\eta_{k}}\right)+c_{\ell}.$

By taking $k\to\infty$ for both sides, we deduce that $A_{i}u_{0}\leq c_{i}$.
Again, by (A.6):

$-b^{T}v_{0}\leq\limsup_{k\to\infty}-b^{T}v_{k}+\eta_{k}\sum_{\ell=1}^{m}\exp\left(\frac{A_{\ell}v_{k}-c_{\ell}}{\eta_{k}}\right)\boldsymbol{\mu}_{\ell}\leq-b^{T}\bar{v}.$

Therefore, we can conclude that $v_{0}$ is also an optimal solution to the
unregularized problem. ∎

##  Appendix B Calculation of the second derivatives of $C(\varepsilon)$

Recall:

$C^{\prime\prime}(\varepsilon)=\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]^{T}\left[-D^{2}_{uu}\Phi(u(\varepsilon),\varepsilon)\right]^{-1}\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]+\frac{{\partial}^{2}}{{\partial}\varepsilon^{2}}\Phi(u(\varepsilon),\varepsilon).$
(B.1)

Let

$e_{rs}=\exp\left[\frac{1}{\eta}(u_{r}-\varepsilon c_{rs})\right]\mu_{r}$

and $e_{rs}^{c}=e_{rs}c_{rs},e_{rs}^{c2}=e_{rs}c_{rs}^{2}$. It is easy to
deduce all the second derivative of $\Phi$:

$\displaystyle\frac{{\partial}^{2}}{{\partial}\varepsilon^{2}}\Phi(u(\varepsilon),\varepsilon)$
$\displaystyle=\frac{1}{\eta}\sum_{s}\left[\frac{\sum_{r}e_{rs}^{c2}}{\sum_{r}e_{rs}}-\frac{\sum_{r}e_{rs}^{c}\cdot\sum_{r}e_{rs}^{c}}{(\sum_{r}e_{rs})^{2}}\right]\nu_{s}$
$\displaystyle\left(\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right)_{i}$
$\displaystyle=-\frac{1}{\eta}\sum_{s}\left[\frac{e_{is}^{c}}{\sum_{r}e_{rs}}-\frac{e_{is}\cdot\sum_{r}e_{rs}^{c}}{(\sum_{r}e_{rs})^{2}}\right]\
\nu_{s}$
$\displaystyle\biggl{(}D^{2}_{uu}\Phi(u(\varepsilon),\varepsilon)\biggr{)}_{i_{1}i_{2}}$
$\displaystyle=\frac{1}{\eta}\sum_{s}\left[\frac{e_{i_{1}s}}{\sum_{r}e_{rs}}\delta_{i_{1}i_{2}}-\frac{e_{i_{1}s}\cdot
e_{i_{2}s}}{(\sum_{r}e_{rs})^{2}}\right]\nu_{s},$

where $\delta_{ip}$ is the Kronecker delta.

Let us calculate the (B.1) separately. For the second part, since
$e_{ij}=\mu_{i}$ when $\varepsilon=0$, then:

$\frac{{\partial}^{2}}{{\partial}\varepsilon^{2}}\Phi(u(0),0)=\frac{1}{\eta}\left({\mathbb{E}}[c^{2}(X,Y)]-{\mathbb{E}}\left[({\mathbb{E}}[c(X,Y)|Y])^{2}\right]\right).$
(B.2)

Let’s deal with the first part of (B.1). With the assumption that $u_{0}\equiv
0$, then when $\varepsilon=0$, for every $i\neq 0$

$\left(\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(0),0)\right)_{i}=\frac{\mu_{i}}{\eta}\bigl{(}{\mathbb{E}}[c(X,Y)]-{\mathbb{E}}[c(X,Y)|X=x_{i}]\bigr{)}$

and

$-D^{2}_{uu}\Phi(u(0),0)=\frac{1}{\eta}\bigl{(}-\text{diag}(\bar{\mu})+\bar{\mu}{\bar{\mu}}^{T}\bigr{)},$

where $\bar{\mu}=\\{\mu_{i}\\}_{i\neq 0}$. Note that, by Sherman-Morrison
formula, we have:

$\displaystyle\left[-D^{2}_{uu}\Phi(u(0),0)\right]^{-1}$
$\displaystyle=\eta\left(-\text{diag}^{-1}(\bar{\mu})-\frac{\text{diag}^{-1}(\bar{\mu})\bar{\mu}{\bar{\mu}}^{T}\text{diag}^{-1}(\bar{\mu})}{1-{\bar{\mu}}^{T}\text{diag}^{-1}(\bar{\mu})\bar{\mu}}\right)$
$\displaystyle=\eta\left(-\text{diag}^{-1}(\bar{\mu})-\frac{\mathbbm{1}\mathbbm{1}^{T}}{\mu_{0}}\right),$

where $\mathbbm{1}$ is a vector with all entries 1 and
$\text{diag}^{-1}(\bar{\mu})$ is a diagonal matrix where each diagonal entry
is the inverse of $\mu_{i}$. Therefore

$\displaystyle\left[-D^{2}_{uu}\Phi(u(0),0)\right]^{-1}\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(0),0)\right]$
$\displaystyle=$
$\displaystyle\eta\left(-\text{diat}^{-1}(\bar{\mu})-\frac{\mathbbm{1}\mathbbm{1}^{T}}{\mu_{0}}\right)\left(\frac{\mu_{i}}{\eta}\bigl{(}{\mathbb{E}}[c(X,Y)]-{\mathbb{E}}[c(X,Y)|X=x_{i}]\bigr{)}\right)_{i\neq
0}$ $\displaystyle=$
$\displaystyle\biggl{(}{\mathbb{E}}[c(X,Y)|X=x_{i}]-{\mathbb{E}}[c(X,Y)]\biggr{)}_{i\neq
0}$ $\displaystyle\qquad\qquad+\left(\frac{1}{\mu_{0}}\sum_{n\neq
0}\left[{\mathbb{E}}[c(X,Y)|X=x_{i}]-{\mathbb{E}}[c(X,Y)]\right]\right)_{i\neq
0}$ $\displaystyle=$
$\displaystyle\bigl{(}{\mathbb{E}}[c(X,Y)|X=x_{i}]-{\mathbb{E}}[c(X,Y)]+{\mathbb{E}}[c(X,Y)]-{\mathbb{E}}[c(X,Y)|X=x_{0}]\bigr{)}_{i\neq
0}$ $\displaystyle=$
$\displaystyle\bigl{(}{\mathbb{E}}[c(X,Y)|X=x_{i}]-{\mathbb{E}}[c(X,Y)|X=x_{0}]\bigr{)}_{i\neq
0}$

and

$\displaystyle\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]^{T}\left[-D^{2}_{uu}\Phi(u(\varepsilon),\varepsilon)\right]^{-1}\left[\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),\varepsilon)\right]$
$\displaystyle=$
$\displaystyle\left(\frac{\mu_{i}}{\eta}\bigl{(}{\mathbb{E}}[c(X,Y)]-{\mathbb{E}}[c(X,Y)|X=x_{i}]\bigr{)}\right)_{i\neq
0}^{T}$
$\displaystyle\qquad\qquad\cdot\bigl{(}{\mathbb{E}}[c(X,Y)|X=x_{i}]-{\mathbb{E}}[c(X,Y)|X=x_{0}]\bigr{)}_{i\neq
0}$ $\displaystyle=$ $\displaystyle\frac{1}{\eta}\sum_{i\neq
0}\bigl{[}{\mathbb{E}}[c(X,Y)]{\mathbb{E}}[c(X,Y)|X=x_{i}]-{\mathbb{E}}[c(X,Y)]{\mathbb{E}}[c(X,Y)|X=x_{0}]$
$\displaystyle\qquad-\left({\mathbb{E}}[c(X,Y)|X=x_{i}]\right)^{2}+{\mathbb{E}}[c(X,Y)|X=x_{0}]{\mathbb{E}}[c(X,Y)|X=x_{i}]\bigr{]}$
$\displaystyle=$
$\displaystyle\frac{1}{\eta}\biggl{\\{}{\mathbb{E}}[c(X,Y)]\sum_{i}{\mathbb{E}}[c(X,Y)|X=x_{i}]\mu_{i}-{\mathbb{E}}[c(X,Y)]\sum_{i}{\mathbb{E}}[c(X,Y)|X=x_{0}]\mu_{i}$
$\displaystyle\qquad-\sum_{i}\left({\mathbb{E}}[c(X,Y)|X=x_{i}]\right)^{2}\mu_{i}+{\mathbb{E}}[c(X,Y)|X=x_{0}]\sum_{i}{\mathbb{E}}[c(X,Y)|X=x_{i}]\mu_{i}\biggr{\\}}$
$\displaystyle\qquad-\frac{1}{\eta}\biggl{\\{}{\mathbb{E}}[c(X,Y)]{\mathbb{E}}[c(X,Y)|X=x_{0}]-{\mathbb{E}}[c(X,Y)]{\mathbb{E}}[c(X,Y)|X=x_{0}]$
$\displaystyle\qquad\qquad-\left({\mathbb{E}}[c(X,Y)|X=x_{0}]\right)^{2}+{\mathbb{E}}[c(X,Y)|X=x_{0}]{\mathbb{E}}[c(X,Y)|X=x_{0}]\biggr{\\}}\mu_{0}$
$\displaystyle=$
$\displaystyle\frac{1}{\eta}\biggl{(}\left({\mathbb{E}}[c(X,Y)]\right)^{2}-{\mathbb{E}}\left[({\mathbb{E}}[c(X,Y)|X])^{2}\right]\biggr{)}.$
(B.3)

Hence, summing (B.2) and (B.3) we get:

$C^{\prime\prime}(0)=\frac{1}{\eta}\biggl{(}\left({\mathbb{E}}[c(X,Y)]\right)^{2}+{\mathbb{E}}[c^{2}(X,Y)]-{\mathbb{E}}\left[({\mathbb{E}}[c(X,Y)|X])^{2}\right]-{\mathbb{E}}\left[({\mathbb{E}}[c(X,Y)|Y])^{2}\right]\biggr{)}.$

###  B.1 The second order derivatives of $\Phi(u,v,\varepsilon)$ for multi
marginal optimal transport

Recall:

$\Phi(u,v,\varepsilon)=-\sum_{r}u_{r}\mu_{r}-\sum_{s}v_{s}\theta_{s}+\eta\sum_{t}\log\left[\sum_{r,s}\exp\left[\frac{1}{\eta}(u_{r}+v_{s}-\varepsilon
c_{rst})\right]\mu_{r}\theta_{s}\right]\nu_{rst}-\eta.$

Let us denote

$e_{ijk}=\exp\left[\frac{1}{\eta}(u_{i}+v_{j}-\varepsilon
c_{ijk})\right],\qquad
e_{ijk}^{c}=\exp\left[\frac{1}{\eta}(u_{i}+v_{j}-\varepsilon
c_{ijk})\right]c_{ijk}.$

The the Hessian matrix of $\Phi$ will be the symmetric matrix:

$D^{2}\Phi=\begin{bmatrix}D^{2}_{uu}\Phi&D^{2}_{uv}\Phi\\\
D^{2}_{vu}\Phi&D^{2}_{vv}\Phi\\\ \end{bmatrix},$

where:

$\displaystyle\biggl{(}D^{2}_{uu}\Phi(u(\varepsilon),v(\varepsilon),\varepsilon)\biggr{)}_{i_{1}i_{2}}$
$\displaystyle=\frac{1}{\eta}\sum_{t}\left[\frac{\sum_{s}e_{i_{1}st}}{\sum_{r,s}e_{rst}}\delta_{i_{1}i_{2}}-\frac{\sum_{s}e_{i_{1}st}\cdot\sum_{s}e_{i_{2}st}}{(\sum_{r,s}e_{rst})^{2}}\right]\nu_{t}$
$\displaystyle\biggl{(}D^{2}_{vv}\Phi(u(\varepsilon),v(\varepsilon),\varepsilon)\biggr{)}_{j_{1}j_{2}}$
$\displaystyle=\frac{1}{\eta}\sum_{t}\left[\frac{\sum_{r}e_{rjt}}{\sum_{r,s}e_{rst}}\delta_{j_{1}j_{2}}-\frac{\sum_{r}e_{rj_{1}t}\cdot\sum_{r}e_{rj_{2}t}}{(\sum_{r,s}e_{rst})^{2}}\right]\nu_{t}$
$\displaystyle\biggl{(}D^{2}_{uv}\Phi(u(\varepsilon),v(\varepsilon),\varepsilon)\biggr{)}_{ij}$
$\displaystyle=\frac{1}{\eta}\sum_{t}\left[\frac{e_{ijt}}{\sum_{r,s}e_{rst}}-\frac{\sum_{s}e_{ist}\cdot\sum_{r}e_{rjt}}{(\sum_{r,s}e_{rst})^{2}}\right]\nu_{t}$

and the mixed derivative is:

$\nabla\frac{{\partial}}{{\partial}\varepsilon}\Phi=\begin{bmatrix}\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\
\nabla_{v}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\ \end{bmatrix},$
(B.4)

where:

$\displaystyle\left(\nabla_{u}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),v(\varepsilon),\varepsilon)\right)_{i}$
$\displaystyle=-\frac{1}{\eta}\sum_{t}\left[\frac{\sum_{s}e_{ist}^{c}}{\sum_{r,s}e_{rst}}-\frac{\sum_{s}e_{ist}\cdot\sum_{r,s}e_{rst}^{c}}{(\sum_{r,s}e_{rst})^{2}}\right]\nu_{t}$
$\displaystyle\left(\nabla_{v}\frac{{\partial}}{{\partial}\varepsilon}\Phi(u(\varepsilon),v(\varepsilon),\varepsilon)\right)_{j}$
$\displaystyle=-\frac{1}{\eta}\sum_{t}\biggl{[}\frac{\sum_{r}e_{rjt}^{c}}{\sum_{r,s}e_{rst}}-\frac{\sum_{r}e_{rjt}\cdot\sum_{r,s}e_{rst}^{c}}{(\sum_{r,s}e_{rst})^{2}}\biggr{]}\nu_{\ell}.$

###  B.2 The second order derivatives of $\Phi(v,g,\varepsilon)$ for
martingale optimal transport

Recall:

$\Phi(v,g,\varepsilon)\coloneqq-\sum_{m}v_{m}\nu_{m}+\eta\sum_{n}\log\left[\sum_{m}\exp\left[\frac{1}{\eta}(v_{m}+g_{n}(y_{m}-x_{n})-\varepsilon
c_{nm})\right]\nu_{m}\right]\mu_{n}-\eta.$

Let denote

$e_{ij}=\exp\left[\frac{1}{\eta}(v_{j}+h_{i}(y_{j}-x_{i})-\varepsilon
c_{ij})\right]\nu_{j}$

and
$e_{ij}^{xy}=e_{ij}(y_{j}-x_{i}),e_{ij}^{xy2}=e_{ij}(y_{j}-x_{i})^{2},e_{ij}^{c}=e_{ij}c_{ij}$
and $e_{ij}^{xyc}=e_{ij}(y_{j}-x_{i})c_{ij}$

The Hessian matrix for $\Phi$ will be the symmetric matrix:

$D^{2}\Phi=\begin{bmatrix}D^{2}_{vv}\Phi&D^{2}_{vh}\Phi\\\
D^{2}_{hv}\Phi&D^{2}_{hh}\Phi\\\ \end{bmatrix},$

where:

$\displaystyle\biggl{(}D^{2}_{vv}\Phi(v(\varepsilon),h(\varepsilon),\varepsilon)\biggr{)}_{j_{1}j_{2}}$
$\displaystyle=\frac{1}{\eta}\sum_{r}\left[\frac{e_{rj}}{\sum_{s}e_{rs}}\delta_{j_{1}j_{2}}-\frac{e_{rj_{1}}\cdot
e_{rj_{2}}}{(\sum_{s}e_{rs})^{2}}\right]\mu_{r}$
$\displaystyle\biggl{(}D^{2}_{hh}\Phi(v(\varepsilon),h(\varepsilon),\varepsilon)\biggr{)}_{i_{1}i_{2}}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{s}e_{i_{1}s}^{xy2}}{\sum_{s}e_{i_{1}s}}-\frac{\sum_{s}e_{i_{1}s}^{xy}\cdot\sum_{s}e_{i_{2}s}^{xy}}{(\sum_{s}e_{i_{1}s})^{2}}\right]\mu_{i_{1}}\delta_{i_{1}i_{2}}$
$\displaystyle\biggl{(}D^{2}_{vh}\Phi(v(\varepsilon),h(\varepsilon),\varepsilon)\biggr{)}_{ij}$
$\displaystyle=\frac{1}{\eta}\left[\frac{e_{ij}}{\sum_{s}e_{is}}-\frac{e_{ij}\cdot\sum_{s}e_{is}^{xy}}{(\sum_{s}e_{is})^{2}}\right]\mu_{i},$

where $\delta_{i\ell}$ is the Kronecker delta. While the mixed derivative:

$\nabla\frac{{\partial}}{{\partial}\varepsilon}\Phi=\begin{bmatrix}\nabla_{v}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\
\nabla_{h}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\ \end{bmatrix},$

where:

$\displaystyle\left(\nabla_{v}\frac{{\partial}}{{\partial}\varepsilon}\Phi(v(\varepsilon),h(\varepsilon),\varepsilon)\right)_{j}$
$\displaystyle=-\frac{1}{\eta}\sum_{r}\left[\frac{e_{rj}^{c}}{\sum_{s}e_{rs}}-\frac{e_{rj}\cdot\sum_{s}e_{rs}^{c}}{(\sum_{s}e_{rs})^{2}}\right]\
\nu_{j}$
$\displaystyle\left(\nabla_{h}\frac{{\partial}}{{\partial}\varepsilon}\Phi(v(\varepsilon),h(\varepsilon),\varepsilon)\right)_{i}$
$\displaystyle=-\frac{1}{\eta}\biggl{[}\frac{\sum_{s}e_{is}^{xyc}}{\sum_{s}e_{is}}-\frac{\sum_{s}e_{is}^{xy}\cdot\sum_{s}e_{is}^{c}}{(\sum_{s}e_{is})^{2}}\biggr{]}\mu_{i}.$

###  B.3 The second order derivatives of $\Phi(v,w,g,h,\varepsilon)$ for
multi-period optimal transport

Recall:

$\displaystyle\Phi(v,w,g,h,\varepsilon)=$
$\displaystyle-\sum_{s}v_{s}\theta_{s}-\sum_{t}w_{t}\nu_{t}+\eta\sum_{r}\log\biggl{[}\sum_{s,t}\exp\biggl{[}\frac{1}{\eta}\bigl{(}v_{s}+w_{t}$
$\displaystyle\qquad+g_{r}(y_{s}-x_{r})+h_{rs}(z_{t}-y_{s})-\varepsilon
c_{rst}\bigr{)}\biggr{]}\theta_{s}\nu_{t}\biggr{]}\nu_{r}-\eta.$

By denoting
$e_{ijk}=\exp[\frac{1}{\eta}(v_{j}+w_{k}+g_{i}(y_{j}-x_{i})+h_{ij}(z_{k}-y_{j})-\varepsilon
c_{ijk})]\theta_{j}\nu_{k}$ and
$e_{ijk}^{xy}=e_{ijk}(y_{j}-x_{i}),e_{ijk}^{yz}=e_{ijk}(z_{k}-y_{j}),e_{ijk}^{xy2}=e_{ijk}(y_{j}-x_{i})^{2},e_{ijk}^{xyz}=e_{ijk}(y_{j}-x_{i})(z_{k}-y_{j}),e_{ijk}^{yz2}=e_{ijk}(z_{k}-y_{j})^{2},e_{ijk}^{c}=e_{ijk}c_{ijk},e_{ijk}^{xyc}=e_{ijk}(y_{j}-x_{i})c_{ijk},e_{ijk}^{yzc}=e_{ijk}(z_{k}-y_{j})c_{ijk},$

the Hessian of $\Phi$ will be the symmetric matrix:

$D^{2}\Phi=\begin{bmatrix}D^{2}_{vv}\Phi&D^{2}_{vw}\Phi&D^{2}_{vg}\Phi&D^{2}_{vh}\Phi\\\
D^{2}_{wv}\Phi&D^{2}_{ww}\Phi&D^{2}_{wg}\Phi&D^{2}_{wh}\Phi\\\
D^{2}_{gv}\Phi&D^{2}_{gw}\Phi&D^{2}_{gg}\Phi&D^{2}_{gh}\Phi\\\
D^{2}_{hv}\Phi&D^{2}_{hw}\Phi&D^{2}_{hg}\Phi&D^{2}_{hh}\Phi\\\ \end{bmatrix},$

where

$\displaystyle\left(D^{2}_{vv}\Phi\right)_{j_{1}j_{2}}$
$\displaystyle=\frac{1}{\eta}\sum_{r}\left[\frac{\sum_{t}e_{rj_{1}t}}{\sum_{s,t}e_{rst}}\delta_{j_{1}j_{2}}-\frac{\sum_{t}e_{rj_{1}t}\cdot\sum_{t}e_{rj_{2}t}}{(\sum_{s,t}e_{rst})^{2}}\right]\mu_{r}$
$\displaystyle\left(D^{2}_{vw}\Phi\right)_{jk}$
$\displaystyle=\frac{1}{\eta}\sum_{r}\left[\frac{e_{rjk}}{\sum_{s,t}e_{rst}}-\frac{\sum_{t}e_{rjt}\cdot\sum_{s}e_{rsk}}{(\sum_{s,t}e_{rst})^{2}}\right]\mu_{r}$
$\displaystyle\left(D^{2}_{vg}\Phi\right)_{ji}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{t}e_{ijt}}{\sum_{s,t}e_{ist}}-\frac{\sum_{t}e_{ijt}\cdot\sum_{s,t}e_{ist}^{xy}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i}$
$\displaystyle\left(D^{2}_{vh}\Phi\right)_{j_{1},(ij_{2})}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{t}e_{ij_{2}}^{yz}}{\sum_{s,t}e_{ist}}\delta_{j_{1}j_{2}}-\frac{\sum_{t}e_{ij_{1}}\cdot\sum_{t}e_{ij_{2}}^{yz}}{(\sum_{s,t}e_{ist})^{2}}\right]\nu_{i}$
$\displaystyle\left(D^{2}_{ww}\Phi\right)_{k_{1}k_{2}}$
$\displaystyle=\frac{1}{\eta}\sum_{r}\left[\frac{\sum_{s}e_{rk_{1}}}{\sum_{s,t}e_{rst}}\delta_{k_{1}k_{2}}-\frac{\sum_{s}e_{rsk_{1}}\cdot\sum_{s}e_{rsk_{2}}}{(\sum_{s,t}e_{rst})^{2}}\right]\mu_{r}$
$\displaystyle\left(D^{2}_{wg}\Phi\right)_{ki}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{s}e_{isk}^{xy}}{\sum_{s,t}e_{ist}}-\frac{\sum_{s}e_{isk}\cdot\sum_{s,t}e_{ist}^{xy}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i}$
$\displaystyle\left(D^{2}_{wh}\Phi\right)_{k,(ij)}$
$\displaystyle=\frac{1}{\eta}\left[\frac{e_{ijk}^{yz}}{\sum_{s,t}e_{ist}}-\frac{\sum_{s}e_{isk}\cdot\sum_{t}e_{ijt}^{yz}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i}$
$\displaystyle\left(D^{2}_{gg}\Phi\right)_{i_{1}i_{2}}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{s,t}e_{i_{1}st}^{xy2}}{\sum_{s,t}e_{i_{1}st}}-\frac{\sum_{s,t}e_{i_{1}st}^{xy}\cdot\sum_{s,t}e_{i_{2}st}^{xy}}{(\sum_{s,t}e_{i_{1}st})^{2}}\right]\mu_{i_{1}}\delta_{i_{1}i_{2}}$
$\displaystyle\left(D^{2}_{gh}\Phi\right)_{i_{1},(i_{2}j)}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{t}e_{i_{1}jt}^{xyz}}{\sum_{s,t}e_{ist}}-\frac{\sum_{s,t}e_{i_{1}st}^{xy}\cdot\sum_{t}e_{i_{2}jt}^{yz}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i_{1}}\delta_{i_{1}i_{2}}$
$\displaystyle\left(D^{2}_{hh}\Phi\right)_{(i_{1}j_{1}),(i_{2}j_{2})}$
$\displaystyle=\frac{1}{\eta}\left[\frac{\sum_{t}e_{i_{1}j_{1}t}^{yz2}}{\sum_{s,t}e_{i_{1}st}}\delta_{j_{1}j_{2}}-\frac{\sum_{t}e_{i_{1}j_{1}t}^{yz}\cdot\sum_{t}e_{i_{2}j_{2}t}^{yz}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i_{1}}\delta_{i_{1}i_{2}}.$

Here we flatten the matrix $\\{h_{ij}\\}$ into a long vector and abuse the
notation that we still call it with the index $(ij)$. For the mixed second
order derivatives:

$\nabla\frac{{\partial}}{{\partial}\varepsilon}\Phi=\begin{bmatrix}\nabla_{v}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\
\nabla_{w}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\
\nabla_{g}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\
\nabla_{h}\frac{{\partial}}{{\partial}\varepsilon}\Phi\\\ \end{bmatrix},$

where

$\displaystyle\left(\nabla_{v}\frac{{\partial}}{{\partial}\varepsilon}\Phi\right)_{j}$
$\displaystyle=-\frac{1}{\eta}\sum_{r}\left[\frac{\sum_{t}e_{rjt}^{c}}{\sum_{s,t}e_{rst}}-\frac{\sum_{t}e_{rjt}\cdot\sum_{s,t}e_{rst}^{c}}{(\sum_{s,t}e_{rst})^{2}}\right]\mu_{r}$
$\displaystyle\left(\nabla_{g}\frac{{\partial}}{{\partial}\varepsilon}\Phi\right)_{i}$
$\displaystyle=-\frac{1}{\eta}\left[\frac{\sum_{s,t}e_{ist}^{xyc}}{\sum_{s,t}e_{ist}}-\frac{\sum_{s,t}e_{ist}^{xy}\cdot\sum_{s,t}e_{ist}^{c}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i}$
$\displaystyle\left(\nabla_{w}\frac{{\partial}}{{\partial}\varepsilon}\Phi\right)_{k}$
$\displaystyle=-\frac{1}{\eta}\sum_{r}\left[\frac{\sum_{s}e_{rsk}^{c}}{\sum_{s,t}e_{rst}}-\frac{\sum_{s}e_{rsk}\cdot\sum_{st}e_{rst}^{c}}{(\sum_{s,t}e_{rst})^{2}}\right]\mu_{r}$
$\displaystyle\left(\nabla_{h}\frac{{\partial}}{{\partial}\varepsilon}\Phi\right)_{(ij)}$
$\displaystyle=-\frac{1}{\eta}\left[\frac{\sum_{t}e_{ijt}^{yzc}}{\sum_{s,t}e_{ist}}-\frac{\sum_{t}e_{ijt}^{yz}\cdot\sum_{s,t}e_{ist}^{c}}{(\sum_{s,t}e_{ist})^{2}}\right]\mu_{i}.$

## References

  * [1] Beatrice Acciaio, Mathias Beiglböck, and Gudmund Pammer.  Weak transport for non-convex costs and model-independence in a fixed-income market.  Mathematical Finance, 31(4):1423–1453, 2021.
  * [2] Martial Agueh and Guillaume Carlier.  Barycenters in the Wasserstein space.  SIAM Journal on Mathematical Analysis, 43(2):904–924, 2011.
  * [3] Julio Backhoff-Veraguas, Daniel Bartl, Mathias Beiglböck, and Manu Eder.  All adapted topologies are equal.  Probability Theory and Related Fields, 178:1125–1172, 2020.
  * [4] Mathias Beiglböck and Claus Griessler.  A land of monotone plenty.  Annali Scuola Normale Superiore-classe Di Scienze, pages 109–127, 2019.
  * [5] Mathias Beiglböck, Benjamin Jourdain, William Margheriti, and Gudmund Pammer.  Stability of the weak martingale optimal transport problem.  The Annals of Applied Probability, 33(6B):5382–5412, 2023.
  * [6] Mathias Beiglböck and Nicolas Juillet.  On a problem of optimal transport under marginal martingale constraints.  The Annals of Probability, 44(1):42 – 106, 2016.
  * [7] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré.  Iterative Bregman projections for regularized transportation problems.  SIAM J. Sci. Comput., 37(2):A1111–A1138, 2015.
  * [8] Espen Bernton, Promit Ghosal, and Marcel Nutz.  Entropic optimal transport: geometry and large deviations.  Duke Math. J., 171(16):3363–3400, 2022.
  * [9] Giuseppe Buttazzo, Luigi De Pascale, and Paola Gori-Giorgi.  Optimal-transport formulation of electronic density-functional theory.  Physical Review A, 85(6):062502, 2012.
  * [10] Guillaume Carlier.  On the linear convergence of the multimarginal Sinkhorn algorithm.  SIAM Journal on Optimization, 32(2):786–794, 2022.
  * [11] Guillaume Carlier, Paul Pegon, and Luca Tamanini.  Convergence rate of general entropic optimal transport costs.  Calc. Var. Partial Differential Equations, 62(4):Paper No. 116, 28, 2023.
  * [12] Maria Colombo, Luigi De Pascale, and Simone Di Marino.  Multimarginal optimal transport maps for one-dimensional repulsive costs.  Canadian Journal of Mathematics, 67(2):350–368, 2015.
  * [13] R Cominetti and J. San Martín.  Asymptotic analysis of the exponential penalty trajectory in linear programming.  Mathematical Programming, 67:169–187, 1994.
  * [14] Giovanni Conforti and Luca Tamanini.  A formula for the time derivative of the entropic cost and applications.  Journal of Functional Analysis, 280(11):108964, 2021.
  * [15] Codina Cotar, Gero Friesecke, and Claudia Klüppelberg.  Density functional theory and optimal transportation with Coulomb cost.  Communications on Pure and Applied Mathematics, 66(4):548–599, 2013.
  * [16] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  Advances in neural information processing systems, 26, 2013.
  * [17] Hadrien De March.  Entropic approximation for multi-dimensional martingale optimal transport.  arXiv preprint arXiv:1812.11104, 2018.
  * [18] Simone Di Marino, Augusto Gerolin, and Luca Nenna.  Optimal transportation theory with repulsive costs.  In Topological optimization and optimal transport, volume 17 of Radon Ser. Comput. Appl. Math., pages 204–256. De Gruyter, Berlin, 2017.
  * [19] Stephan Eckstein and Gudmund Pammer.  Computational methods for adapted optimal transport.  The Annals of Applied Probability, 34(1A):675–713, 2024.
  * [20] Joshua Zoen-Git Hiew, Tongseok Lim, Brendan Pass, and Marcelo Cruz de Souza.  Geometry of vectorial martingale optimal transport and robust option pricing.  arXiv preprint arXiv:2309.04947, 2023.
  * [21] Rémi Lassalle.  Causal transport plans and their Monge–Kantorovich problems.  Stochastic Analysis and Applications, 36(3):452–484, 2018.
  * [22] Rodrigue Lelotte.  Sur les gaz de Coulomb et le transport optimal multimarge.  Theses, Université Paris sciences et lettres, October 2023.
  * [23] Christian Léonard.  A survey of the Schrödinger problem and some of its connections with optimal transport.  Discrete and Continuous Dynamical Systems, 34(4):1533–1574, 2014.
  * [24] Robert J. McCann.  A convexity principle for interacting gases.  Adv. Math., 128(1):153–179, 1997.
  * [25] Luca Nenna and Brendan Pass.  An ODE characterisation of multi-marginal optimal transport.  arXiv preprint arXiv:2212.12492, 2022.
  * [26] Luca Nenna and Paul Pegon.  Convergence rate of entropy-regularized multi-marginal optimal transport costs.  Canadian Journal of Mathematics, page 1–22, 2024.
  * [27] Marcel Nutz, Florian Stebegg, and Xiaowei Tan.  Multiperiod martingale transport.  Stochastic Processes and their Applications, 130(3):1568–1615, 2020.
  * [28] Marcel Nutz and Johannes Wiesel.  Entropic optimal transport: Convergence of potentials.  Probability Theory and Related Fields, 184(1):401–424, 2022.
  * [29] J. M. Ortega and W. C. Rheinboldt.  Iterative solution of nonlinear equations in several variables.  Society for Industrial and Applied Mathematics, 2000.
  * [30] Soumik Pal.  On the difference between entropic cost and the optimal transport cost.  Ann. Appl. Probab., 34(1B):1003–1028, 2024.
  * [31] Brendan Pass.  Multi-marginal optimal transport: theory and applications.  ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique, 49(6):1771–1790, 2015.
  * [32] E. Posner.  Random coding strategies for minimum entropy.  IEEE Transactions on Information Theory, 21(4):388–391, 1975.
  * [33] Filippo Santambrogio.  Optimal transport for applied mathematicians, volume 87 of Progress in nonlinear differential equations and their applications.  Birkhäuser/Springer, Cham, 2015.  Calculus of variations, PDEs, and modeling.
  * [34] Cédric Villani.  Topics in optimal transportation, volume 58 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2003.
  * [35] Cédric Villani.  Optimal transport: old and new, volume 338.  Springer, 2009.
  * [36] Danila A Zaev.  On the Monge-Kantorovich problem with additional linear constraints.  Mathematical Notes, 98:725–741, 2015.