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

Joshua Hiew, Luca Nenna, Brendan Pass

► To cite this version:

Joshua Hiew, Luca Nenna, Brendan Pass. An ordinary differential equation for entropic optimal transport and its linearly constrained variants. 2024. hal-04525602

HAL Id: hal-04525602

https://inria.hal.science/hal-04525602

Preprint submitted on 28 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An ordinary differential equation for entropic optimal transport and its linearly constrained variants

Joshua Hiew, Luca Nenna, Brendan Pass

► To cite this version:

Joshua Hiew, Luca Nenna, Brendan Pass. An ordinary differential equation for entropic optimal transport and its linearly constrained variants. 2024. hal-04525602

HAL Id: hal-04525602

https://inria.hal.science/hal-04525602v1

Preprint submitted on 28 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An ordinary differential equation for entropic optimal transport and its linearly constrained variants

Joshua Zoen-Git Hiew,\* Luca Nenna† and Brendan Pass ‡

March 28, 2024

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 \to \mathbb{R}$, the Monge-Kantorovich optimal transport problem is to

<footnote>\*Department of Mathematical and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta, Canada, T6G 2G1. email: joshuazo@ualberta.ca</footnote>

<footnote>†Université Paris-Saclay, CNRS, Laboratoire de mathématiques d'Orsay, ParMA, Inria Saclay, 91405, Orsay, France. email: luca.nenna@universite-paris-saclay.fr</footnote>

<footnote>‡Department of Mathematical and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta, Canada, T6G 2G1. email: pass@ualberta.ca</footnote>
# An ordinary differential equation for entropic optimal transport and its linearly constrained variants

Joshua Zoen-Git Hiew,\* Luca Nenna† and Brendan Pass ‡

March 28, 2024

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 \to \mathbb{R}$, the Monge-Kantorovich optimal transport problem is to

<footnote>\*Department of Mathematical and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta, Canada, T6G 2G1. email: joshuazo@ualberta.ca</footnote>

<footnote>†Université Paris-Saclay, CNRS, Laboratoire de mathématiques d'Orsay, ParMA, Inria Saclay, 91405, Orsay, France. email: luca.nenna@universite-paris-saclay.fr</footnote>

<footnote>‡Department of Mathematical and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta, Canada, T6G 2G1. email: pass@ualberta.ca</footnote>
minimize

$$ \int_{X^1 \times X^2} c(x^1, x^2) d\gamma(x^1, x^2) \qquad (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 \left( \frac{d\gamma}{d(\mu^1 \otimes \mu^2)} \right) 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 [7, 16]. 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 [8, 11, 14, 23, 26, 30].

In this paper, we show that, for discrete marginals<fnref n="1" />, 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 [18, 31] 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, 20, 21]. 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

<footnote n="1">¹General marginals will be dealt with in a separate, forthcoming paper.</footnote>
minimize

$$ \int_{X^1 \times X^2} c(x^1, x^2) d\gamma(x^1, x^2) \qquad (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 \left( \frac{d\gamma}{d(\mu^1 \otimes \mu^2)} \right) 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 [7, 16]. 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 [8, 11, 14, 23, 26, 30].

In this paper, we show that, for discrete marginals<fnref n="1" />, 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 [18, 31] 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, 20, 21]. 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

<footnote n="1">¹General marginals will be dealt with in a separate, forthcoming paper.</footnote>
(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 pre-
sented here to compute solutions to a variety of problems. A key computational advan-
tage 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 es-
sentially 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 func-
tion 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 \to \infty$). Taylor expansions of the cost around the fully unregularized,
optimal transport limit ($\varepsilon \to \infty$, or $\eta \to 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''(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.
(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 pre-
sented here to compute solutions to a variety of problems. A key computational advan-
tage 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 es-
sentially 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 func-
tion 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 \to \infty$). Taylor expansions of the cost around the fully unregularized,
optimal transport limit ($\varepsilon \to \infty$, or $\eta \to 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''(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 measures<fnref n="2" /> on respective bounded domains $X^1, ..., X^n$ in $\mathbb{R}^d$. Set $X := 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, \qquad (2.1) $$

where $\Pi^Q(\mu^1, ..., \mu^n) \subset \Pi(\mu^1, ..., \mu^n)$ is the set of couplings such that $\int q d\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 := \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, \qquad (2.2) $$

where $\Psi(Q, c) := \{(\psi, q) | \sum_{i=1}^{n} \psi^i + q \le 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):

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

<footnote n="2">²From 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.</footnote>
## 2. Optimal Transport with Extra Linear constraints

### 2.1. The basic problem

Let $\mu^1, ..., \mu^n$ be probability measures<fnref n="2" /> on respective bounded domains $X^1, ..., X^n$ in $\mathbb{R}^d$. Set $X := 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, \qquad (2.1) $$

where $\Pi^Q(\mu^1, ..., \mu^n) \subset \Pi(\mu^1, ..., \mu^n)$ is the set of couplings such that $\int q d\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 := \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, \qquad (2.2) $$

where $\Psi(Q, c) := \{(\psi, q) | \sum_{i=1}^{n} \psi^i + q \le 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):

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

<footnote n="2">²From 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.</footnote>
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:

$$\sup_{(\psi,q) \in \Psi(Q,c)} \sum_{i=1}^{n} \int_{X^i} \psi^i d\mu^i - \eta \int_{X} \exp \left( \frac{\sum_{i=1}^{n} \psi^i + q - c}{\eta} \right) d(\otimes_{i=1}^{n} \mu^i). \quad (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 \to 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 nontrivial $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 \to \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 \ge 0$ to obtain the problem

$$\inf_{\gamma \in \Pi^Q(\mu^1, \dots, \mu^n)} \int_X \varepsilon c d\gamma + \eta H_{\otimes_{i=1}^n \mu^i}(\gamma). \quad (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, \dots, \mu^n)} H_{\otimes_{i=1}^n \mu^i}(\gamma).$$

Note that in the unconstrained problem $Q = \{0\}$, a trivial solution arises—the product

<footnote n="3">³One could of course more simply set $\varepsilon = 1/\eta$ in the original problem (2.3), which corresponds to taking $\eta = 1$ in (2.5). However, we find it more convenient to retain the parameter $\eta$ here was well. In particular, below we will characterize solutions by an ODE in $\varepsilon$. Keeping $\eta$ will allow us to solve (2.3) by solving the ODE up to $\varepsilon = 1$ for any value of $\eta$, instead of having to consider all values of $\varepsilon$ from 0 up to $\infty$.</footnote>
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:

$$\sup_{(\psi,q) \in \Psi(Q,c)} \sum_{i=1}^{n} \int_{X^i} \psi^i d\mu^i - \eta \int_{X} \exp \left( \frac{\sum_{i=1}^{n} \psi^i + q - c}{\eta} \right) d(\otimes_{i=1}^{n} \mu^i). \quad (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 \to 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 nontrivial $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 \to \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 \ge 0$ to obtain the problem

$$\inf_{\gamma \in \Pi^Q(\mu^1, \dots, \mu^n)} \int_X \varepsilon c d\gamma + \eta H_{\otimes_{i=1}^n \mu^i}(\gamma). \quad (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, \dots, \mu^n)} H_{\otimes_{i=1}^n \mu^i}(\gamma).$$

Note that in the unconstrained problem $Q = \{0\}$, a trivial solution arises—the product

<footnote n="3">³One could of course more simply set $\varepsilon = 1/\eta$ in the original problem (2.3), which corresponds to taking $\eta = 1$ in (2.5). However, we find it more convenient to retain the parameter $\eta$ here was well. In particular, below we will characterize solutions by an ODE in $\varepsilon$. Keeping $\eta$ will allow us to solve (2.3) by solving the ODE up to $\varepsilon = 1$ for any value of $\eta$, instead of having to consider all values of $\varepsilon$ from 0 up to $\infty$.</footnote>
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 minimizer $\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} \exp \left( \frac{\sum_{i=1}^{n} \psi^i + q - \varepsilon c}{\eta} \right) d(\otimes_{i=1}^{n} \mu^i). \quad (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_{x_i}^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 \ge 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 $\mu = (\mu_x)_{x \in X} = (\otimes_{i=1}^n \mu_{x_i}^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}{ll} \min & \sum_{x \in X} \varepsilon c_x \gamma_x + \eta \sum_{x \in X} \log \left( \frac{\gamma_x}{\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, \dots, n \\ & \sum_{x \in X} q_x^j \gamma_x = 0 & j = 1, \dots, K \\ & \gamma_x \ge 0 & \forall x \in X, \end{array} \quad (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) \mu_x. \quad (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 $\mu := \otimes_{i=1}^n \mu^i$, so that $\mu_\ell = \otimes_{i=1}^n \mu_{x_i}^i$ for all $x_l = (x_l^1, x_l^2, ..., x_l^n) \in X$. We then note that $\sum_{i=1}^n \psi_{x_\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}. \quad (2.9) $$

It is easy to see that the above primal problem (2.7) can be expressed as:
measure $\gamma \in \Pi^Q(\mu^1, ..., \mu^n)$ itself.

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

The dual of (2.5) is simply (2.4) with the cost function $\epsilon 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} \exp \left( \frac{\sum_{i=1}^{n} \psi^i + q - \epsilon c}{\eta} \right) d(\otimes_{i=1}^{n} \mu^i). \quad (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_{x_i}^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 \ge 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 $\mu = (\mu_x)_{x \in X} = (\otimes_{i=1}^n \mu_{x_i}^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}{ll} \min & \sum_{x \in X} \epsilon c_x \gamma_x + \eta \sum_{x \in X} \log \left( \frac{\gamma_x}{\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, \dots, n \\ & \sum_{x \in X} q_x^j \gamma_x = 0 & j = 1, \dots, K \\ & \gamma_x \ge 0 & \forall x \in X, \end{array} \quad (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 - \epsilon c_x}{\eta} \right) \mu_x. \quad (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 $\mu := \otimes_{i=1}^n \mu^i$, so that $\mu_\ell = \otimes_{i=1}^n \mu_{x_i}^i$ for all $x_l = (x_l^1, x_l^2, ..., x_l^n) \in X$. We then note that $\sum_{i=1}^n \psi_{x_\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}. \quad (2.9) $$

It is easy to see that the above primal problem (2.7) can be expressed as:
$$
\begin{array}{ll}
\min & \varepsilon c^T \gamma + \eta \sum_{\ell=1}^{m} \log(\frac{\gamma_{\ell}}{\mu_{\ell}}) \gamma_{\ell} \\
\text{subject to} & A^T \gamma = b \\
& \gamma \ge 0
\end{array}
\qquad (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) \mu_{\ell}, \quad (2.11)
$$

where $b = [\mu, 0]^T = [\mu^1 \mu^2 \dots \mu^n, 0, \dots, 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, \dots, \mu^n)$ we have $\mu = \gamma^T B$ and $0 = \gamma^T C$.

**Lemma 2.3.** Assume that $\Pi^Q(\mu^1, \mu^2, \ldots, \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) \mu_{\ell} \quad (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 (2.11) and (2.12) have the same range provided the matrices

$$
P = \begin{bmatrix} b^T \\ A \end{bmatrix} = \begin{bmatrix} \mu & 0 \\ B & C \end{bmatrix} \quad \text{and} \quad \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
$[\hat{\mu}]$ has the same range as $[\hat{\mu}]$. 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, \dots, x_{\ell}^n)$, so if $B\psi = 0$, $\sum_{i=1}^{n} \psi_{x_{\ell}^i} = 0$ for
each $l$.
$$
\begin{array}{ll}
\min & \varepsilon c^T \gamma + \eta \sum_{\ell=1}^{m} \log(\frac{\gamma_{\ell}}{\mu_{\ell}})\gamma_{\ell} \\
\text{subject to} & A^T \gamma = b \\
& \gamma \ge 0
\end{array}
\qquad (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) \mu_{\ell}, \qquad (2.11)
$$

where $b = [\mu, 0]^T = [\mu^1 \mu^2 \dots \mu^n, 0, \dots, 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, \dots, \mu^n)$ we have $\mu = \gamma^T B$ and $0 = \gamma^T C$.

**Lemma 2.3.** Assume that $\Pi^Q(\mu^1, \mu^2, \dots, \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) \mu_{\ell} \quad (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 (2.11) and (2.12) have the same range provided the matrices

$$
P = \begin{bmatrix} b^T \\ A \end{bmatrix} = \begin{bmatrix} \mu & 0 \\ B & C \end{bmatrix} \quad \text{and} \quad \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
$[\hat{\mu}]$ has the same range as $[\hat{\mu}]$. 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, \dots, x_{\ell}^n)$, so if $B\psi = 0$, $\sum_{i=1}^{n} \psi_{x_{\ell}^i} = 0$ for
each $l$.
Multiplying by $\mu_\ell$ and summing over $l$, we have

$$0 = \sum_{l=1}^{m} \mu_{\ell} \sum_{i=1}^{n} \psi_{x_{\ell}^{i}} = \sum_{i=1}^{n} \sum_{l:x_{\ell}^{i}=x^{i}} \mu_{\ell} \psi_{x^{i}} = \sum_{i=1}^{n} \mu_{x^{i}}^{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 $[\hat{\mu} \ 0]$ and $[\hat{\mu} \ 0]$ 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 belong to the column space of $\hat{B}$, we set $\hat{C}^1 = q^1$. Otherwise, we claim that $[\hat{\mu}]$ and $[\hat{\mu} \ 0]$ have the same range. If $\hat{B}\varphi = q^1$ and $\gamma \in \Pi^Q(\mu^1, \dots, \mu^n)$, then $\hat{B}^T\gamma = \hat{\mu}^T$ and $q^{1T}\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 $[u \ v]$ of the appropriate dimension, we have

$$[\hat{\mu} \ 0] [u \ v] = [\hat{\mu}u \ \hat{B}u + q^1 v] = [\hat{\mu}(u + \varphi v)] = [\hat{B}] [u + \varphi v], \quad (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. \quad (2.14)$$

Then with the similar argument as above, $[\hat{\mu} \ 0]$ and $[\hat{\mu} \ 0 \ 0]$ 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 continue 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} \quad \text{and} \quad \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. $\square$
Multiplying by $\mu_\ell$ and summing over $l$, we have

$$0 = \sum_{l=1}^{m} \mu_{\ell} \sum_{i=1}^{n} \psi_{x_{\ell}^{i}} = \sum_{i=1}^{n} \sum_{l:x_{\ell}^{i}=x^{i}} \mu_{\ell} \psi_{x^{i}} = \sum_{i=1}^{n} \mu_{x^{i}}^{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 $[\hat{\mu} \ 0]$ and $[\hat{\mu} \ 0]$ 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 belong to the column space of $\hat{B}$, we set $\hat{C}^1 = q^1$. Otherwise, we claim that $[\hat{\mu}]$ and $[\hat{\mu} \ 0]$ 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^{1T}\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 $[u \ v]$ of the appropriate dimension, we have

$$[\hat{\mu} \ 0] [u \ v] = [\hat{\mu}u \ \hat{B}u + q^1v] = [\hat{\mu}(u + \varphi v)] = [\hat{B}] [u + \varphi v], \quad (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. \quad (2.14)$$

Then with the similar argument as above, $[\hat{\mu} \ 0]$ and $[\hat{\mu} \ 0 \ 0]$ 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 continue 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} \quad \text{and} \quad \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. $\square$
# 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}. \quad (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). \quad (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. \quad (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=1}^{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:

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

$$ D_{\varphi\varphi}^2 \Phi(\varphi, \varepsilon) = \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 D A, \quad (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 D A$ is positive definite, and as the entries are continuous, uniformly bounded below $A^T D A \ge d > 0$ on the compact domain $V \times T$. $\square$
# 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 $η > 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) \mu_{\ell}. \quad (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). \quad (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. \quad (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=1}^{m} \exp(\frac{1}{\eta}(A_{\ell} \varphi - \varepsilon c_{\ell}))\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:

$$ \nabla_{\varphi} \Phi(\varphi, \varepsilon) = -b^T + \sum_{\ell=1}^{m} A_{\ell}^T \exp(\frac{1}{\eta}(A_{\ell} \varphi - \varepsilon c_{\ell})) \mu_{\ell} $$

$$ D_{\varphi\varphi}^2 \Phi(\varphi, \varepsilon) = \frac{1}{\eta} \sum_{\ell=1}^{m} A_{\ell}^T \exp(\frac{1}{\eta}(A_{\ell}\varphi - \varepsilon c_{\ell})) A_{\ell} \mu_{\ell}. $$

The Hessian can be rewritten as:

$$ D_{\varphi\varphi}^2 \Phi(\varphi, \varepsilon) = \frac{1}{\eta} A^T D A, \quad (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}))\mu_{\ell}$. This implies that $D$ and therefore $A^T D A$ is positive definite, and as the entries are continuous, uniformly bounded below $A^T D A \ge d > 0$ on the compact domain $V \times T$. $\square$
Differentiating (3.3) with respect to $\epsilon$ we obtain:

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

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

$$ \left\{ \begin{array}{l} \frac{d\varphi}{d\epsilon}(\epsilon) = -[D_{\varphi\varphi}^2 \Phi(\varphi(\epsilon), \epsilon)]^{-1} \frac{\partial}{\partial \epsilon} \nabla_{\varphi} \Phi(\varphi(\epsilon), \epsilon) \\ \varphi(0) = \phi_0. \end{array} \right. \qquad (3.5) $$

Here, $\phi_0$ represents the initial condition of the ODE corresponding to the optimal $\varphi$ solving (3.2) when $\epsilon = 0$. This corresponds to finding the measure $\gamma$ which minimizes the relative entropy among the admissible class $\Pi^Q(\mu^1, \mu^2, \dots, \mu^n)$, and is often easier to find than the solution to (3.2) for $\epsilon > 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 \epsilon} \nabla_{\varphi} \Phi(\varphi, \epsilon) = -\frac{1}{\eta} \sum_{\ell=1}^{m} c_{\ell} A_{\ell}^{T} \exp\left(\frac{1}{\eta}(A_{\ell}\varphi - \epsilon c_{\ell})\right) \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(\epsilon)$ is clearly continuous, there exist an $M > 0$ such that $||\varphi(\epsilon)|| \le M, \epsilon \in [0, 1]$. Therefore, we aim to prove the well-posedness of (3.2) over $V \times [0, 1]$ where:

$$ V := \{\varphi \in \mathbb{R}^e \mid ||\varphi|| \le M\}. $$

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

$$ \left\{ \begin{array}{l} \frac{d\varphi}{d\epsilon}(\epsilon) = -[D_{\varphi\varphi}^2 \Phi(\varphi(\epsilon), \epsilon)]^{-1} \frac{\partial}{\partial \epsilon} \nabla_{\varphi} \Phi(\varphi(\epsilon), \epsilon) \\ \varphi(0) = \phi_0 \end{array} \right. . \qquad (3.6) $$

Here, $\phi_0$ is the solution to the problem:

$$ \min_{\varphi} \Phi(\varphi, 0), \qquad (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(\epsilon)$ and the fact that it satisfies the ODE is demonstrated
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$:

$$ \left\{ \begin{array}{l} \frac{d\varphi}{d\varepsilon}(\varepsilon) = -[D_{\varphi\varphi}^2 \Phi(\varphi(\varepsilon), \varepsilon)]^{-1} \frac{\partial}{\partial \varepsilon} \nabla_{\varphi} \Phi(\varphi(\varepsilon), \varepsilon) \\ \varphi(0) = \phi_0. \end{array} \right. \qquad (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 admissible class $\Pi^Q(\mu^1, \mu^2, \dots, \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\left(\frac{1}{\eta}(A_{\ell}\varphi - \varepsilon c_{\ell})\right) \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)|| \le M, \varepsilon \in [0, 1]$. Therefore, we aim to prove the well-posedness of (3.2) over $V \times [0, 1]$ where:

$$ V := \{\varphi \in \mathbb{R}^e \mid ||\varphi|| \le M\}. $$

**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:

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

Here, $\phi_0$ is the solution to the problem:

$$ \min_{\varphi} \Phi(\varphi, 0), \qquad (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 - [D_{\varphi\varphi}^2 \Phi(\varphi(\varepsilon), \varepsilon)]^{-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}$). $\square$

**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 \ge 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 \ge 0$. Together with the differentiability of $\Phi$ in $\varphi$, by Theorem 14.6.7 of [29], the result follows. $\square$

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) := \inf_{\gamma \in \Pi(\mu, \nu)} \int_{X \times Y} \varepsilon c(x, y) d\gamma(x, y) + \eta H_{\mu \otimes \nu}(\gamma), \quad (4.1)
$$
in section 3.1.

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

$$
(\varphi, \varepsilon) \mapsto - [D_{\varphi\varphi}^2 \Phi(\varphi(\varepsilon), \varepsilon)]^{-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}$). $\square$

**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 \ge 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 \ge 0$. Together with the differentiability of $\Phi$ in $\varphi$, by Theorem 14.6.7 of [29], the result follows. $\square$

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) := \inf_{\gamma \in \Pi(\mu, \nu)} \int_{X \times Y} \varepsilon c(x, y) d\gamma(x, y) + \eta H_{\mu \otimes \nu}(\gamma), \quad (4.1)
$$
while the dual is:

$$D(\varepsilon) := \sup_{u,v} \int_X u d\mu + \int_Y v d\nu - \eta \int_{X \otimes Y} e^{\frac{u+v-\varepsilon c}{\eta}} d\mu d\nu. \quad (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) := \inf_{u,v} \Phi(u, v, \varepsilon) = -D(\varepsilon) = -P(\varepsilon), \quad (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}. \quad (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 = -v_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]. \quad (4.5)$$

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

$$\Phi(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, \quad (4.6)$$

where we have abused the notation for $\Phi$; note that we can now write $C(\varepsilon) = \inf_u \Phi(u, \varepsilon)$.
while the dual is:

$$D(\varepsilon) := \sup_{u,v} \int_X u d\mu + \int_Y v d\nu - \eta \int_{X \otimes Y} e^{\frac{u+v-\varepsilon c}{\eta}} d\mu d\nu. \quad (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) := \inf_{u,v} \Phi(u, v, \varepsilon) = -D(\varepsilon) = -P(\varepsilon), \quad (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}. \quad (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 = -v_j + \sum_r \exp \left[ \frac{1}{\eta} (u_r + v_j - \varepsilon c_{rj}) \right] \mu_r v_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]. \quad (4.5)$$

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

$$\Phi(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, \quad (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,

$$ 
\begin{aligned}
-P'(\varepsilon) &= -D'(\varepsilon) = C'(\varepsilon) = \frac{d}{d\varepsilon}\Phi(u(\varepsilon), \varepsilon) \\
&= \nabla_u\Phi(u(\varepsilon), \varepsilon)\frac{du}{d\varepsilon} + \frac{\partial}{\partial\varepsilon}\Phi(u(\varepsilon), \varepsilon) \\
&= \frac{\partial}{\partial\varepsilon}\Phi(u(\varepsilon), \varepsilon) \\
&= -\sum_s \left[ \frac{\sum_r e_{rs}c_{rs}}{\sum_r e_{rs}} \right] \nu_s,
\end{aligned} 
\quad (4.7) $$

where we have used the first order optimality condition $\nabla_u\Phi(u(\varepsilon), \varepsilon) = 0, \forall \varepsilon \ge 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'(0) = -\sum_s \left[ \frac{\sum_r c_{rs}\mu_r}{\sum_r \mu_r} \right] \nu_s = -\mathbb{E}[c(X, Y)]. \quad (4.8) $$

The expectation in (4.8) is taken with respect to the product measure $\mu \otimes \nu$.<fnref n="4" /> For the second order derivative,

$$ 
\begin{aligned}
C''(\varepsilon) &= \frac{d}{d\varepsilon} \left( \frac{\partial}{\partial\varepsilon} \Phi(u(\varepsilon), \varepsilon) \right) \\
&= \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) \\
&= \left[ \nabla_u \frac{\partial}{\partial\varepsilon} \Phi(u(\varepsilon), \varepsilon) \right]^T \left[ -D_{uu}^2 \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). 
\end{aligned} 
\quad (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''(0) = \frac{1}{\eta} \left( (\mathbb{E}[c(X, Y)])^2 + \mathbb{E}[c^2(X, Y)] - \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2] - \mathbb{E}[(\mathbb{E}[c(X, Y)|Y])^2] \right), \quad (4.10) $$

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

<footnote n="4">⁴Note 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.</footnote>
For simplicity, we denote $e_{ij} = \exp[\frac{1}{\eta}(u_i - \varepsilon c_{ij})]\mu_i$. Note that,

$$ 
\begin{aligned}
-P'(\varepsilon) &= -D'(\varepsilon) = C'(\varepsilon) = \frac{d}{d\varepsilon}\Phi(u(\varepsilon), \varepsilon) \\
&= \nabla_u\Phi(u(\varepsilon), \varepsilon)\frac{du}{d\varepsilon} + \frac{\partial}{\partial\varepsilon}\Phi(u(\varepsilon), \varepsilon) \\
&= \frac{\partial}{\partial\varepsilon}\Phi(u(\varepsilon), \varepsilon) \\
&= -\sum_s \left[ \frac{\sum_r e_{rs}c_{rs}}{\sum_r e_{rs}} \right] \nu_s,
\end{aligned} 
\quad (4.7) $$

where we have used the first order optimality condition $\nabla_u\Phi(u(\varepsilon), \varepsilon) = 0, \forall \varepsilon \ge 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'(0) = -\sum_s \left[ \frac{\sum_r c_{rs}\mu_r}{\sum_r \mu_r} \right] \nu_s = -\mathbb{E}[c(X, Y)]. \quad (4.8) $$

The expectation in (4.8) is taken with respect to the product measure $\mu \otimes \nu$.<fnref n="4" /> For the second order derivative,

$$ 
\begin{aligned}
C''(\varepsilon) &= \frac{d}{d\varepsilon} \left( \frac{\partial}{\partial\varepsilon} \Phi(u(\varepsilon), \varepsilon) \right) \\
&= \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) \\
&= \left[ \nabla_u \frac{\partial}{\partial\varepsilon} \Phi(u(\varepsilon), \varepsilon) \right]^T \left[ -D_{uu}^2 \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).
\end{aligned} 
\quad (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''(0) = \frac{1}{\eta} \left( (\mathbb{E}[c(X, Y)])^2 + \mathbb{E}[c^2(X, Y)] - \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2] - \mathbb{E}[(\mathbb{E}[c(X, Y)|Y])^2] \right), \quad (4.10) $$

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

<footnote n="4">⁴Note 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.</footnote>
$$C''(0) = \frac{1}{\eta} \left( \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) - \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) \right).$$

We note that higher order derivatives can be evaluated similarly; for instance, one can find $C'''(\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'''(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''(0) = -P''(0) \ge 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

$$\begin{aligned} \mathrm{Var}[\mathbb{E}[c(X, Y)|X]] &= \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2] - (E[c])^2 \\ \mathbb{E}[\mathrm{Var}[c(X, Y)|X]] &= \mathbb{E}[c^2] - \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2], \end{aligned}$$

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}[\mathrm{Var}[c(X, Y)|X]] \ge \mathrm{Var}[\mathbb{E}[c(X, Y)|X]].$$

However, it is easy to calculate that $\mathbb{E}[\mathrm{Var}[c(X, Y)|X]] = 0$ but $\mathrm{Var}[\mathbb{E}[c(X, Y)|X]] = 4\mathrm{Var}[X] \ne 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 method$^5$. Furthermore, we have chosen examples for which explicit solutions $\gamma_0$ of the unregularized problem (2.1) are known. This then yields a

<footnote>$^5$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</footnote>
$$C''(0) = \frac{1}{\eta} \left( \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) \right. \\ \left. - \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) \right).$$

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

**Remark 4.1.** $P(\epsilon)$ as defined by (4.1) is clearly concave, as an infimum of affine functions. Therefore, we must have $C''(0) = -P''(0) \ge 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

$$\begin{aligned} \mathrm{Var}[\mathbb{E}[c(X, Y)|X]] &= \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2] - (E[c])^2 \\ \mathbb{E}[\mathrm{Var}[c(X, Y)|X]] &= \mathbb{E}[c^2] - \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2], \end{aligned}$$

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}[\mathrm{Var}[c(X, Y)|X]] \ge \mathrm{Var}[\mathbb{E}[c(X, Y)|X]].$$

However, it is easy to calculate that $\mathbb{E}[\mathrm{Var}[c(X, Y)|X]] = 0$ but $\mathrm{Var}[\mathbb{E}[c(X, Y)|X]] = 4\mathrm{Var}[X] \ne 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 $\epsilon$ 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 method<fnref n="5" />. Furthermore, we have chosen examples for which explicit solutions $\gamma_0$ of the unregularized problem (2.1) are known. This then yields a

<footnote n="5">⁵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(\epsilon)$ is in fact the main object of</footnote>
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 \le \inf_{\gamma \in \Pi^Q(\mu^1, \dots, \mu^n)} \int_X cd\gamma + \eta H_{\otimes_{i=1}^n \mu^i}(\gamma) \le \int_X cd\gamma_0 + \eta H_{\otimes_{i=1}^n \mu^i}(\gamma_0). \quad (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), \quad (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. \quad (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:

$$ v_j = -\eta \log \left[ \sum_r \exp \left[ \frac{1}{\eta} (u_r - \varepsilon c_{rj}) \right] \mu_r \right] =: v_j(u). \quad (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. \quad (5.5) $$

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

<footnote>interest there, rather than the final value $\gamma(1)$, and this curve of marginals is not as readily obtained by the Sinkhorn method.</footnote>
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 \le \inf_{\gamma \in \Pi^Q(\mu^1, \dots, \mu^n)} \int_X cd\gamma + \eta H_{\otimes_{i=1}^n \mu^i}(\gamma) \le \int_X cd\gamma_0 + \eta H_{\otimes_{i=1}^n \mu^i}(\gamma_0). \quad (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), \quad (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. \quad (5.3) $$

It is 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:

$$ v_j = -\eta \log \left[ \sum_r \exp \left[ \frac{1}{\eta} (u_r - \varepsilon c_{rj}) \right] \mu_r \right] =: v_j(u). \quad (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. \quad (5.5) $$

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

<footnote>interest there, rather than the final value $\gamma(1)$, and this curve of marginals is not as readily obtained by the Sinkhorn method.</footnote>
**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:

$$ 
\begin{aligned}
\nabla \bar{\Phi}(u, \varepsilon) &= \nabla_u \Phi(u, v(u), \varepsilon) + \nabla_v \Phi(u, v(u), \varepsilon) Dv(u) \\
D^2 \bar{\Phi}(u, \varepsilon) &= 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) \\
&\quad + Dv(u)^T D_{vv}^2 \Phi(u, v(u), \varepsilon) Dv(u) + \nabla_v \Phi(u, v(u), \varepsilon) D^2 v(u).
\end{aligned}
$$

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}. \qquad (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, \quad 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:

$$ \left( D_{uu}^2 \Phi(u(\varepsilon), \varepsilon) \right)_{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, \qquad (5.7) $$

where $\delta_{i_1 i_2}$ is the Kronecker delta and the mixed derivative is:

$$ \left( \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right)_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. \qquad (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).
**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:

$$ 
\begin{aligned}
\nabla \bar{\Phi}(u, \varepsilon) &= \nabla_u \Phi(u, v(u), \varepsilon) + \nabla_v \Phi(u, v(u), \varepsilon) Dv(u) \\
D^2 \bar{\Phi}(u, \varepsilon) &= 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) \\
&\quad + Dv(u)^T D_{vv}^2 \Phi(u, v(u), \varepsilon) Dv(u) + \nabla_v \Phi(u, v(u), \varepsilon) D^2 v(u).
\end{aligned}
$$

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}. \qquad (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:

$$ \left( D_{uu}^2 \Phi(u(\varepsilon), \varepsilon) \right)_{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, \qquad (5.7) $$

where $\delta_{i_1 i_2}$ is the Kronecker delta and the mixed derivative is:

$$ \left( \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right)_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. \qquad (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><img src="image_2.png" /><figcaption>Figure 1: Evolution of the primal coupling at 16 evenly distributed values of $\epsilon$ between 0 and 1 calculated by the ODE method for two marginal optimal transport with attractive cost $c(x, y) = |y - x|^2$.</figcaption></figure>

<table><caption>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$.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.0050</td><td>0.0052</td></tr><tr><td>Optimal unregularized primal value</td><td>0</td><td>0</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.0092</td><td>0.0092</td></tr><tr><td>Iterations</td><td>100</td><td>324</td></tr><tr><td>CPU time (sec)</td><td>1.40</td><td>0.04</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the primal coupling at 16 evenly distributed values of $\epsilon$ between 0 and 1 calculated by the ODE method for two marginal optimal transport with attractive cost $c(x, y) = |y - x|^2$.</figcaption></figure>

<table><caption>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$.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.0050</td><td>0.0052</td></tr><tr><td>Optimal unregularized primal value</td><td>0</td><td>0</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.0092</td><td>0.0092</td></tr><tr><td>Iterations</td><td>100</td><td>324</td></tr><tr><td>CPU time (sec)</td><td>1.40</td><td>0.04</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of the primal coupling at 16 evenly distributed values of $\epsilon$ 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|)$.</figcaption></figure>

## 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(\epsilon)$ 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
<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of the primal coupling at 16 evenly distributed values of $\epsilon$ 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|)$.</figcaption></figure>

## 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(\epsilon)$ 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
<table><caption>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|)$.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.5033</td><td>0.5080</td></tr><tr><td>Optimal unregularized primal value</td><td>0.5024</td><td>0.5024</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.5117</td><td>0.5117</td></tr><tr><td>Iterations</td><td>100</td><td>684</td></tr><tr><td>CPU time (sec)</td><td>1.49</td><td>0.15</td></tr></tbody></table>

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}. \quad (5.9) $$

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

$$ w_k = -\eta \log \left[ \sum_{rs} \exp \left[ \frac{1}{\eta} (u_r + v_s - \varepsilon c_{rsk}) \right] \mu_r \theta_s \right], \quad (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. \quad (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.
<table><caption>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|)$.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.5033</td><td>0.5080</td></tr><tr><td>Optimal unregularized primal value</td><td>0.5024</td><td>0.5024</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.5117</td><td>0.5117</td></tr><tr><td>Iterations</td><td>100</td><td>684</td></tr><tr><td>CPU time (sec)</td><td>1.49</td><td>0.15</td></tr></tbody></table>

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}. \quad (5.9) $$

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

$$ w_k = -\eta \log \left[ \sum_{rs} \exp \left[ \frac{1}{\eta} (u_r + v_s - \varepsilon c_{rsk}) \right] \mu_r \theta_s \right], \quad (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. \quad (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><img src="image_4.png" /><figcaption>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</figcaption></figure>

<table><caption>Table 3: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for three marginal optimal transport problem.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>1.9163</td><td>1.9193</td></tr><tr><td>Optimal unregularized primal value</td><td>1.9137</td><td>1.9137</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>1.9647</td><td>1.9647</td></tr><tr><td>Iterations</td><td>100</td><td>1006</td></tr><tr><td>CPU time (sec)</td><td>12.40</td><td>47.46</td></tr></tbody></table>

### 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)\}. \qquad (5.12)$$

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

$$\int_Y y d\pi_x(y) = x, \qquad (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. \qquad (5.14)$$
<figure><img src="image_4.png" /><figcaption>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</figcaption></figure>

<table><caption>Table 3: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for three marginal optimal transport problem.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>1.9163</td><td>1.9193</td></tr><tr><td>Optimal unregularized primal value</td><td>1.9137</td><td>1.9137</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>1.9647</td><td>1.9647</td></tr><tr><td>Iterations</td><td>100</td><td>1006</td></tr><tr><td>CPU time (sec)</td><td>12.40</td><td>47.46</td></tr></tbody></table>

### 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)\}. \qquad (5.12)$$

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

$$\int_Y y d\pi_x(y) = x, \qquad (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. \qquad (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), \qquad (5.15) $$

and the corresponding dual problem:

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

subject to $u(x) + v(y) + g(x)(y - x) \le 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), \qquad (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}. \qquad (5.18) $$

The first order optimality condition implies

$$ u_i = -\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] \qquad (5.19) $$

$$ v_j = -\eta \log \left[ \sum_r \exp \left[ \frac{1}{\eta} (u_r + g_r(y_r - x_r) - \varepsilon c_{rj}) \right] \mu_r \right] \qquad (5.20) $$

$$ 0 = \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. \qquad (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) := - \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. \qquad (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
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), \qquad (5.15) $$

and the corresponding dual problem:

$$ \sup \int_X u(x) d\mu(x) + \int_Y v(y) d\nu(y), \qquad (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), \qquad (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. \qquad (5.18) $$

The first order optimality condition implies

$$ u_i = -\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] \qquad (5.19) $$

$$ v_j = -\eta \log \left[ \sum_r \exp \left[ \frac{1}{\eta} (u_r + g_r(y_r - x_r) - \varepsilon c_{rj}) \right] \mu_r \right] \qquad (5.20) $$

$$ 0 = \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. \qquad (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) := - \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. \qquad (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 $\epsilon$. 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 $\epsilon = 0$, by the Sinkhorn algorithm. Note that by our numerical experiment, the running time for the Sinkhorn algorithm is much faster for $\epsilon = 0$ then $\epsilon = 1$. This keeps the overall calculation time for the ODE method competitive with the pure Sinkhorn method to minimize the $\epsilon = 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 $\epsilon = 1$ for this numerical example. We employ 25 steps to discretize $\epsilon$ 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><img src="image_5.png" /><figcaption>Figure 4: The optimal measure for martingale optimal transport with cost $c(x, y) = \exp(-x) \cdot y^2$ generated by ODE</figcaption></figure>

## 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
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><img src="image_5.png" /><figcaption>Figure 4: The optimal measure for martingale optimal transport with cost $c(x, y) = \exp(-x) \cdot y^2$ generated by ODE</figcaption></figure>

## 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
<table><caption>Table 4: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for martingale optimal transport problem.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.2990</td><td>0.2990</td></tr><tr><td>Optimal unregularized primal value</td><td>0.2964</td><td>0.2964</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.3211</td><td>0.3211</td></tr><tr><td>Iterations</td><td>25</td><td>152</td></tr><tr><td>CPU time (sec)</td><td>1.44</td><td>3.38</td></tr></tbody></table>

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, \quad \mathbb{E}_{\pi}[Z|X, Y] = Y \quad \Longleftrightarrow \quad \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

$$\begin{aligned} \Phi(u, v, w, g, h, \varepsilon) = & - \sum_r u_r \mu_r - \sum_s v_s \theta_s - \sum_t w_t \nu_t \quad (5.23) \\ & + \eta \sum_{r,s,t} \exp \left[ \frac{1}{\eta} \left( u_r + v_s + w_t + g_r(y_s - x_r) \right. \right. \\ & \qquad \left. \left. + h_{rs}(z_t - y_s) - \varepsilon c_{rst} \right) \right] \mu_r \theta_s \nu_t. \end{aligned}$$

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

$$u_i = -\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] \quad (5.24)$$

$$v_j = -\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]$$

$$w_k = -\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]$$
<table><caption>Table 4: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for martingale optimal transport problem.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.2990</td><td>0.2990</td></tr><tr><td>Optimal unregularized primal value</td><td>0.2964</td><td>0.2964</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.3211</td><td>0.3211</td></tr><tr><td>Iterations</td><td>25</td><td>152</td></tr><tr><td>CPU time (sec)</td><td>1.44</td><td>3.38</td></tr></tbody></table>

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, \quad \mathbb{E}_{\pi}[Z|X, Y] = Y \quad \Longleftrightarrow \quad \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

$$\begin{aligned} \Phi(u, v, w, g, h, \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} \left( u_r + v_s + w_t + g_r(y_s - x_r) \right. \right. \\ & \qquad \left. \left. + h_{rs}(z_t - y_s) - \varepsilon c_{rst} \right) \mu_r \theta_s \nu_t \right]. \end{aligned} \quad (5.23)$$

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

$$u_i = -\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] \quad (5.24)$$

$$v_j = -\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]$$

$$w_k = -\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]$$
$$0 = \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 \\ 0 = \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:

$$\Phi(v, w, g, h, \varepsilon) = - \sum_s v_s \theta_s - \sum_t w_t \nu_t + \eta \sum_r \log \left[ \sum_{s,t} \exp \left[ \frac{1}{\eta} (v_s + w_t + g_r(y_s - x_r) + h_{rs}(z_t - y_s) - \varepsilon c_{rst}) \right] \theta_s \nu_t \right] \mu_r + \eta, \quad (5.25)$$

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 with the 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.
$$0 = \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 \\ 0 = \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:

$$\Phi(v, w, g, h, \varepsilon) = - \sum_s v_s \theta_s - \sum_t w_t \nu_t + \eta \sum_r \log \left[ \sum_{s,t} \exp \left[ \frac{1}{\eta} (v_s + w_t + g_r(y_s - x_r) + h_{rs}(z_t - y_s) - \varepsilon c_{rst}) \right] \theta_s \nu_t \right] \mu_r + \eta, \quad (5.25)$$

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 with the 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><img src="image_6.png" /><figcaption>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</figcaption></figure>

<table><caption>Table 5: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for multi-period martingale transport problem.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.3807</td><td>0.3807</td></tr><tr><td>Optimal unregularized primal value</td><td>0.3767</td><td>0.3767</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.4127</td><td>0.4127</td></tr><tr><td>Iterations</td><td>25</td><td>60</td></tr><tr><td>CPU time (sec)</td><td>8.71</td><td>7.40</td></tr></tbody></table>

## 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
<figure><img src="image_6.png" /><figcaption>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</figcaption></figure>

<table><caption>Table 5: Comparison of the performance between 4-th order Runge-Kutta ODE method and Sinkhorn algorithm for multi-period martingale transport problem.</caption><thead><tr><th></th><th>ODE</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>Computed regularized primal value</td><td>0.3807</td><td>0.3807</td></tr><tr><td>Optimal unregularized primal value</td><td>0.3767</td><td>0.3767</td></tr><tr><td>Optimal unregularized primal value + Entropy</td><td>0.4127</td><td>0.4127</td></tr><tr><td>Iterations</td><td>25</td><td>60</td></tr><tr><td>CPU time (sec)</td><td>8.71</td><td>7.40</td></tr></tbody></table>

### 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 \mu^2}(\gamma).$$

⁶ 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. In a very similar way one can extend this approach

<figure><img src="image_7.png" /><figcaption>Figure 6: Entropic interpolant at different time steps $\varepsilon$</figcaption></figure>

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, \dots, 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) \ge 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,

<footnote n="6">⁶Note 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}$.</footnote>
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 U \otimes \mu^2}(\gamma).$$

⁶ 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. In a very similar way one can extend this approach

<figure><img src="image_7.png" /><figcaption>Figure 6: Entropic interpolant at different time steps $\varepsilon$</figcaption></figure>

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, \dots, 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) \ge 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,

<footnote n="6">⁶Note 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 $U$.</footnote>
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], \quad (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.

## 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 q d\gamma_\infty = \lim_{k \to \infty} \int_X q d\gamma_k = 0, \quad \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 c d\gamma + \eta_k H_{\otimes_{i=1}^n \mu^i}(\gamma). \quad (A.1) $$
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], \quad (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.

## 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 q d\gamma_\infty = \lim_{k \to \infty} \int_X q d\gamma_k = 0, \quad \forall q \in Q$$

since $Q \subset C_b(X)$. Hence $\gamma_\infty \in \Pi^Q$.

<figure><img src="image_8.png" /></figure>

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 c d\gamma + \eta_k H_{\otimes_{i=1}^n \mu^i}(\gamma). \quad (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 \operatorname*{argmin}_{\gamma \text{ solves (2.1)}} 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}) \le \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 \le \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}) \le \int_X cd\bar{\gamma} + \eta_k H_{\otimes_{i=1}^n \mu^i}(\bar{\gamma}). \quad (A.2) $$

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

$$ \int_X cd\gamma_{\eta_k} \ge \int_X cd\bar{\gamma}. \quad (A.3) $$

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

$$ H_{\otimes_{i=1}^n \mu^i}(\gamma_{\eta_k}) \le 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) \le H_{\otimes_{i=1}^n \mu^i}(\bar{\gamma}). $$

This means $\gamma_0$ has the minimum relative entropy among all optimal measures. $\square$
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 \operatorname*{argmin}_{\gamma \text{ solves (2.1)}} 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}) \le \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 \le \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}) \le \int_X cd\bar{\gamma} + \eta_k H_{\otimes_{i=1}^n \mu^i}(\bar{\gamma}). \quad (A.2) $$

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

$$ \int_X cd\gamma_{\eta_k} \ge \int_X cd\bar{\gamma}. \quad (A.3) $$

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

$$ H_{\otimes_{i=1}^n \mu^i}(\gamma_{\eta_k}) \le 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) \le H_{\otimes_{i=1}^n \mu^i}(\bar{\gamma}). $$

This means $\gamma_0$ has the minimum relative entropy among all optimal measures. $\square$
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, \qquad (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 \qquad (A.5) $$

subject to $A_\ell v \le c_\ell$ and

$$ V_0 := \{v \in \mathbb{R}^d | A_\ell v \le c_\ell, -b^T v = D_0\}, $$

then there exist an $M > 0$ such that $||v|| < M, \forall v \in V_0$. Moreover we assume that $\forall \ell, c_\ell \ge 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 \le -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 \le -b^T \bar{v} + \eta_k \le \bar{M} \qquad (A.6) $$

for some $\bar{M} > 0$. The second inequality comes from the fact that $A_\ell \bar{v} \le 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 \le \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} \ne 0$, then

$$ A_\ell \frac{v_k}{||v_k||} \le \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} \le 0$. Moreover, the positive of exponential function and $\mu_i$ also implies that $-b^T \hat{u} \le 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
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}, \qquad (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 \qquad (A.5) $$

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

$$ V_0 := \{v \in \mathbb{R}^d | A_{\ell}v \le c_{\ell}, -b^T v = D_0\}, $$

then there exist an $M > 0$ such that $||v|| < M, \forall v \in V_0$. Moreover we assume that $\forall \ell, c_{\ell} \ge 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} \le -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} \le -b^T \bar{v} + \eta_k \le \bar{M} \qquad (A.6) $$

for some $\bar{M} > 0$. The second inequality comes from the fact that $A_{\ell}\bar{v} \le 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 \le \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} \ne 0$, then

$$ A_{\ell} \frac{v_k}{||v_k||} \le \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} \le 0$. Moreover, the positive of exponential function and $\mu_i$ also implies that $-b^T \hat{u} \le 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

$$ \eta_k \sum_{\ell} \exp \left( \frac{A_{\ell} v_k - \varepsilon c_{\ell}}{\eta_k} \right) \boldsymbol{\mu}_{\ell} \leq \hat{M} \\ A_{\ell} v_k \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 \le c_i$. Again, by (A.6):

$$ -b^T v_0 \le \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} \le -b^T \bar{v}. $$

Therefore, we can conclude that $v_0$ is also an optimal solution to the unregularized problem. $\square$

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

Recall:

$$ C''(\varepsilon) = \left[ \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right]^T \left[ -D_{uu}^2 \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). \quad (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$:

$$ \frac{\partial^2}{\partial \varepsilon^2} \Phi(u(\varepsilon), \varepsilon) = \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 \\ \left( \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right)_i = -\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 \\ \left( D_{uu}^2 \Phi(u(\varepsilon), \varepsilon) \right)_{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, $$

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}[(\mathbb{E}[c(X, Y)|Y])^2] \right). \quad (B.2) $$

Let's deal with the first part of (B.1). With the assumption that $u_0 \equiv 0$, then when
$\hat{M} > 0$ such that

$$ \eta_k \sum_{\ell} \exp \left( \frac{A_{\ell} v_k - \varepsilon c_{\ell}}{\eta_k} \right) \boldsymbol{\mu}_{\ell} \leq \hat{M} \\ A_{\ell} v_k \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 \le c_i$. Again, by (A.6):

$$ -b^T v_0 \le \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} \le -b^T \bar{v}. $$

Therefore, we can conclude that $v_0$ is also an optimal solution to the unregularized problem. $\square$

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

Recall:

$$ C''(\varepsilon) = \left[ \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right]^T \left[ -D_{uu}^2 \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). \quad (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$:

$$ \frac{\partial^2}{\partial \varepsilon^2} \Phi(u(\varepsilon), \varepsilon) = \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 \\ \left( \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right)_i = -\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 \\ \left( D_{uu}^2 \Phi(u(\varepsilon), \varepsilon) \right)_{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, $$

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}[(\mathbb{E}[c(X, Y)|Y])^2] \right). \quad (B.2) $$

Let's deal with the first part of (B.1). With the assumption that $u_0 \equiv 0$, then when
$\epsilon = 0$, for every $i \neq 0$

$$ \left( \nabla_u \frac{\partial}{\partial \epsilon} \Phi(u(0), 0) \right)_i = \frac{\mu_i}{\eta} \left( \mathbb{E}[c(X, Y)] - \mathbb{E}[c(X, Y)|X = x_i] \right) $$

and

$$ -D_{uu}^2 \Phi(u(0), 0) = \frac{1}{\eta} (-\mathrm{diag}(\bar{\mu}) + \bar{\mu} \bar{\mu}^T), $$

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

$$ \begin{aligned} [-D_{uu}^2 \Phi(u(0), 0)]^{-1} &= \eta \left( -\mathrm{diag}^{-1}(\bar{\mu}) - \frac{\mathrm{diag}^{-1}(\bar{\mu})\bar{\mu}\bar{\mu}^T \mathrm{diag}^{-1}(\bar{\mu})}{1 - \bar{\mu}^T \mathrm{diag}^{-1}(\bar{\mu})\bar{\mu}} \right) \\ &= \eta \left( -\mathrm{diag}^{-1}(\bar{\mu}) - \frac{\mathbf{1}\mathbf{1}^T}{\mu_0} \right), \end{aligned} $$

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

$$ \begin{aligned} & [-D_{uu}^2 \Phi(u(0), 0)]^{-1} \left[ \nabla_u \frac{\partial}{\partial \epsilon} \Phi(u(0), 0) \right] \\ &= \eta \left( -\mathrm{diag}^{-1}(\bar{\mu}) - \frac{\mathbf{1}\mathbf{1}^T}{\mu_0} \right) \left( \frac{\mu_i}{\eta} \left( \mathbb{E}[c(X, Y)] - \mathbb{E}[c(X, Y)|X = x_i] \right) \right)_{i \neq 0} \\ &= \left( \mathbb{E}[c(X, Y)|X = x_i] - \mathbb{E}[c(X, Y)] \right)_{i \neq 0} \\ &\quad + \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} \\ &= \left( \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] \right)_{i \neq 0} \\ &= \left( \mathbb{E}[c(X, Y)|X = x_i] - \mathbb{E}[c(X, Y)|X = x_0] \right)_{i \neq 0} \end{aligned} $$
$\epsilon = 0$, for every $i \neq 0$

$$ \left( \nabla_u \frac{\partial}{\partial \epsilon} \Phi(u(0), 0) \right)_i = \frac{\mu_i}{\eta} \left( \mathbb{E}[c(X, Y)] - \mathbb{E}[c(X, Y)|X = x_i] \right) $$

and

$$ -D_{uu}^2 \Phi(u(0), 0) = \frac{1}{\eta} (-\mathrm{diag}(\bar{\mu}) + \bar{\mu} \bar{\mu}^T), $$

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

$$ \begin{aligned} [-D_{uu}^2 \Phi(u(0), 0)]^{-1} &= \eta \left( -\mathrm{diag}^{-1}(\bar{\mu}) - \frac{\mathrm{diag}^{-1}(\bar{\mu})\bar{\mu}\bar{\mu}^T \mathrm{diag}^{-1}(\bar{\mu})}{1 - \bar{\mu}^T \mathrm{diag}^{-1}(\bar{\mu})\bar{\mu}} \right) \\ &= \eta \left( -\mathrm{diag}^{-1}(\bar{\mu}) - \frac{\mathbf{1}\mathbf{1}^T}{\mu_0} \right), \end{aligned} $$

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

$$ \begin{aligned} & [-D_{uu}^2 \Phi(u(0), 0)]^{-1} \left[ \nabla_u \frac{\partial}{\partial \epsilon} \Phi(u(0), 0) \right] \\ &= \eta \left( -\mathrm{diag}^{-1}(\bar{\mu}) - \frac{\mathbf{1}\mathbf{1}^T}{\mu_0} \right) \left( \frac{\mu_i}{\eta} \left( \mathbb{E}[c(X, Y)] - \mathbb{E}[c(X, Y)|X = x_i] \right) \right)_{i \neq 0} \\ &= \left( \mathbb{E}[c(X, Y)|X = x_i] - \mathbb{E}[c(X, Y)] \right)_{i \neq 0} \\ &\quad + \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} \\ &= \left( \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] \right)_{i \neq 0} \\ &= \left( \mathbb{E}[c(X, Y)|X = x_i] - \mathbb{E}[c(X, Y)|X = x_0] \right)_{i \neq 0} \end{aligned} $$
and

$$
\begin{align*}
& \left[ \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right]^T \left[ -D_{uu}^2 \Phi(u(\varepsilon), \varepsilon) \right]^{-1} \left[ \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right] \\
&= \left( \frac{\mu_i}{\eta} (\mathbb{E}[c(X, Y)] - \mathbb{E}[c(X, Y)|X=x_i]) \right)_{i \neq 0}^T \\
&\quad \cdot (\mathbb{E}[c(X, Y)|X=x_i] - \mathbb{E}[c(X, Y)|X=x_0])_{i \neq 0} \\
&= \frac{1}{\eta} \sum_{i \neq 0} \left[ \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] \right. \\
&\qquad \left. - (\mathbb{E}[c(X, Y)|X=x_i])^2 + \mathbb{E}[c(X, Y)|X=x_0] \mathbb{E}[c(X, Y)|X=x_i] \right] \\
&= \frac{1}{\eta} \left\{ \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 \right. \\
&\qquad \left. - \sum_i (\mathbb{E}[c(X, Y)|X=x_i])^2 \mu_i + \mathbb{E}[c(X, Y)|X=x_0] \sum_i \mathbb{E}[c(X, Y)|X=x_i] \mu_i \right\} \\
&\qquad - \frac{1}{\eta} \left\{ \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] \right. \\
&\qquad \left. - (\mathbb{E}[c(X, Y)|X=x_0])^2 + \mathbb{E}[c(X, Y)|X=x_0] \mathbb{E}[c(X, Y)|X=x_0] \right\} \mu_0 \\
&= \frac{1}{\eta} \left( (\mathbb{E}[c(X, Y)])^2 - \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2] \right). \tag{B.3}
\end{align*}
$$

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

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

B.1. The second order derivatives of Φ(u, v, ε) 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], \quad e_{ijk}^{c} = \exp \left[ \frac{1}{\eta} (u_i + v_j - \varepsilon c_{ijk}) \right] c_{ijk}.
$$
and

$$
\begin{align*}
& \left[ \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right]^T \left[ -D_{uu}^2 \Phi(u(\varepsilon), \varepsilon) \right]^{-1} \left[ \nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), \varepsilon) \right] \\
&= \left( \frac{\mu_i}{\eta} \left( \mathbb{E}[c(X, Y)] - \mathbb{E}[c(X, Y)|X=x_i] \right) \right)_{i \neq 0}^T \\
&\quad \cdot \left( \mathbb{E}[c(X, Y)|X=x_i] - \mathbb{E}[c(X, Y)|X=x_0] \right)_{i \neq 0} \\
&= \frac{1}{\eta} \sum_{i \neq 0} \left[ \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] \right. \\
&\qquad \left. - (\mathbb{E}[c(X, Y)|X=x_i])^2 + \mathbb{E}[c(X, Y)|X=x_0] \mathbb{E}[c(X, Y)|X=x_i] \right] \\
&= \frac{1}{\eta} \left\{ \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 \right. \\
&\qquad \left. - \sum_i (\mathbb{E}[c(X, Y)|X=x_i])^2 \mu_i + \mathbb{E}[c(X, Y)|X=x_0] \sum_i \mathbb{E}[c(X, Y)|X=x_i] \mu_i \right\} \\
&\quad - \frac{1}{\eta} \left\{ \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] \right. \\
&\qquad \left. - (\mathbb{E}[c(X, Y)|X=x_0])^2 + \mathbb{E}[c(X, Y)|X=x_0] \mathbb{E}[c(X, Y)|X=x_0] \right\} \mu_0 \\
&= \frac{1}{\eta} \left( (\mathbb{E}[c(X, Y)])^2 - \mathbb{E}[(\mathbb{E}[c(X, Y)|X])^2] \right). \tag{B.3}
\end{align*}
$$

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

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

B.1. The second order derivatives of Φ(u, v, ε) 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], \quad e_{ijk}^{c} = \exp \left[ \frac{1}{\eta} (u_i + v_j - \varepsilon c_{ijk}) \right] c_{ijk}.
$$
The Hessian matrix of $\Phi$ will be the symmetric matrix:

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

where:

$$\left(D_{uu}^2 \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_{i_1 i_2} = \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$$

$$\left(D_{vv}^2 \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_{j_1 j_2} = \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$$

$$\left(D_{uv}^2 \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_{ij} = \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}, \qquad (B.4)$$

where:

$$\left(\nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_i = -\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$$

$$\left(\nabla_v \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_j = -\frac{1}{\eta} \sum_t \left[ \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} \right] \nu_t.$$

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

Recall:

$$\Phi(v, g, \varepsilon) := - \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_{vv}^2\Phi & D_{vh}^2\Phi \\ D_{hv}^2\Phi & D_{hh}^2\Phi \end{bmatrix},$$
The Hessian matrix of $\Phi$ will be the symmetric matrix:

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

where:

$$ \begin{aligned} \left(D_{uu}^2 \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_{i_1 i_2} &= \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 \\ \left(D_{vv}^2 \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_{j_1 j_2} &= \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 \\ \left(D_{uv}^2 \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_{ij} &= \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 \end{aligned} $$

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}, \qquad (B.4)$$

where:

$$ \begin{aligned} \left(\nabla_u \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_i &= -\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 \\ \left(\nabla_v \frac{\partial}{\partial \varepsilon} \Phi(u(\varepsilon), v(\varepsilon), \varepsilon)\right)_j &= -\frac{1}{\eta} \sum_t \left[ \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} \right] \nu_\ell. \end{aligned} $$

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

Recall:

$$\Phi(v, g, \varepsilon) := - \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_{vv}^2\Phi & D_{vh}^2\Phi \\ D_{hv}^2\Phi & D_{hh}^2\Phi \end{bmatrix},$$
where:

$$
\begin{align*}
\left( D_{vv}^2 \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_{j_1 j_2} &= \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 \\
\left( D_{hh}^2 \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_{i_1 i_2} &= \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} \\
\left( D_{vh}^2 \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_{i j} &= \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,
\end{align*}
$$

where $\delta_{il}$ 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:

$$
\begin{align*}
\left( \nabla_v \frac{\partial}{\partial \varepsilon} \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_j &= -\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 \\
\left( \nabla_h \frac{\partial}{\partial \varepsilon} \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_i &= -\frac{1}{\eta} \left[ \frac{\sum_s e_{is}^{yc}}{\sum_s e_{is}} - \frac{\sum_s e_{is}^{xy} \cdot \sum_s e_{is}^c}{(\sum_s e_{is})^2} \right] \mu_i.
\end{align*}
$$

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

Recall:

$$
\Phi(v, w, g, h, \varepsilon) = - \sum_s v_s \theta_s - \sum_t w_t \nu_t + \eta \sum_r \log \left[ \sum_{s,t} \exp \left[ \frac{1}{\eta} (v_s + w_t + g_r(y_s - x_r) + h_{rs}(z_t - y_s) - \varepsilon c_{rst}) \right] \theta_s \nu_t \right] \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}^{yc} = 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_{vv}^2 \Phi & D_{vw}^2 \Phi & D_{vg}^2 \Phi & D_{vh}^2 \Phi \\ D_{wv}^2 \Phi & D_{ww}^2 \Phi & D_{wg}^2 \Phi & D_{wh}^2 \Phi \\ D_{gv}^2 \Phi & D_{gw}^2 \Phi & D_{gg}^2 \Phi & D_{gh}^2 \Phi \\ D_{hv}^2 \Phi & D_{hw}^2 \Phi & D_{hg}^2 \Phi & D_{hh}^2 \Phi \end{bmatrix},
$$
where:

$$
\begin{align*}
\left( D_{vv}^2 \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_{j_1 j_2} &= \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 \\
\left( D_{hh}^2 \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_{i_1 i_2} &= \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} \\
\left( D_{vh}^2 \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_{i j} &= \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,
\end{align*}
$$

where $\delta_{il}$ 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:

$$
\begin{align*}
\left( \nabla_v \frac{\partial}{\partial \varepsilon} \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_j &= -\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 \\
\left( \nabla_h \frac{\partial}{\partial \varepsilon} \Phi(v(\varepsilon), h(\varepsilon), \varepsilon) \right)_i &= -\frac{1}{\eta} \left[ \frac{\sum_s e_{is}^{yc}}{\sum_s e_{is}} - \frac{\sum_s e_{is}^{xy} \cdot \sum_s e_{is}^c}{(\sum_s e_{is})^2} \right] \mu_i.
\end{align*}
$$

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

Recall:

$$
\Phi(v, w, g, h, \varepsilon) = - \sum_s v_s \theta_s - \sum_t w_t \nu_t + \eta \sum_r \log \left[ \sum_{s,t} \exp \left[ \frac{1}{\eta} (v_s + w_t + g_r(y_s - x_r) + h_{rs}(z_t - y_s) - \varepsilon c_{rst}) \right] \theta_s \nu_t \right] \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}^{yc} = 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_{vv}^2 \Phi & D_{vw}^2 \Phi & D_{vg}^2 \Phi & D_{vh}^2 \Phi \\
D_{wv}^2 \Phi & D_{ww}^2 \Phi & D_{wg}^2 \Phi & D_{wh}^2 \Phi \\
D_{gv}^2 \Phi & D_{gw}^2 \Phi & D_{gg}^2 \Phi & D_{gh}^2 \Phi \\
D_{hv}^2 \Phi & D_{hw}^2 \Phi & D_{hg}^2 \Phi & D_{hh}^2 \Phi
\end{bmatrix},
$$
where

$$ (D_{vv}^2 \Phi)_{j_1 j_2} = \frac{1}{\eta} \sum_r \left[ \frac{\sum_t e_{rj_1t}}{\sum_{s,t} e_{rst}} \delta_{j_1 j_2} - \frac{\sum_t e_{rj_1t} \cdot \sum_t e_{rj_2t}}{(\sum_{s,t} e_{rst})^2} \right] \mu_r $$

$$ (D_{vw}^2 \Phi)_{jk} = \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 $$

$$ (D_{vg}^2 \Phi)_{ji} = \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 $$

$$ (D_{vh}^2 \Phi)_{j_1, (ij_2)} = \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 $$

$$ (D_{ww}^2 \Phi)_{k_1 k_2} = \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 $$

$$ (D_{wg}^2 \Phi)_{ki} = \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 $$

$$ (D_{wh}^2 \Phi)_{k, (ij)} = \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 $$

$$ (D_{gg}^2 \Phi)_{i_1 i_2} = \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} $$

$$ (D_{gh}^2 \Phi)_{i_1, (ij_2)} = \frac{1}{\eta} \left[ \frac{\sum_t e_{i_1 j_1 t}^{xyz}}{\sum_{s,t} e_{ist}} - \frac{\sum_{s,t} e_{i_1 st}^{xy} \cdot \sum_t e_{i_2 j t}^{yz}}{(\sum_{s,t} e_{ist})^2} \right] \mu_{i_1} \delta_{i_1 i_2} $$

$$ (D_{hh}^2 \Phi)_{(i_1 j_1), (i_2 j_2)} = \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

$$ (D_{vv}^2 \Phi)_{j_1 j_2} = \frac{1}{\eta} \sum_r \left[ \frac{\sum_t e_{rj_1t}}{\sum_{s,t} e_{rst}} \delta_{j_1 j_2} - \frac{\sum_t e_{rj_1t} \cdot \sum_t e_{rj_2t}}{(\sum_{s,t} e_{rst})^2} \right] \mu_r $$

$$ (D_{vw}^2 \Phi)_{jk} = \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 $$

$$ (D_{vg}^2 \Phi)_{ji} = \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 $$

$$ (D_{vh}^2 \Phi)_{j_1, (ij_2)} = \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 $$

$$ (D_{ww}^2 \Phi)_{k_1 k_2} = \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 $$

$$ (D_{wg}^2 \Phi)_{ki} = \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 $$

$$ (D_{wh}^2 \Phi)_{k, (ij)} = \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 $$

$$ (D_{gg}^2 \Phi)_{i_1 i_2} = \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} $$

$$ (D_{gh}^2 \Phi)_{i_1, (ij_2)} = \frac{1}{\eta} \left[ \frac{\sum_t e_{i_1 j_1 t}^{xyz}}{\sum_{s,t} e_{ist}} - \frac{\sum_{s,t} e_{i_1 st}^{xy} \cdot \sum_t e_{i_2 j t}^{yz}}{(\sum_{s,t} e_{ist})^2} \right] \mu_{i_1} \delta_{i_1 i_2} $$

$$ (D_{hh}^2 \Phi)_{(i_1 j_1), (i_2 j_2)} = \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

$$ \left(\nabla_v \frac{\partial}{\partial \varepsilon} \Phi\right)_j = -\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 $$

$$ \left(\nabla_g \frac{\partial}{\partial \varepsilon} \Phi\right)_i = -\frac{1}{\eta} \left[ \frac{\sum_{s,t} e_{ist}^{yc}}{\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 $$

$$ \left(\nabla_w \frac{\partial}{\partial \varepsilon} \Phi\right)_k = -\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 $$

$$ \left(\nabla_h \frac{\partial}{\partial \varepsilon} \Phi\right)_{(ij)} = -\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.
where

$$ \left(\nabla_v \frac{\partial}{\partial \varepsilon} \Phi\right)_j = -\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 $$

$$ \left(\nabla_g \frac{\partial}{\partial \varepsilon} \Phi\right)_i = -\frac{1}{\eta} \left[ \frac{\sum_{s,t} e_{ist}^{yc}}{\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 $$

$$ \left(\nabla_w \frac{\partial}{\partial \varepsilon} \Phi\right)_k = -\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 $$

$$ \left(\nabla_h \frac{\partial}{\partial \varepsilon} \Phi\right)_{(ij)} = -\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.
[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.
[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.