# Entropic Semi-Martingale Optimal Transport

Jean-David Benamou, Guillaume Chazareix, Marc Hoffmann, Gregoire

Loeper, François-Xavier Vialard

► To cite this version:

Jean-David Benamou, Guillaume Chazareix, Marc Hoffmann, Gregoire Loeper, François-Xavier Vialard. Entropic Semi-Martingale Optimal Transport. 2024. hal-04673273v2

HAL Id: hal-04673273

https://hal.science/hal-04673273v2

Preprint submitted on 16 Dec 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
# Entropic Semi-Martingale Optimal Transport

Jean-David Benamou<sup>1,2*</sup>, Guillaume Chazareix<sup>2,3</sup>,
Marc Hoffmann<sup>1,4</sup>, Grégoire Loeper<sup>3</sup>, François-Xavier Vialard<sup>5</sup>

<sup>1</sup>CEREMADE, Université Paris Dauphine-PSL.

<sup>2</sup>MOKAPLAN, INRIA Paris.

<sup>3</sup>BNP Paribas Global Markets.

<sup>4</sup>Institut Universitaire de France.

<sup>5</sup>LIGM, Université Gustave Eiffel.

*Corresponding author(s). E-mail(s): Jean-David.Benamou@inria.fr;

Contributing authors: Guillaume.Chazareix@inria.fr;

hoffmann@ceremade.dauphine.fr; gregoire.loeper@bnpparibas.com;

francois-xavier.vialard@univ-eiffel.fr;

## Abstract

Entropic Optimal Transport (EOT), also known as the Schrödinger problem, aims to find random processes with given initial and final marginals, minimizing the relative entropy (RE) with respect to a reference measure. Both processes (the reference and the controlled one) necessarily share the same diffusion coefficients to ensure finiteness of the RE. This initially suggests that controlled-diffusion Semi-Martingale Optimal Transport (SMOT) problems may be incompatible with entropic regularization. However, when the process is observed at discrete times, forming a Markov chain, the RE remains finite even with variable diffusion coefficients. In this case, discrete semi-martingales can emerge as solutions to multi-marginal EOT problems. For smooth semi-martingales, the scaled limit of the relative entropy of their time discretizations converges to the “specific relative entropy”, a convex functional of the variance. This observation leads to an entropy regularized time discretization of the continuous SMOT problems, enabling the computation of discrete approximations via a multi-marginal Sinkhorn algorithm. We prove convergence of the time-discrete entropic problem to the continuous case, provide an implementation, and present numerical experiments supporting the theoretical results.

**Keywords:** Semi-Martingale Optimal Transport, Multi-Marginal Optimal Transport, Entropic Penalisation, Specific Relative Entropy, Sinkhorn Algorithm
<table><tr><td>4</td><td>Dual formulation and Sinkhorn algorithm</td><td>12</td></tr><tr><td>4.1</td><td>Fenchel-Rockafellar duality</td><td>12</td></tr><tr><td>4.2</td><td>Multi-Marginal Sinkhorn Algorithm</td><td>14</td></tr><tr><td>4.3</td><td>Numerical experiments</td><td>16</td></tr><tr><td>5</td><td>Appendix</td><td>23</td></tr><tr><td>5.1</td><td>Proof of proposition 2 [Specific Relative Entropy]</td><td>23</td></tr><tr><td>5.2</td><td>Proof of theorem 10 [Convergence of Markov Chain to diffusion]</td><td>27</td></tr><tr><td>5.3</td><td>Convergence of Markov chains to diffusions</td><td>27</td></tr><tr><td>5.4</td><td>Proof of Theorem 3 [Well posedness for the discrete problem]</td><td>29</td></tr><tr><td>5.5</td><td>Proof of Lemma 5 [Diffusion Coefficients rescaling]</td><td>30</td></tr><tr><td>5.6</td><td>Proofs of Lemma 6 [Regularization time continuous case]</td><td>32</td></tr><tr><td>5.7</td><td>Proof of Lemma 7</td><td>32</td></tr><tr><td>5.8</td><td>Proof of Lemma 8 [Regularization time discrete case]</td><td>33</td></tr><tr><td>5.9</td><td>Additional Lemmas</td><td>36</td></tr></table>

# 1 Introduction

The distinction between Entropy-Regularized Optimal Transport (EOT) and Semi-Martingale Optimal Transport (SMOT), further discussed in sections 2.1 and 2.2, hinges on the underlying stochastic processes being optimized. EOT [1-4] is a well-established method in computational optimal transport, equivalent to relative entropy (RE) minimization with respect to a given reference diffusion, typically the Wiener measure. To ensure the finiteness of the RE, both the reference and controlled processes must have the same diffusion coefficients. This characteristic makes EOT particularly suitable for applications where controlling the drift is adequate and computational efficiency is paramount, such as in machine learning, statistics, and physics.

In contrast, SMOT [5] introduces greater complexity by allowing control over both the drift and the diffusion coefficients. This expanded control framework enables one to model more intricate dynamics that cannot be adequately captured by a drift control alone, with model calibration in finance serving as a prime example.

When discretizing continuous diffusions processes in time and formulating them as Markov chains, the RE may remain finite even for different diffusion coefficients that imply orthogonal laws on the space of continuous trajectories. Notably, as the time discretization step approaches zero, the asymptotic behavior of the discrete RE can be analyzed. Properly scaled by the time step, the first-order term converges to a well-defined limit known in the literature as the specific relative entropy (SRE), which represents a divergence between the respective diffusion coefficients. This concept was apparently initially proposed in [6].

In section 3.1, we use SRE as a regularization term for the time-continuous SMOT problem. It can also be interpreted as a fidelity term pertaining to the diffusion coefficient of the optimal associated semimartingale. The natural time discretization optimizes over Markov chains, and the SRE becomes a “correctly” scaled discrete RE regularization. In this framework, arbitrary time discretizations of semimartingales can be realized as solutions to multi-marginal EOT problems. Our first contribution, Theorem 4, demonstrates that the sequence of Markov chains indexed by the time step converges weakly (similar to Gamma-convergence) to a minimizer of the time-continuous ESMOT problem (V0). The convergence proof is sketched in section 3.2, with detailed elaboration provided in Appendix 5.

The numerical solution relies on the dual formulation of the problem, as presented in equations (22-24-25). We formally address the duality in section 4.1, where we apply Fenchel-Rockafellar convex duality, common in optimal transport, by using a linear change of variables
involving moments rather than conditional moments. The integrand is interpreted as a convex lower semicontinuous perspective function (45). Our implementation is based on the Sinkhorn algorithm, which is detailed in section 4.2. While the convergence analysis of Sinkhorn in both continuous and discrete cases remains to be thoroughly examined, we present the implementation and several test case experiments in section 4.3.

Our novel regularization approach illustrates that EOT can be extended to more general diffusion processes, bridging the classical EOT framework and broader SMOT problems, thus giving rise to the term Entropic Semi-Martingale Optimal Transport (ESMOT). The significance of this new approach is not solely theoretical. Traditional methods often involve maximizing the dual problem through gradient ascent or primal-dual techniques, which require solving a fully nonlinear Hamilton-Jacobi-Bellman equation at each iteration [7] and [8, 9]. In contrast, our approach using the Multi-Marginal Sinkhorn algorithm, which builds on [10] and [11], is straightforward to implement and comes with convergence guarantees. This method has been successfully applied to model calibration using derivative price observations [12], yielding smoother volatility surfaces compared to alternative parametric or PDE-based methods. Potential applications are probably not limited to finance: for instance, our approach could improve on the recent OT modeling for genes' developmental trajectories of cells [13].

There are very few related works in the literature, [14] highlights the challenges of using relative entropy between singular measures and discusses a discrete version of SRE. However, this work remains entirely at the time discrete level and does not connect with EOT. Additionally, the characterization of specific relative entropy for smooth Itô diffusions can also be found in [15] and [16], where the authors apply time discretization and scaling techniques, substituting relative entropy with p-Wasserstein distance and investigating the properties of this new divergence between martingales.

# 2 Background

## 2.1 Entropic Optimal Transport (EOT)

EOT is a well-studied problem. Given an initial distribution $\rho_0$ on $\mathbb{R}^d$ and $\bar{\sigma}^2 > 0$, call $\bar{P}$ the distribution of a Brownian motion with variance $\bar{\sigma}^2$ and initial random condition distributed according to $\rho_0$. Call $\bar{P}_{0,1}$ its joint (initial, final) distribution at times 0 and 1. The EOT problem can be formulated in a static form as

$$ \inf_{\bar{P}_{0,1} \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d), \bar{P}_0 = \rho_0, \bar{P}_1 = \rho_1} \bar{\sigma}^2 \mathcal{H}(\bar{P}_{0,1} | \bar{P}_{0,1}). \quad (1) $$

The equivalent time continuous dynamic form version of (1), known as the Schrödinger problem, is

$$ \inf_{\bar{P} \in \mathcal{P}, \bar{P}_0 = \rho_0, \bar{P}_1 = \rho_1} \bar{\sigma}^2 \mathcal{H}(\bar{P} | \bar{P}). \quad (2) $$

The equivalence is easily seen from the additive property of the relative entropy applied to the disintegration of the measures $\bar{P}$ and $\bar{P}$

$$ \mathcal{H}(\bar{P} | \bar{P}) = \mathcal{H}(\bar{P}_{0,1} | \bar{P}_{0,1}) + \mathbb{E}_{\bar{P}_{0,1}} (\mathcal{H}(\bar{P}(\cdot | X_0, X_1) | \bar{P}(\cdot | X_0, X_1))). $$

In particular, the optimal $\bar{P}$ and $\bar{P}$ necessarily share the same point-to-point Brownian bridges. It is well-known (see [2] for a survey) that the minimizers of the Schrödinger problem are semi-martingales with fixed volatility $\bar{\sigma}$ and local drift converging to the OT map as $\bar{\sigma} \to 0$. There is a large literature based on this property as, unlike the classical OT problem, the EOT is a strictly convex optimization problem and can be solved efficiently through the celebrated Sinkhorn algorithm (see [1, 3]). The Entropic regularisation and Sinkhorn solving approaches are not limited to the static two-marginals problem (1). They can be applied to multimarginal
measures. Let $\mathbb{P}$ be the law of a continuous semi-martingale, and $\mathbb{P}^h = (e_{t_0, t_1, \dots, t_N})_{\#\mathbb{P}}$ its time discretization (and using a similar notation for $\bar{\mathbb{P}}$). Remark that even tough $\mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h)$ is finite, a discete time multi-marginal entropic $\bar{\sigma}^2 \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h)$ regularization has the strong consequence that in the limit $h \searrow 0$, $\mathbb{P}^h$ must have asymptotically the same volatility $\bar{\sigma}$ otherwise the continuous time entropy $\mathcal{H}(\mathbb{P} | \bar{\mathbb{P}})$ is infinite. It was used to optimize over constant prescribed diffusion model in [11, 17] in particular.

## 2.2 Semi Martingale OT (SMOT)

SMOT allows for local in time and space diffusion coefficients. Our goal is to extend the Entropic regularisation approach in this setting.

We closely follow the framework proposed by [5] and also in [18]. On the path space $\Omega = C^0([0, 1], \mathbb{R}^d)$, we consider the canonical process $X_t(\omega) = \omega(t)$ and its associated canonical filtration, $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \in [0, 1]})$. Let $\mathcal{P}^0$ be the set of probability measures $\mathbb{P}$ adapted to $\mathcal{F}$, such that $X$ is an $(\mathcal{F}, \mathbb{P})$-semimartingale on $[0, 1]$. Thus, $X_t$ can be decomposed as

$$X_t = X_0 + B_t^\mathbb{P} + M_t^\mathbb{P}, \quad \langle X \rangle_t = \langle M^\mathbb{P} \rangle_t = A_t^\mathbb{P},$$

where $M^\mathbb{P}$ is an $(\mathcal{F}, \mathbb{P})$-martingale on $[0, 1]$. In addition, $(B^\mathbb{P}, A^\mathbb{P})$ is $\mathcal{F}$-adapted and $\mathbb{P}$-as absolutely continuous with respect to time. In particular, $\mathbb{P}$ is said to have characteristics $(a^\mathbb{P}, b^\mathbb{P})$, which are defined in the following way,

$$a_t^\mathbb{P} = \frac{dA_t^\mathbb{P}}{dt}, \quad b_t^\mathbb{P} = \frac{dB_t^\mathbb{P}}{dt}.$$

Note that $(a^\mathbb{P}, b^\mathbb{P})$ is $\mathcal{F}$-adapted and determined $d\mathbb{P} \otimes dt$ almost everywhere. We might use $a, b$ instead of $a^\mathbb{P}, b^\mathbb{P}$ for simplicity, but it is fundamental to note that the object we are optimising is the measure $\mathbb{P}$, from which $a, b$ are derived. Letting $S^d$ be the set of symmetric matrices and $S_+^d$ the set of positive semidefinite matrices, $(a^\mathbb{P}, b^\mathbb{P})$ takes values in the space $S_+^d \times \mathbb{R}^d$. For $a_1, a_2 \in S_+^d$, we define $a_1 : a_2 := \text{trace}(a_1^\top a_2)$.

Denote by $\mathcal{P}^1 \subset \mathcal{P}^0$ the set of probability measures $\mathbb{P}$ whose characteristics $(a^\mathbb{P}, b^\mathbb{P})$ are $\mathbb{P} \otimes dt$-integrable:

$$\mathbb{E} \left( \int_0^1 |a^\mathbb{P}| + |b^\mathbb{P}| dt \right) < \infty.$$

Recall that $\mathbb{P}_t = X_t \#\mathbb{P}$, denotes the marginal of $\mathbb{P}$ at time $t$.

The SMOT stochastic control problem (see [5, 18–20]) takes the form

$$\inf_{\mathbb{P} \in \mathcal{P}^1} \mathbb{E}_{\mathbb{P}} \left( \int_0^1 F(t, X_t, b_t^\mathbb{P}, a_t^\mathbb{P}) dt \right) + \mathcal{D}(\mathbb{P}_0, \rho_0) + \mathcal{D}(\mathbb{P}_1, \rho_1), \quad (3)$$

where $F$ has to be convex in $(a, b)$ at every $(t, x)$ and also enforces the non-negativity of $a$ the diffusion coefficient, $\mathcal{D}(\mathbb{P}_t, \rho)$ is convex in $\mathbb{P}_t$, for a given $\rho$. Typically $\mathcal{D}$ will be a penalty function to enforce that $\mathbb{P}_t$ is either equal or close to $\rho$ at initial and final times. We will use hard constraints as $\iota_{\rho_t}(\mathbb{P}_t)$, as in the usual optimal transport problem in the numerical section. Still, the convergence theorem relies on a soft constraint for instance $W_2(\mathbb{P}_t, \rho_t)$ the quadratic Wasserstein distance. We will assume that there exists a positive constant $C$ such that for every pair $\rho_1, \rho_2$ (denoting $W_2(\rho_0, \rho_1)$ the classical quadratic Wasserstein distance)

**Hypothesis H1 (H1).** We assume the following hypotheses on $\mathcal{D}$

1. $\mathcal{D}$ is strictly convex on its domain,

2. $\forall \rho_1, \rho_2 \in \mathcal{P}(\mathbb{R}^d)$, such that $\mathcal{D}(\rho_1, \rho_2) < +\infty$, then $\int |x|^4 d\rho_1(x) < \infty$,
3. $\forall \rho_2 \in \mathcal{P}(\mathbb{R}^d)$, there exists $C(\rho_2) > 0$ such that $|\mathcal{D}(\rho_1, \rho_2) - \mathcal{D}(\rho_3, \rho_2)| \le CW_2(\rho_1, \rho_3)$, when all the quantities are finite.

As a consequence, if $\eta$ is a space convolution kernel with variance $v_\eta$, then $|\mathcal{D}(\eta \star \rho_1, \rho_2) - \mathcal{D}(\eta \star \rho_1, \rho_2)| \le Cv_\eta$ when $\mathcal{D}(\eta \star \rho_1, \rho_2) < +\infty$.

**Example 1.** These constraints are satisfied by the family of discrepancies indexed by $M > 0$,

$$ \mathcal{D}_M(\rho_1, \rho_2) = \iota_{C(M)}(\rho_1) + W_2(\rho_1, \rho_2), \qquad (4) $$

where $\iota_{C(M)}$ is the convex indicator function of the convex set $C(M) := \{\rho \in \mathcal{P}(\mathbb{R}^d), \int |x|^4 d\rho(x) < M\}$.

**Remark 1.** • Note that the discrepancy $\mathcal{D}$ can be nonsymmetric. In the first argument, we will consider the marginals of the probability measure on the path space on which the optimisation is done. In the second argument, we consider the data which can be any probability measure.

• Assuming $\rho_0$ and $\rho_1$ are in convex order, time-dependent martingale constrained processes (see for example [21]) correspond to (3) with $F = F(a)$ if $b \equiv 0, +\infty$ otherwise. In particular, the Bass Martingale Problem (see [22] for example) corresponds to the choice $F = |a - Id|^2$.

• The Schrödinger bridge problem (2) corresponds to $F = |b|^2$ if $a \equiv \varepsilon, +\infty$ otherwise, for a positive constant $\varepsilon$.

Usual coercivity assumptions on $F$ are

$$ F \ge C(1 + |a|^p + |b|^p) \quad \text{for } C > 0, p > 1. \qquad (5) $$

We will restrict to a strongly coercive case:

$$ F = G + \iota_{[\lambda, \Lambda]}(a) + \iota_{[0, B]}(|b|) \qquad (H2) $$

with $\lambda, \Lambda, B$ given positive parameters, and $G$ satisfying

$$ G \text{ is strictly convex and Lipschitz on } [0, B] \times [\lambda, \Lambda]. \qquad (H3) $$

The bounds in (H2) are technical assumptions that allow us to later handle the convergence of the discrete entropy to the specific entropy.

As explained in [20], when $F$ is a function of $(t, X, b, a)$ the optimal $a$ and $b$ will automatically be functions of $(t, X)$, even if one does not restrict such local processes in the minimisation procedure. This is due to the convexity of $F$ which favors local functions. We state here a general and standard existence result for (3), that can be found in [18]:

**Theorem 1.** Under conditions (H1, H2, H3), if there exists an admissible solution (i.e. a probability $P$ with finite energy in (3)) there exists a minimiser to (3).

This result follows from standard arguments on the compactness of any minimizing sequence, and the lower semi-continuity of the energy functional.

## 2.3 Time discretisation and specific relative entropy

The concept of specific relative entropy between diffusion processes seems to go back to [6], see also [23] and [15].

Let us consider $P \in \mathcal{P}^1$ and the reference measure $\bar{P}$ as section 2.1. We denote by $\mathbb{P}^h = (e_{t_1, t_2, \dots})_{\mathbb{P}}$ and $\bar{\mathbb{P}}^h = (e_{t_1, t_2, \dots})_{\bar{\mathbb{P}}}$ their time discretisation with time mesh $h$ (see the notation section). Then under regularity conditions on the drift and diffusion coefficients
specified in Proposition 2 below, the specific relative entropy between $\mathbb{P}$ and $\bar{\mathbb{P}}$ is defined as the limit:

$$S^{\mathcal{H}}(\mathbb{P}|\bar{\mathbb{P}}) := \lim_{h \searrow 0} h \mathcal{H}(\mathbb{P}^h|\bar{\mathbb{P}}^h).$$

Under additional smoothness and boundedness assumptions on the characteristics $b^{\mathbb{P}}, a^{\mathbb{P}}$, the specific entropy can be characterized (see proposition 2) explicitly:

$$S^{\mathcal{H}}(\mathbb{P}|\bar{\mathbb{P}}) = \lim_{h \searrow 0} h \mathcal{H}(\mathbb{P}^h|\bar{\mathbb{P}}^h) = \frac{1}{2} \mathbb{E}_{\mathbb{P}} \left( \int_0^1 S^{\mathcal{I}}(a_t(X_t)|\bar{a}) dt \right). \quad (6)$$

The integrand $S^{\mathcal{I}}$, is given (in dimension 1) by

$$S^{\mathcal{I}}(a|\bar{a}) := \frac{a}{\bar{a}} - 1 - \log \left( \frac{a}{\bar{a}} \right). \quad (7)$$

In the case of a positive definite symmetric matrix in dimension $d$, the definition becomes

$$S^{\mathcal{I}}(a|\bar{a}) := \mathrm{Tr}(\bar{a}^{-1}(a - \bar{a})) - \log \left( \frac{\det(a)}{\det \bar{a}} \right).$$

This function is strictly convex with minimum at $\bar{a}$, a barrier for vanishing $a$ and strictly increasing but just sub-linearly as $a \to +\infty$.

Thanks to the discretisation and renormalisation by $h$, we will see that the limit (6) is well defined even when the diffusions matrices $a$ and $\bar{a}$ differ, in which case the probability measures $\mathbb{P}$ and $\bar{\mathbb{P}}$ are singular as probability measures on the space of continuous functions, which entails that the relative entropy blows up. This is easily understood in the following simplified setting: instead of discretizing a time continuous process, directly assume that $\mathbb{P}^h$ and $\bar{\mathbb{P}}^h$ are a Markov chains with transistions following normal distribution with space dependent $(\mu_i^h, v_i^h)$ coefficients such that:

$$\mathbb{P}_{i \to i+1}^h = \mathcal{N}(x_i + h \mu_i^h(x_i), h v_i^h(x_i))$$

and likewise $\bar{\mathbb{P}}_{i \to i+1}^h = \mathcal{N}(x_i, h \bar{a})$. A direct computation expresses the discrete relative entropy as

$$h \mathcal{H}(\mathbb{P}^h|\bar{\mathbb{P}}^h) = h \mathcal{H}(\mathbb{P}_0^h|\rho_0) + \frac{h}{2} \sum_{i=0}^{N-1} \mathbb{E}_{\mathbb{P}_i^h} (S^{\mathcal{I}}(v_i^h(X_i^h)|\bar{a}) + \frac{h}{\bar{a}} \|\mu_i^h(X_i^h)\|^2), \quad (8)$$

where $(X_i^h)_{i=0..N}$ is the canonical discrete time process associated to $\mathbb{P}^h$. Assuming convergence as $h \searrow 0$ of the piecewise constant in time interpolation of the coefficients and without the $h$ scaling in (8) the zero-order term in (8) blows up consistently with the definition of relative entropy between singular continuous-time diffusion processes. With the scaling, we recover Definition 6.

Considering now an arbitrary Markov chain $\mathbb{P}^h$ we define

$$\left\{ \begin{aligned} b_i^h(x_i) &:= \frac{1}{h} \mathbb{E}_{\mathbb{P}_{i \to i+1}^h} [X_{i+1}^h - x_i], \\ a_i^h(x_i) &:= \frac{1}{h} \mathbb{E}_{\mathbb{P}_{i \to i+1}^h} [(X_{i+1}^h - x_i)(X_{i+1}^h - x_i)^*]. \end{aligned} \right. \quad (9)$$

Remark that compared to the previous Gaussian example $b = \mu$ but $a = v + hb^2$. This choice of characteristics is only equivalent when $h \to 0$, but has the fundamental property that $\mathbb{P}_i^h a_i^h, \mathbb{P}_i^h b_i^h$ are linear quantities with respect to $\mathbb{P}^h$. Using this definition we can show in particular that the discrete relative entropy controls the discrete version of the specific relative entropy, we show (Proposition 2 ii)) that

$$h \mathcal{H}(\mathbb{P}^h|\bar{\mathbb{P}}^h) \ge h \mathcal{H}(\mathbb{P}_0^h|\rho_0) + \frac{h}{2} \sum_{i=0}^{N-1} \mathbb{E}_{\mathbb{P}_i^h} (S^{\mathcal{I}}(a_i^h(X_i^h)|\bar{a})). \quad (10)$$

We gather these results in:
**Proposition 2** (Specific Relative Entropy). i) Let $\mathbb{P} \in \mathcal{P}_1$, see (3), assume that both functions $(t, x) \mapsto b_t(x)$ and $(t, x) \mapsto a_t(x)$ are twice continuously differentiable, and that the matrix $a_t(x)$ is invertible for every $(t, x) \in [0, 1] \times \mathbb{R}^d$. Assume further that $(b, a)$ are bounded above and a below by a positive constant (in the sense of strong ellipticity for a). Then (6) holds.

ii) Let $\mathbb{P}^h$ be the law of a discrete-time Markov chain and $(b_i^h, a_i^h)$ given by (9). We assume further that $d\mathbb{P}^h = \rho.d\bar{\mathbb{P}}^h$ with a continuous density $\rho$, then (10) holds.

*Proof.* See appendix 5.1. □

# 3 Entropic Semi-Martingale Optimal Transport

## 3.1 ESMOT primal formulation and discretization

We introduce the $\mathcal{S}^\mathcal{I}$-regularized time continuous functional

$$ \mathcal{I}^0(\mathbb{P}) := \mathbb{E}_{\mathbb{P}} \left( \int_0^1 F(b_t^\mathbb{P}(X_t), a_t^\mathbb{P}(X_t)) + \mathcal{S}^\mathcal{I}(a_t^\mathbb{P}(X_t)|\bar{a}) dt \right) + \mathcal{D}(\mathbb{P}_0, \rho_0) + \mathcal{D}(\mathbb{P}_1, \rho_1) \quad (11) $$

where $\sqrt{\bar{a}}$ is a given reference volatility target, $b^\mathbb{P}$, $a^\mathbb{P}$ are the characteristic coefficients associated to $\mathbb{P} \in \mathcal{P}^1$ (see section 2.2). Note that there are, formally at least, no obstacles to considering a local in time and space coefficient $\bar{a}$. The fidelity terms $\mathcal{D}(\mathbb{P}_t, \rho)$ properties are summarized in (H1).

The time continuous SRE regularisation of (3) is:

$$ \inf_{\mathbb{P} \in \mathcal{P}^1} \mathcal{I}^0(\mathbb{P}). \qquad (\mathcal{V}^0) $$

Penalization/regularisation of the volatility is often necessary to deal with underdetermined problems like volatility calibration of pricing models in finance (see [14, 18, 19, 24, 25]). Here, the choice of the SRE regularisation is also driven by its discrete-time formulation, which will allow to rely on a Sinkhorn-like algorithm for numerical resolution. Indeed, based on (6) and Theorem 10 we will use the natural time discretisation of (11):

$$ \mathcal{I}^h(\mathbb{P}^h) := h \sum_{i=0}^{N-1} \mathbb{E}_{\mathbb{P}_i^h} \left( F(b_i^h(X_i^h), a_i^h(X_i^h)) \right) + h \mathcal{H}(\mathbb{P}^h|\bar{\mathbb{P}}^h) + \mathcal{D}(\mathbb{P}_0^h, \rho_0) + \mathcal{D}(\mathbb{P}_1^h, \rho_1), \quad (12) $$

with $(b_i^h, a_i^h)$, the discrete drift and quadratic variation increments defined in (9). We thus consider

$$ \left\{ \begin{array}{ll} \mathbb{P}^h \in \mathcal{P}(\bigotimes_{i=0}^{N-1} \mathbb{R}_{t_i}) \text{ s.t.} & \mathcal{I}^h(\mathbb{P}^h). \qquad (\mathcal{V}^h) \\ h \sum_{i=0}^{N-1} \mathbb{E}_{\mathbb{P}_i^h} \left( c_i^h(X_i^h) \right) \le M & \end{array} \right. $$

In this formulation, we have introduced an additional constraint on a moment of order larger than 4 ($\alpha > 0$) referred to abusively as *kurtosis* (kurtosis is $\alpha = 0$, see also (42):

$$ c_i^h(x_i) := h^{-(2+\alpha)} \mathbb{E}_{\mathbb{P}_{i \to i+1}^h} (\|X_{i+1} - x_i\|^{4+2\alpha}). $$

This additional constraint seems necessary to apply Theorem 10 and ensures that the sequence $(\mathbb{P}^h)_h$ converges in $\mathcal{P}^1$. Controlling the specific entropy through (10) is not sufficient to guarantee (as in [26] Chap.11 or [27]) that the limit of $(\mathbb{P}^h)_h$ is a jump free diffusion process. We start with an existence result for the discretised problem:
**Theorem 3.** Under the hypotheses (H1-H2-H3) problem ($V^h$) has a unique Markovian solution $\mathbb{P}^h$.

*Proof.* Appendix 5.4 □

Before studying the convergence of the discretized problem to the continuous one, let us first comment on the difference between the two formulations: there is no explicit bound on the kurtosis of the continuous process. This bound comes from the boundedness of the drift and diffusion coefficients as an immediate consequence of Lemma 11, section 5.9. Therefore, choosing the constant $M$ large enough in the discretised problem will be sufficient to ensure that the continuous limit of this constant is not saturated.

Our main result is the convergence of the values of the discretized problems to the continuous one:

**Theorem 4.** Let $\rho_0, \rho_1 \in \mathcal{P}_2(\mathbb{R}^d)$ be given. Under the hypotheses (H1-H2-H3), let $\mathbb{P}^0$ be a minimizer of ($V^0$) and ($\mathbb{P}^h$)$_h$ a sequence of minimizers of ($V^h$). Then,

$$ \lim_{h \searrow 0} \mathcal{I}^h(\mathbb{P}^h) = \mathcal{I}^0(\mathbb{P}^0). $$

*Proof.* See Section 3.2. □

**Remark 2.** Theorem 4 does not elaborate on the convergence of the minimizing sequences $(b^h, a^h, \mathbb{P}^h)$, but in the course of the proof we construct from $(\mathbb{P}^h)_h$ a family of smooth measures on Markov chains $(\mathbb{P}_\varepsilon^h)_{h,\varepsilon}$, depending on a small regularisation parameter $\varepsilon$, such that

i) $\mathcal{I}^h(\mathbb{P}_\varepsilon^h)$ is close to $\mathcal{I}^h(\mathbb{P}^h)$ depending on $\varepsilon$,

ii) for all $\varepsilon > 0$, $\mathbb{P}_\varepsilon^h$ converges as $h \to 0$ (in a sense to be defined) to $\mathbb{P}_\varepsilon$, a smooth measure on $\Omega$,

iii) $\mathcal{I}^0(\mathbb{P}_\varepsilon)$ is close to $\mathcal{I}^0(\mathbb{P}^0)$ depending on $\varepsilon$.

**Remark 3.** We point out that this result is not a classic $\Gamma$-convergence result ($\mathcal{I}^h \to \mathcal{I}^0$) as we deal with discrete time Markov chains on the $\mathcal{I}^h$ side and continuous diffusion processes as arguments of $\mathcal{I}^0$. Our interest is in the behavior of numerical solutions to ($V^h$).

**Remark 4.** The discrete approximation of ($V^0$) corresponds to a relative entropic regularised problem and will be solved in Sections 4 and 4.3 using a Sinkhorn algorithm. Note that under stronger convexity assumptions on $F$, we can apply this strategy to any stochastic control problem of the type (3) simply by writing $F = F - S^\mathcal{I} + S^\mathcal{I}$ applying (6) only to the last term.

## 3.2 Convergence of the time discretisation

The proof shares similarities with $\Gamma$-convergence and relies on well-chosen regularization of the minimizers. Since it is long, we highlight the main steps and associated intermediate lemmas. Details are deferred to appendix 5.

We start with three technical lemmas. The first lemma concerns the scaling of the drift and diffusion coefficients, which is needed to satisfy the hard constraints in (H2) when using regularisations of the minimizers.

**Lemma 5 (Diffusion coefficients rescaling).** Let $X_t$ have law $\mathbb{P} \in \mathcal{P}^1$ with characteristics $b$, a satisfying $\lambda \mathrm{Id} < a < \Lambda \mathrm{Id}$ and $|b| < B$. For all $1 \gg \delta > 0$, let

$$ X_{\delta} = \sqrt{\alpha X} + \sqrt{\varepsilon B(t)}, $$

$\varepsilon = \delta \frac{\lambda+\Lambda}{\Lambda-\lambda}$ and $\alpha = 1 - \frac{2\varepsilon}{\lambda+\Lambda}$ depending on $\delta$, and $B(t)$ an independent standard Brownian motion. Then, $X_\delta$ has law $\mathbb{P}_\delta \in \mathcal{P}^1$ with characteristics $b_\delta, a_\delta$, and we have that:
i) Letting $T_{\sqrt{\alpha}}(x) = \sqrt{\alpha}x$, and $\mathbb{P}_{t,\delta} := (e_t)_{\#}\mathbb{P}_{\delta}$,

$$ \mathbb{P}_{t,\delta} := \mathcal{N}(0, \varepsilon t) \star [T_{\sqrt{\alpha}}]_{\#} \mathbb{P}_t, $$

ii) Scaling of the coefficients: $(\lambda + \delta) \mathrm{Id} < a_{\delta} < (\Lambda - \delta) \mathrm{Id}$, and $|b_{\delta}| \le (1 - \delta)B$.

iii) $\mathbb{P}_{\delta}$ narrowly converges to $\mathbb{P}$ as $\delta \searrow 0$ and $\mathcal{I}^0(\mathbb{P}_{\delta}) \le \mathcal{I}^0(\mathbb{P}) + O(\delta)$, where the notation $O(\delta)$ hides a constant that depends on the Lipschitz constant of $G$ and $\lambda, \Lambda, B$.

Moreover, the same holds for a discrete process $X^h$ with law $\mathbb{P}^h$, considering $\mathbb{P}_{\delta}^h$ the probability of $X_{\delta}^h(t) = \sqrt{\alpha}X^h(t) + \sqrt{\varepsilon}B(t)$ as above. We have, for $\alpha$ sufficiently close to 1,

iv) Scaling of the coefficients: $(\lambda + \delta) \mathrm{Id} < a_{\delta}^h < (\Lambda - \delta) \mathrm{Id}$, and $|b_{\delta}^h| \le (1 - \delta)B$.

v) $\mathbb{P}_{\delta}^h$ narrowly converges to $\mathbb{P}^h$ as $\delta \searrow 0$ and

$$ \mathcal{I}^h(\mathbb{P}_{\delta}^h) \le \mathcal{I}^h(\mathbb{P}^h) + O(\delta). \qquad (13) $$

*Proof.* Appendix 5.5. $\square$

The next lemma concerns the regularization technique for the time-continuous formula-
tion.

**Lemma 6** (Time (and space) regularization - the continuous case). Let $(\mathbb{P}^0, b^0, a^0)$ representing a Markov process $X(t)$ solving, for a Brownian motion $B(t)$,

$$ dX(t) = b^0(t, X)dt + a^0(t, X)dB(t) \qquad (14) $$

such that $\mathcal{I}(\mathbb{P}^0) < +\infty$. Fix $\sigma > 0$ and let $\varepsilon > 0$ such that $2\varepsilon\Lambda < \sigma$. We further assume that $|b^0| \le M/(1 + 2\varepsilon)$ and $a^0 \le \frac{\Lambda}{(1+2\varepsilon)^2}$ Id. Then, there exists another Markov process $X_{\varepsilon,\sigma}(t)$, solving the SDE as above for the quantities $b_{\varepsilon,\sigma}(t,x)$, $a_{\varepsilon,\sigma}(t,x)$ which are smooth in $t, x$ such that $\mathcal{I}(\mathbb{P}_{\varepsilon,\sigma}) \le \mathcal{I}(\mathbb{P}^0) + O(\varepsilon) + O(\sigma)$.

*Proof.* Appendix 5.6. $\square$

**Main steps of the proof:** Up to step 4, we prove that the value of the continuous problem is lower bounded by the lim sup of the values of the discrete problems. For that, we consider a probability with finite value for the continuous problem, rescale and regularize it, and then discretize it to obtain a candidate for the discrete problem.

* Step 0: Recall that $(\mathcal{V}^h)$ is well posed and $\mathbb{P}^h$ is markovian (see Theorem 3). We also have well-posedness and existence of a solution $\mathbb{P}^0$ of the time continuous problem $(\mathcal{V}^0)$ is known (see Theorem 1).

* Step 1: Given $\mathbb{P}^0 \in \mathcal{P}^1$ a minimizer of $(\mathcal{V}^0)$ we construct a regularised version $\mathbb{P}_{\varepsilon,\sigma}^0$ by applying lemma 6 ($\sigma$ depends on $\delta$). We get $\mathcal{I}^0(\mathbb{P}_{\varepsilon,\sigma}) \le \mathcal{I}^0(\mathbb{P}) + O(\sigma) + O(\varepsilon)$.

* Step 2: Set $\mathbb{P}_{\varepsilon,\sigma}^h = (e_{t_0,t_1,\dots,t_N})_{\#}\mathbb{P}_{\varepsilon,\sigma}$ the time discretisation of $\mathbb{P}_{\varepsilon,\sigma}$, the output of step 1. We want to compare $\mathcal{I}^h(\mathbb{P}_{\varepsilon,\sigma}^h)$ and $\mathcal{I}^0(\mathbb{P}_{\varepsilon,\sigma})$.

**Lemma 7.** For all $h$ sufficiently small w.r.t. $\delta, \varepsilon$, we have

i) The family $(\mathbb{P}_{\delta,\varepsilon}^h)_h$ satisfies the hypothesis of proposition 2 i) and in particular (see (36)) we have:

$$ h \mathcal{H}(\mathbb{P}_{\varepsilon,\sigma}^h | \bar{\mathbb{P}}^h) = \frac{1}{2} \sum_{i=1}^{N} h \mathbb{E}_{\mathbb{P}_{\delta,\varepsilon}} [S^{\mathcal{I}}(a^{\mathbb{P}_{\varepsilon,\sigma}}(ih, X_{ih})|\bar{a}] + O(h^{1/4}). $$
ii) It holds $\mathcal{I}^h(\mathbb{P}_{\varepsilon,\sigma}^h) \le \mathcal{I}^0(\mathbb{P}_{\varepsilon,\sigma}) + O(h^{\frac{1}{4}})$.

Proof. Appendix 5.7 □

* Step 3: For any minimizer of $(\mathcal{V}^h)$, denoted $\mathbb{P}^h$, and gathering the results in steps 1 and 2 we get

$$ \mathcal{I}^h(\mathbb{P}^h) \le \mathcal{I}^h(\mathbb{P}_{\varepsilon,\sigma}^h) \le \mathcal{I}^0(\mathbb{P}^0) + O(\sigma) + O(\varepsilon) + O(h^{\frac{1}{4}}). $$

We pass to the limit first in $h$ and then in $\varepsilon, \sigma$ to obtain $\limsup_{h \searrow 0} \mathcal{I}^h(\mathbb{P}^h) \le \mathcal{I}^0(\mathbb{P}^0)$. The next steps show that the lim inf of the values of the discrete problems is lower bounded by the value of the continuous problem.

* Step 4: We need the following definitions:

For $\mathbb{P}^h \in \mathcal{P}(\bigotimes_{i=0}^N \mathbb{R}_{t_i})$ Markovian, we define piecewise constant in time interpolants

$$ \tilde{\rho}_t^h = \sum_{i=0}^{N} \mathbf{1}_{[t_i, t_{i+1})} \mathbb{P}_i^h, \quad \tilde{b}_t^h = \sum_{i=0}^{N} \mathbf{1}_{[t_i, t_{i+1})} b_i^h, \quad \tilde{a}_t^h = \sum_{i=0}^{N} \mathbf{1}_{[t_i, t_{i+1})} a_i^h \quad (15) $$

$(b_i^h, a_i^h)$ defined in (41). We use the moment notations

$$ \tilde{m}_t^h = \tilde{\rho}_t^h \tilde{a}_t^h, \quad \tilde{n}_t^h = \tilde{\rho}_t^h \tilde{b}_t^h. \quad (16) $$

We define the proxy functional for $\mathcal{I}^0$

$$ \mathcal{J}(\rho, m, n) = \int_{t,x} F(m/\rho, n/\rho) + S^{\mathcal{I}}(n/\rho) d\rho_t(x)dt. \quad (17) $$

In particular we have for $\mathbb{P} \in \mathcal{P}^1$ of finite $\mathcal{I}^0$ energy: $\mathcal{I}^0(\mathbb{P}) = \mathcal{J}(\mathbb{P}_t, \mathbb{P}_t b^\mathbb{P}, \mathbb{P}_t a^\mathbb{P})$, and using Proposition 2 ii) : $\mathcal{J}(\tilde{\rho}^h, \tilde{m}^h, \tilde{n}^h) \le \mathcal{I}^h(\mathbb{P}^h)$ We need the following lemma:

**Lemma 8** (Time (and space) regularization - the discrete case). Let $(\mathbb{P}^h)_h$ be sequence of Markov Chains such that $\sup_h \mathcal{I}^h(\mathbb{P}^h) < +\infty$ and such that the fourth moment of $\mathbb{P}_0^h$ is uniformly bounded in $h$. Let $\delta, \varepsilon > 0$. Then, there exist $\mathbb{P}_{\delta, \varepsilon}^h \in \mathcal{P}(\mathbb{R}^{dN})$ and $\mathbb{P}_{\delta, \varepsilon} \in \mathcal{P}^1$ such that, up to a subsequence in $h$,

i) we have

$$ \mathcal{I}^h(\mathbb{P}_{\delta, \varepsilon}^h) \le \mathcal{I}^h(\mathbb{P}^h) + O(\varepsilon) + O(h) + O(\delta) \quad (18) $$

and

$$ \mathcal{I}^h(\mathbb{P}_{\delta, \varepsilon}^h) \ge \mathcal{J}(\tilde{\rho}_{\delta, \varepsilon}^h, \tilde{m}_{\delta, \varepsilon}^h, \tilde{n}_{\delta, \varepsilon}^h). \quad (19) $$

where $(\tilde{\rho}_{\delta, \varepsilon}^h, \tilde{m}_{\delta, \varepsilon}^h, \tilde{n}_{\delta, \varepsilon}^h)$ is the change of variable (15-16) associated to $\mathbb{P}_{\delta, \varepsilon}^h$

ii) $\tilde{\mathbb{P}}_{\delta, \varepsilon}^h$ narrowly converges to $\mathbb{P}_{\delta, \varepsilon} \in \mathcal{P}^1$ with characteristics $a_t^{\mathbb{P}_{\delta, \varepsilon}}, b_t^{\mathbb{P}_{\delta, \varepsilon}}$.

iii) In addition, we have

$$ \mathcal{I}^0(\mathbb{P}_{\delta, \varepsilon}) = \mathcal{J}(\mathbb{P}_{\delta, \varepsilon, t}, \tilde{\mathbb{P}}_{\delta, \varepsilon, t} b_t^{\tilde{\mathbb{P}}_{\delta, \varepsilon}}, \tilde{\mathbb{P}}_{\delta, \varepsilon, t} a_t^{\tilde{\mathbb{P}}_{\delta, \varepsilon}}) \le \liminf_{h \searrow 0} \mathcal{J}(\tilde{\rho}_{\delta, \varepsilon}^h, \tilde{m}_{\delta, \varepsilon}^h, \tilde{n}_{\delta, \varepsilon}^h). $$
Proof. Appendix 5.8. Note that it is possible to apply this lemma due to Hypothesis (H1) which ensures that the initial density’s fourth moment is finite. □

* **Step 5:** Lemma 8 i) and ii) gives $\mathcal{J}(\tilde{\rho}_{\delta,\varepsilon}^h, \tilde{m}_{\delta,\varepsilon}^h, \tilde{n}_{\delta,\varepsilon}^h) \le \mathcal{I}^h(\mathbb{P}^h) + O(\varepsilon) + O(\delta) + O(h)$ and passing to the limit in $h$ using Lemma 8 iii) and for a minimizer $\mathbb{P}^0$ of $(\mathcal{V}^0)$ we have

$$ \mathcal{I}^0(\mathbb{P}^0) \le \mathcal{I}^0(\mathbb{P}_{\delta,\varepsilon}) \le \liminf_{h\searrow 0} \mathcal{I}^h(\mathbb{P}^h) + O(\varepsilon) + O(\delta). \quad (20) $$

* **Last Step:** Applying **Step 5** to a sequence of minimizers $(\mathbb{P}^h)_h$ of $(\mathcal{V}^h)$ and using **Step 3**, we get $\mathcal{I}^0(\mathbb{P}^0) \le \lim_{\delta,\varepsilon\searrow 0} \mathcal{I}^0(\mathbb{P}_{\delta,\varepsilon}) \le \liminf_{h\searrow 0} \mathcal{I}^h(\mathbb{P}^h) \le \limsup_{h\searrow 0} \mathcal{I}^h(\mathbb{P}^h) \le \mathcal{I}^0(\mathbb{P}^0)$ which concludes the proof of Theorem 4 i).

# 4 Dual formulation and Sinkhorn algorithm

In this section the support of $\mathbb{P}^h$ is a finite product of compact subsets of $\mathbb{R}^d$: $\Omega^h := \bigotimes_{i=0}^N \chi_i$. This restriction is consistent with the space discretisation and truncation of our implementation (Section 4.3). In order to simplify the presentation and the notations we restrict to dimension $d=1$ and consider a dependence of the payoff $F$ just on the drift $b$, generalisations to dependence in $a$ follow the same lines. We also drop all additional constraints on the coefficients and the Kurtosis constraints since it does not seems necessary to obtain convergence of the Sinkhorn algorithm in the next section and also to let $h$ go to 0 together with the space discretisation step. We present the duality without rigorous proofs that are left for further studies. Let us just mentioned that Multi-Marginal OT on a product of compact space is usually easier to deal with as in [28] or [29] for instance.

The simplified discrete primal problem becomes:

$$ \inf_{\mathbb{P}^h \in \mathcal{P}(\Omega^h)} \mathbb{E}_{\mathbb{P}^h} \left( h \sum_{i=0}^{N-1} F(b_i^h(X_i^h)) \right) + \mathcal{D}(\mathbb{P}_0^h, \rho_0) + \mathcal{D}(\mathbb{P}_1^h, \rho_1) + h \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) \quad (21) $$

where $b_i^h(x_i) := \frac{1}{h} E_{\mathbb{P}_{i \to i+1}^h} (X_{i+1}^h - x_i)$ is given as a function of $\mathbb{P}^h$ as in (9). We also recall that the reference measure is defined by $\bar{\mathbb{P}}_{i \to i+1}^h = \mathcal{N}(x_i + h b_i^h(x_i), h \bar{a})$ and $\bar{\mathbb{P}}_0^h = \rho_0$.

## 4.1 Fenchel-Rockafellar duality

We remark that the drifts $b_i^h$ can be written using local linear functions of $\mathbb{P}^h$ (abusing again notations and using $\mathbb{P}^h$ for the distribution and their densities):

$$ b_i^h(x_i) = \frac{1}{h} \frac{\mathbb{E}_{\mathbb{P}_{i, i+1}^h}(x_{i, \cdot})(X_{i+1}^h - x_i)}{\mathbb{P}_i^h(x_i)} $$

where $\mathbb{P}_{i, i+1}^h(x_{i, \cdot}) := \mathbb{P}_i^h \mathbb{P}_{i \to i+1}^h$ is to be understood as the measure on $\mathbb{R}_{t_{i+1}}$ obtained by freezing the first variable $x_i$ in the joint $\mathbb{R}_{t_i} \times \mathbb{R}_{t_{i+1}}$ probability $\mathbb{P}_{i, i+1}^h$.

The primal problem (21) can be rewritten

$$ \inf_{\{\mathbb{P}^h \in \mathcal{P}(\Omega^h)\}} \mathcal{F}(\Delta^\dagger \mathbb{P}^h) + h \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) $$
using the linear change of variable:

$$ \Delta^{\dagger} \mathbb{P}^h := \{ (\mathbb{P}_i^h)_{i=0}^N, (x_i \to \frac{1}{h} \mathbb{E}_{\mathbb{P}_{i,i+1}^h}(x_i, \dots, (X_{i+1}^h - x_i))_{i=0}^{N-1} \}, $$

and modified functional:

$$ \mathcal{F}((m_{0,i})_{i=0}^N, (m_{1,i})_{i=0}^{N-1}) := \mathcal{D}(m_{0,0}, \rho_0) + \mathcal{D}(m_{0,N}, \rho_1) + h \sum_{i=0}^{N-1} \mathbb{E}_{m_{0,i}}(F(\frac{m_{1,i}}{m_{0,i}})). $$

In this section (again for simplicity) we use hard constraints for the loss function on initial/final marginals. It corresponds to the characteristic function $\mathcal{D}(\mu, \rho) := 0$ if $\mu = \rho$ and $+\infty$ otherwise. The (jointly convex) perspective functions $(m_0, m_1) \to h \mathbb{E}_{m_0}(F(m_1/m_0))$ are naturally extended to $+\infty$ if $m_0 \le 0$ and $0$ if $m_1 = m_0 = 0$. The $\dagger$ notation denotes that $\Delta^\dagger$ is the adjoint of the linear operator (23).

Fenchel Rockafellar duality yields the equivalent dual problem:

$$ \sup_{\Phi^h := (\phi_{0,i}^h)_{i=0}^N, (\phi_{1,i}^h)_{i=0}^{N-1}} -\mathcal{F}^*(-\Phi) - h \mathbb{E}_{\bar{\mathbb{P}}^h} \left( \exp\left(\frac{\Delta \Phi^h}{h}\right) - 1 \right). \quad (22) $$

The dual variable is a vector of continuous functions on the spaces $(\mathcal{X}_i)_{i=0}^N \Phi^h = (\phi_{0,i}^h)_{i=0}^N, (\phi_{1,i}^h)_{i=0}^{N-1}$. The second term in (22) is the Legendre Fenchel transform of the relative entropy and

$$ \Delta \Phi = \bigoplus_{i=0}^{N-1} \phi_{0,i} + \frac{1}{h}(x_{i+1} - x_i) \phi_{1,i} + \phi_{0,N}. \quad (23) $$

The Legendre Fenchel transform $\mathcal{F}^*$ of $\mathcal{F}$ is explicitly given, using the separability in space and in “times” $i$ by:

$$ \mathcal{F}^*(\Phi) = \mathbb{E}_{\rho_0}(\phi_{0,0} + h F^*(\frac{\phi_{1,0}}{h})) + \sum_{i=1}^{N-1} \chi_0(\phi_{0,i} + h F^*(\frac{\phi_{1,i}}{h})) + \mathbb{E}_{\rho_1}(\phi_{0,N}) $$

where $\chi_0(f) := 0$ if $f$ is the null function in $C(\mathcal{X})$ and $+\infty$ else. It is convenient to eliminate the $(\phi_{0,i})_{i=1}^{N-1}$ variables from $\Phi$ (we keep the same notations) and simplify problem (22) using

$$ \Delta \Phi = \phi_{0,0} + \frac{1}{h}(x_1 - x_0) \phi_{1,0} + \bigoplus_{i=1}^{N-1} h F^*(-\frac{1}{h}\phi_{1,i}) + \frac{1}{h}(x_{i+1} - x_i) \phi_{1,i} + \phi_{0,N} \quad (24) $$

and

$$ \mathcal{F}^*(\Phi) = \mathbb{E}_{\rho_0}(\phi_{0,0} + h F^*(\frac{\phi_{1,0}}{h})) + \mathbb{E}_{\rho_1}(\phi_{0,N}). \quad (25) $$

Applying classical Fenchel-Rockafellar duality ([30] theorem 6.3)

**Proposition 9** (Solutions of (22)). i) Problem (22) has a unique Markovian solution which can be expressed, in terms of $\Phi^h = (\phi_{0,i=0}^h, \phi_{0,i=N}^h), (\phi_{1,i}^h)_{i=0}^{N-1}$ a vector of continuous functions on the spaces $(\mathbb{R}_{t_i})$, the solution of the dual problem (22)-(24)-(25):

$$ \mathbb{P}^h = \exp \frac{1}{h} \Delta \Phi^h \bar{\mathbb{P}}, \\ \Delta \Phi = \phi_{0,0} + \frac{1}{h}(x_1 - x_0) \phi_{1,0} + \bigoplus_{i=1}^{N-1} h F^*(-\frac{1}{h}\phi_{1,i}) + \frac{1}{h}(x_{i+1} - x_i) \phi_{1,i} + \phi_{0,N}. $$
ii) $\mathbb{P}^h$ satisfies, for all $i$, the factorisation $\mathbb{P}^h$ can be tensorized with densities:

$$ \mathbb{P}^h = \rho_0 \prod_{i=0}^{N-1} \mathbb{P}_{i \to i+1}^h \qquad \mathbb{P}_{i \to i+1}^h = \frac{\mathbb{P}_{i,i+1}^h}{\mathbb{P}_i^h} $$

and $\forall i$

$$ \mathbb{P}_{i,i+1}^h = \mathcal{U}_i \exp \left( \frac{1}{h^2} (x_{i+1} - x_i) \phi_{1,i} \right) \bar{\mathbb{P}}_{i,i+1}^h \mathcal{D}_{i+1}. \qquad (26) $$

Each $(\mathcal{U}_i)$ and $(\mathcal{D}_i)$ are functions defined recursively Upward and Downward using $\Phi^h$ and $\bar{\mathbb{P}}^h$ by

$$ \left\{ \begin{array}{l} \mathcal{U}_0 = \exp\left(\frac{1}{h}\phi_{0,0}\right), \\ \mathcal{U}_{i+1} = \exp\left(F^*(-\frac{1}{h}\phi_{1,i+1})\right) \int \mathcal{U}_i \exp\left(\frac{1}{h^2}(x_{i+1}-x_i)\phi_{1,i}\right) \bar{\mathbb{P}}_{i,i+1}^h dx_i, \quad (27) \\ \text{for } i=0, N-2. \end{array} \right. $$

$$ \left\{ \begin{array}{l} \mathcal{D}_N = \exp\left(\frac{1}{h}\phi_{0,N}\right), \\ \mathcal{D}_i = \exp\left(F^*(-\frac{1}{h}\phi_{1,i})\right) \int \bar{\mathbb{P}}_{i,i+1}^h \exp\left(\frac{1}{h^2}(x_{i+1}-x_i)\phi_{1,i}\right) \mathcal{D}_{i+1} dx_{i+1}, \quad (28) \\ \text{for } i = N-1, 1. \end{array} \right. $$

Please recall that all variable subscripted by $i$ are to be understood as function of $x_i$.

**Remark 5.** The modifications to (24) and (25) in the more general case of time and diffusion dependent functions $(t_i, b_i, a_i) \to F_i(b_i, a_i)$ are simply $\Phi^h = (\phi_{0,i=0}^h, \phi_{0,i=N}^h), (\phi_{1,i}^h)_{i=0}^{N-1}, (\phi_{2,i}^h)_{i=0}^{N-1}$ the additional potentials $\phi_{2,i}^h$ are dual to $\mathbb{P}_i^h a_i^h$ and

$$ \begin{aligned} \Delta \Phi = & \phi_{0,0} + \frac{1}{h}(x_1 - x_0)\phi_{1,0} + \frac{1}{h}(x_1 - x_0)^2\phi_{2,0} + \\ & \bigoplus_{i=1}^{N-1} h F_i^*(-\frac{1}{h}\phi_{1,i}, -\frac{1}{h}\phi_{2,i}) + \frac{1}{h}(x_{i+1} - x_i)\phi_{1,i} + \frac{1}{h}(x_{i+1} - x_i)^2\phi_{2,i} \\ & + \phi_{0,N} \end{aligned} $$

and

$$ \mathcal{F}^*(\Phi) = \mathcal{D}_{\rho_0}^*(\phi_{0,0} + h F_0^*\left(\frac{\phi_{1,0}}{h}, \frac{\phi_{2,0}}{h}\right)) + \mathcal{D}_{\rho_1}^*(\phi_{0,N}) $$

where $\mathcal{D}_{\rho_i}^*$ is the Legendre Fenchel Transform of the (convex) marginal fidelity term $\rho \to \mathcal{D}(\rho, \rho_i)$.

## 4.2 Multi-Marginal Sinkhorn Algorithm

We are now working with a fixed $h$ and drop the dependence in the notations for clarity.
The simplest interpretation of Sinkhorn Algorithm (see [1]) is to perform an iterative coordinate-wise (in the components of $\Phi$) ascent to maximize the concave dual problem (22-24-25):

$$ \Phi := (\phi_{0,0}, (\phi_{1,i})_{i=0}^{N-1}, \phi_{0,N}) \qquad \sup_{\Phi} \quad -\mathcal{F}^*(-\Phi) - h \mathbb{E}_{\bar{\mathbb{P}}} (\exp\left(\frac{1}{h}\Delta\Phi\right) - 1). $$

One cycle of dual variable optimisation will be indexed by $k$ and the potentials updated “à la Gauss Seidel” in the inner loop over the $N+2$ functions (Kantorovich potentials) in $\Phi$ (the order is not important). Notations are a little bit more involved than in the classical two marginals problem, to ease the presentation we set:

$$ \left\{ \begin{array}{l} \Phi^{k,0} = (\phi_{0,0}, (\phi_{1,i}^k)_{i=0}^{N-1}, \phi_{0,N}^k), \\ \Phi^{k,i+1} = (\phi_{0,0}^{k+1}, (\phi_{1,0}^{k+1} \cdots \phi_{1,i-1}^{k+1}, \phi_{1,i}, \phi_{1,i+1}^k \cdots \phi_{1,N}^k), \phi_{0,N}^k), \quad i=0 \cdots N-1, \quad (29) \\ \Phi^{k,N+2} = (\phi_{0,0}^{k+1}, (\phi_{1,i}^{k+1})_{i=0}^{N-1}, \phi_{0,N}^k). \end{array} \right. $$

One Sinkhorn cycle $k \to k+1$, the updates from $\Phi^{k,0}$ to $\Phi^{k+1,0}$ can be written in compact form as a loop on its components:

$$ \left\{ \begin{aligned} \phi_{0,0}^{k+1} &= \arg\sup_{\phi_{0,0}} -\mathcal{F}^*(-\Phi^{k,0}) - h \mathbb{E}_{\bar{\mathbb{P}}} \left( \exp\left(\frac{1}{h}\Delta\Phi^{k,0}\right) - 1 \right), \\ \phi_{1,j+1}^{k+1} &= \arg\sup_{\phi_{1,j}} -\mathcal{F}^*(-\Phi^{k,j+1}) - h \mathbb{E}_{\bar{\mathbb{P}}} \left( \exp\left(\frac{1}{h}\Delta\Phi^{k,j+1}\right) - 1 \right), \quad i=0 \cdots N-1, \\ \phi_{0,N}^{k+1} &= \arg\sup_{\phi_{0,N}} -\mathcal{F}^*(-\Phi^{k,N+1}) - h \mathbb{E}_{\bar{\mathbb{P}}} \left( \exp\left(\frac{1}{h}\Delta\Phi^{k,N+2}\right) - 1 \right). \end{aligned} \right. \quad (30) $$

Each of these maximization problems is strictly concave and sufficiently smooth (depending on $F^*$). They are also separable in space and (30) amounts to solve in sequence the following set of equations ($\partial F^*$ is the Frechet derivative, a gradient in finite dimension):

$$ \left\{ \begin{array}{l} \rho_0 = \exp\left(\frac{1}{h}\phi_{0,0}^{k+1}(x_0)\right) \int \exp\left(\frac{1}{h^2}(x_1 - x_0)\phi_{1,0}^k\right) \bar{P}_{0,1}^h D_1^k dx_1 \\ \rho_0 \partial F^*(-\frac{1}{h}\phi_{1,0}^{k+1}) = \exp\left(\frac{1}{h}\phi_{0,0}^{k+1}\right) \int \frac{1}{h}(x_1 - x_0) \exp\left(\frac{1}{h^2}(x_1 - x_0)\phi_{1,0}^{k+1}\right) \bar{P}_{0,1}^h D_1^k dx_1 \\ 0 = U_i^{k+1} \int \left( \partial F^*(-\frac{1}{h}\phi_{1,i}^{k+1}) + \frac{1}{h}(x_{i+1} - x_i) \right) \exp\left( F^*(-\frac{1}{h}\phi_{1,i}^{k+1}) + \frac{1}{h^2}(x_{i+1} - x_i) \phi_{1,0}^{k+1} \right) \\ \qquad \cdots \bar{P}_{i,i+1}^h D_{i+1}^k dx_{i+1} \quad \text{for all } i=1, N-1 \\ \rho_1 = U_{N-1}^{k+1} \exp\left(\frac{1}{h}\phi_{0,N}^{k+1}\right) \end{array} \right. $$

Please recall that all variable subscripted by $i$ are to be understood as function of $x_i$. Each line is to be understood as point-wise in space. The functions $(\mathcal{U}_i^k, \mathcal{D}_i^k)$ are defined by the recursions (27-28) applied to the iterative Gauss-Seidel update of the potentials: (29)

$$ \left\{ \begin{array}{l} \mathcal{U}_0^k = \exp\left(\frac{1}{h}\phi_{0,0}^k\right) \\ \mathcal{U}_{i+1}^k = \exp\left(F^*(-\frac{1}{h}\phi_{1,i+1}^k)\right) \int \mathcal{U}_i^k \exp\left(\frac{1}{h^2}(x_{i+1} - x_i)\phi_{1,i}^k\right) \bar{P}_{i,i+1}^h dx_i \quad (31) \\ \text{for } i=0, N-2 \end{array} \right. $$
$$ \left\{ \begin{array}{l} \mathcal{D}_N^k = \exp\left(\frac{1}{h}\phi_{0,N}^k\right) \\ \mathcal{D}_i^k = \exp\left(F^*(-\frac{1}{h}\phi_{1,i}^k)\right) \int \bar{\mathbb{P}}_{i,i+1}^h \exp\left(\frac{1}{h^2}(x_{i+1} - x_i)\phi_{1,i}^k\right) \mathcal{D}_{i+1}^k dx_{i+1} \quad (32) \\ \text{for } i = N-1, 1 \end{array} \right. $$

The numerical study in Section 4.3 gives experimental convergence curves. In the case of finite dimension in space, i.e, $\mathbb{P}^h$ as support on a truncated grid in space. The analysis of the convergence (in $k$) of this algorithm should rely on a multi-marginal extension of [31], see also [29] or [28] and the references therein. Note however that they apply to “standard” multi-marginal constraints while we use a convex payoff of the moments of $\mathbb{P}^h$. This problem is closer in spirit to the “weak” OT [32] setting or rather the entropic regularization of a weak OT problem.

## 4.3 Numerical experiments

In all experiments, we use the function $F = \gamma \|\cdot\|^2$ with a large parameter $\gamma$ (typically $1.E06$), this is to enforce a soft martingale constraint (no drift). The space discretisation is chosen as $dx = (10\bar{\sigma})/K\sqrt{h}$ the parameter $K$ corresponds to the number of “active” points in the reference measure Gaussian kernel (typically $K = 64$). Note that, we extend the test cases to reference measures $\bar{\mathbb{P}}^h$ with the variable in time $\bar{\sigma}_i = \sqrt{a_i}$. The finest computed discretisation in time corresponds to $1/h = N_t = 513$ but we will also provide convergence curves in $N_t$.

The first test case, see figure 1, captures a diffusion between Gaussian source and target measures in convex order. The means are $\mu = 0.5$ and respective standard deviation $\sigma_0 = 0.05$ and $\sigma_1 = 0.1$ (Figure 1 (a)). The reference measure volatility (b) is constant and satisfies $0.009 = \bar{\sigma} = \sqrt{a} > \sqrt{\sigma_1^2 - \sigma_0^2} = 0.0075$. It is over-diffusive with respect to the target measure and not an admissible solution for the problem. Figure 1 (e) shows the method recovers the constant $\sigma = 0.0075$ solution. The volatility drops at the boundary because we have truncated the domain and strongly penalized the drift (f), the density is close to 0 there, see (c) for the curve in time of marginal densities. Figure 1 (h) shows that the scaled relative entropy and discrete specific entropy get closer as $N_T$ the number of timesteps increase. Figure 1 (g) confirms that under our choice of truncation and space discretisation the cost of the algorithm is of order $O(N_T^{3/2})$.

The second test case, see figure 2 is similar to figure 1 but the target measure is a mixture of Gaussians (a). $\rho_1(x) = p\mathcal{N}(x; 0.5, 0.1) + \frac{1-p}{2}\mathcal{N}(x; 0.5 - d_1, 0.05) + \frac{1-p}{2}\mathcal{N}(x; 0.5 + d_1, 0.05)$ ($p = 0.6$). The optimal volatility surface The results are consistent with those of the first test case, we have a nice convergence in $N_T$ for the entropies (h). The volatility surface (e) is not simple.

In the third test case, figure 3, the setting is similar to figure 2 except for the reference volatility which is discontinuous in time with a drop around time $0.5$: $\bar{\sigma} = 0.03$ if $|t-0.5| \le 0.2$ and $0.01$ else. Again results are consistent with the previous test cases. We can also observe the regularizing effect of the specific entropy on (e) (which is asymptotically close to a weighted $\mathcal{L}^2$ norm).

For the final test figure 4, the solution cannot be a martingale, the source and targets (a) are not in convex order. $\rho_0(x) = \mathcal{N}(x; 0.5, 0.1)$ if $x \le 0.5$ and $\mathcal{N}(x; 0.5, 0.3)$ if $x > 0.5$, $\rho_1(x) = \mathcal{N}(x; 0.5, 0.5)$ if $x \le 0.5$ and $\mathcal{N}(x; 0.5, 0.2)$ if $x > 0.5$. The optimal solutions (c) (e) (f) generate diffusion where possible ($x < 0.5$) and use the drift in the non-convex order side ($x > 0.5$).
<figure><img src="image_2.png" /><figcaption>(a) Source and target distributions</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) Marginal surface</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(e) Final volatility surface</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(g) Runtime curve</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Reference measure volatility surface</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Sinkhorn convergence of the relative $L^{\infty}$ error</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f) Final drift surface</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(h) Convergence of the entropy and specific entropy</figcaption></figure>

<caption>Fig. 1: Gaussian to gaussian with overdiffusive reference volatility</caption>
<figure><img src="image_10.png" /><figcaption>Fig. 2: A Gaussian to a mixture of gaussian</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) Source and target distributions</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Reference measure volatility surface</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Marginal surface</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Sinkhorn convergence of the relative $L^{\infty}$ error</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e) Final volatility surface</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) Final drift surface</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(g) Runtime curve</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(h) Convergence of the entropy and specific entropy</figcaption></figure>

<caption>**Fig. 3:** A Gaussian to sum of gaussians with a discontinuous in time volatility</caption>
<figure><img src="image_19.png" /><figcaption>(a) Source and target distributions</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) Marginal surface</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(e) Final volatility surface</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(g) Runtime curve</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Reference measure volatility surface</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(d) Convergence of the relative $L^\infty$ error</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(f) Final drift surface</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(h) Convergence of the entropy and specific entropy</figcaption></figure>

<caption>Fig. 4: Marginals not in convex order</caption>
References

[1] Peyré, G., Cuturi, M.: Computational optimal transport. Found. Trends Mach. Learning **11**, 355–607 (2019)

[2] Léonard, C.: A survey of the Schrödinger problem and some of its connections with optimal transport. Discrete Contin. Dyn. Syst. **34**, 1533–1574 (2014)

[3] Galichon, A.: Optimal Transport Methods in Economics, 1st edn. Princeton University Press. https://EconPapers.repec.org/RePEc:pup:pbooks:10870

[4] Benamou, J.-D.: Optimal transportation, modelling and numerical simulation. Acta Numerica **30**, 249–325 (2021) https://doi.org/10.1017/S0962492921000040

[5] Tan, X., Touzi, N.: Optimal transportation under controlled stochastic dynamics. The Annals of Probability **41**(5) (2013) https://doi.org/10.1214/12-aop797

[6] Gantert, N.: Some large deviations of Brownian motion. Universitat Bonn (1991)

[7] Benamou, J.-D., Brenier, Y.: A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem. Numer. Math. **84**, 375–393 (2000)

[8] Papadakis, N., Peyré, G., Oudet, E.: Optimal transport with proximal splitting. SIAM J. Imaging Sci. **7**, 212–238 (2014) https://doi.org/10.1137/130920058

[9] Guo, I., Loeper, G., Wang, S.: Calibration of local-stochastic volatility models by optimal transport. Mathematical Finance **32**(1), 46–77 (2022) https://doi.org/10.1111/mafi.12335 https://onlinelibrary.wiley.com/doi/pdf/10.1111/mafi.12335

[10] Benamou, J.-D., Carlier, G., Nenna, L.: Generalized incompressible flows, multi-marginal transport and Sinkhorn algorithm. Numer. Math. **142**, 33–54 (2019) https://doi.org/10.1007/s00211-018-0995-x

[11] Benamou, J.-D., Carlier, G., Di Marino, S., Nenna, L.: An entropy minimization approach to second-order variational mean-field games. Math. Models Methods Appl. Sci. **29**, 1553–1583 (2019)

[12] Benamou, J.-D., Chazareix, G., Loeper, G.: From entropic transport to martingale transport, and applications to model calibration (2024). https://arxiv.org/abs/2406.11537

[13] Bunne, C., Schiebinger, G., Krause, A., Regev, A., Cuturi, M.: Optimal transport for single-cell and spatial omics (2024). https://doi.org/10.1038/s43586-024-00334-2

[14] Marco Avellaneda, R.H. Craig Friedman, Samperi, D.: Calibrating volatility surfaces via relative-entropy minimization. Applied Mathematical Finance **4**(1), 37–64 (1997)
[15] Backhoff-Veraguas, J., Unterberger, C.: On the specific relative entropy between martingale diffusions on the line (2023)

[16] Backhoff-Veraguas, J., Zhang, X.: Specific Wasserstein divergence between continuous martingales (2024). https://arxiv.org/abs/2404.19672

[17] Peyré, G.: Entropic approximation of Wasserstein gradient flows. SIAM J. Imaging Sci. 8, 2323-2351 (2015)

[18] Guo, I., Loeper, G.: Path dependent optimal transport and model calibration on exotic derivatives. The Annals of Applied Probability 31(3), 1232-1263 (2021)

[19] Guo, I., Loeper, G., Oblój, J., Wang, S.: Joint modeling and calibration of spx and vix by optimal transport. SIAM Journal on Financial Mathematics 13(1), 1-31 (2022) https://doi.org/10.1137/20M1375905

[20] Guo, I., Loeper, G., Wang, S.: Calibration of local-stochastic volatility models by optimal transport. Mathematical Finance 32(1), 46-77 (2022) https://doi.org/10.1111/mafi.12335 https://onlinelibrary.wiley.com/doi/pdf/10.1111/mafi.12335

[21] Huesmann, M., Trevisan, D.: A Benamou-Brenier formulation of martingale optimal transport. Bernoulli 25(4A), 2729-2757 (2019) https://doi.org/10.3150/18-BEJ1069

[22] Backhoff-Veraguas, J., Beiglböck, M., Schachermayer, W., Tschiderer, B.: The structure of martingale Benamou-Brenier in $R^d$ (2023)

[23] Föllmer, H.: Doob Decomposition, Dirichlet Processes, and Entropies on Wiener Space. In: Festschrift in Honour of Masatoshi Fukushima's Beiju. Springer Proceedings in Mathematics & Statistics, (2022). https://hal.science/hal-03589442

[24] Henry-Labordere, P.: From (Martingale) Schrodinger bridges to a new class of Stochastic Volatility Models. arXiv (2019). https://doi.org/10.48550/ARXIV.1904.04554. https://arxiv.org/abs/1904.04554

[25] Benamou, J.-D., Chazareix, G., Loeper, G.: From entropic transport to martingale transport, and applications to model calibration (2024). https://arxiv.org/abs/2406.11537

[26] Stroock, D.W., Varadhan, S.R.S.: Multidimensional Diffusion Processes. Grundlehren der mathematischen Wissenschaften. Springer. https://books.google.de/books?id=DuDsmoyqCy4C

[27] Kushner, H.J.: On the Weak Convergence of Interpolated Markov Chains to a Diffusion. The Annals of Probability 2(1), 40-50 (1974) https://doi.org/10.1214/aop/1176996750

[28] Di Marino, S., Gerolin, A.: An optimal transport approach for the schrödinger
bridge problem and convergence of sinkhorn algorithm. Journal of Scientific Computing 85, 27 (2020) https://doi.org/10.1007/s10915-020-01325-7

[29] Carlier, G.: On the linear convergence of the multimarginal sinkhorn algorithm. SIAM Journal on Optimization 32(2), 786-794 (2022) https://doi.org/10.1137/21M1410634

[30] Carlier, G.: Classical and Modern Optimization, (2021)

[31] Cominetti, R., Martín, J.S.: Asymptotic analysis of the exponential penalty trajectory in linear programming. Math. Program. 67, 169-187 (1994)

[32] Backhoff-Veraguas, J., Pammer, G.: Applications of weak transport theory. Bernoulli 28(1), 370-394 (2022) https://doi.org/10.3150/21-BEJ1346

[33] Azencott, R.: Densité des diffusions en temps petit: développements asymptotiques. In: Séminaire de Probabilités XVIII 1982/83: Proceedings, pp. 402-498 (1984). Springer

[34] Aronson, D.G.: Bounds for the fundamental solution of a parabolic equation. Bulletin of the American Mathematical Society, 73, 890-903 (1967)

[35] Ambrosio, L., Gigli, N., Savaré, G.: Gradient Flows in Metric Spaces and in the Space of Probability Measures, 2. ed edn. Lectures in Mathematics ETH Zürich. Birkhäuser. OCLC: 254181287

[36] Ambrosio, L., Stra, F., Trevisan, D.: Weak and strong convergence of derivations and stability of flows with respect to MGH convergence (2016)

[37] Trevisan, D.: Well-posedness of multidimensional diffusion processes with weakly differentiable coefficients. Electronic Journal of Probability 21 (2015) https://doi.org/10.1214/16-EJP4453

# 5 Appendix

## 5.1 Proof of proposition 2 [Specific Relative Entropy]

**Proof of Prop. 2 i)** We have $X_t^h = X_{ih}$, where $(X_t)_{t \in [0,1]}$ solves $dX_t = b_t(X_t)dt + a_t(X_t)^{1/2}dB_t$.

**Step 1:** We recall the ingredient we need for small-time asymptotics of the heat kernel. Let $p_{a,b}^h(t,x,y)$ denote the transition density associated with the diffusion process $dX_t = b(t, X_t)dt + \sigma(t, x)dB_t$ with values in $\mathbb{R}^d$ for some $d$-dimensional standard Brownian motion, and $a(t, x) = \sigma(t, x)\sigma(t, x)^*$. We assume that both functions $(t, x) \mapsto b(t, x)$ and $(t, x) \mapsto a(t, x)$ are twice continuously differentiable, and that the matrix $a(t, x)$ is invertible for every $(t, x) \in [0, 1] \times \mathbb{R}^d$. For $x, y \in \mathbb{R}^d$ and $0 \le t \le t+h \le 1$, let

$$p_{a,b}(t, h, x, y)dy = \mathbb{P}(X_{t+h} \in dy | X_t = x)$$

denote the transition probability associated with the diffusion process with drift $b$ and diffusion matrix $a$. From [33], the following asymptotic expansion is known: for every compact
$U \subset \mathbb{R}^d$, there exists $\tau > 0$ such that for $\|x - y\| \le \tau$, we have

$$p_{a,b}(t, h, x, y) = \frac{1}{(2\pi h)^{d/2} \det(a(t, x)))^{1/2}} e^{\left(\frac{1}{2h}(y-x)^* a(t, x)^{-1}(y-x)\right)} (1 + hR(t, h, x, y)) \quad (33)$$

where $\sup_{0 \le t \le t+h \le 1, x, y \in U} |R(t, h, x, y)| \le C$ and $R$ only depends on $\partial_t^i \partial_x^j \sigma$ and $\partial_t^i \partial_x^j b$, for $0 \le i + j \le 2$, computed on $x$ and on $f_u$, $0 \le u \le 1$, where $f_u$ is the unit length geodesic joining $x$ and $y$ for the Riemannian metric associated with $a(t, x)^{-1}$.

**Step 2: We separate the basic Kullback-Leibler divergence into a local and a global part.** Define $q_a(t, h, x, y) = \frac{1}{(2\pi h)^{d/2} \det(a(t, x))^{1/2}} \exp(-\frac{1}{2h}(y-x)^* a(t, x)^{-1}(y-x))$. From (33), a first-order expansion in $h$ yields $\log(\frac{p_{a,b}(t,h,x,y)}{p_{\bar{a},0}(t,h,x,y)}) = \log(\frac{q_a(t,h,x,y)}{q_a(t,h,x,y)}) + h\tilde{R}(t,h,x,y)$ where $\tilde{R}$ has the same properties as $R$ with dependence upon $a, b$ and $\bar{a}$, i.e. it is valid for $\|x - y\| \le \tau$. For $\|x - y\| \ge \tau$, we always have the Aaronson's type estimate [34] that reads

$$c_{-}h^{-d/2}e^{-C_{-}h^{-1}\|x-y\|} \le p_{a,b}(t, h, x, y) \le c_{+}h^{-d/2}e^{-C_{+}h^{-1}\|x-y\|} \quad (34)$$

where $c_\pm, C_\pm$ depend on $\inf_{t,x} \|a(t,x)\| > 0$ and $\sup_{t,x} (\|b(t,x)\|, \|a(t,x)\|)$, thus requiring uniform boundedness. Thus

$$\log \left( \frac{p_{a,b}(t, h, x, y)}{p_{\bar{a},0}(t, h, x, y)} \right) = \left( \log \left( \frac{q_a(t, h, x, y)}{q_{\bar{a}}(t, h, x, y)} \right) + h \tilde{R}(t, h, x, y) \right) \mathbf{1}_{\|x-y\| \le \eta} + V_{\eta}(t, h, x, y),$$

where $V_\eta(t, h, x, y) = \log(\frac{p_{a,b}(t,h,x,y)}{p_{\bar{a},0}(t,h,x,y)}) \mathbf{1}_{\|x-y\| \ge \eta}$. The estimate (34) yields

$$\begin{align*}
\mathbb{E}_{\mathbb{P}}[V_{\eta}(ih, h, X_{ih}, X_{(i+1)h})]
&\le \frac{c_{a,b}^{+}}{c_{\bar{a},0}^{-}} \mathbb{P}([\|X_{(i+1)h} - X_{ih}\| \ge \eta] - (C_{a,b}^{+} - C_{a,0}^{-})h^{-1} \mathbb{E}_{\mathbb{P}}[\|X_{(i+1)h} - X_{ih}\|^2 \mathbf{1}_{\|X_{(i+1)h} - X_{ih}\| \ge \eta}]) \\
&\le \left(\frac{c_{a,b}^{+}}{c_{\bar{a},0}^{-}} - (C_{a,b}^{+} - C_{a,0}^{-})h^{-1} \mathbb{E}_{\mathbb{P}}[\|X_{(i+1)h} - X_{ih}\|^4]^{1/2}\right) \eta^{-1/2} \mathbb{E}[\|X_{(i+1)h} - X_{ih}\|]^{1/2} = O(h^{1/4}).
\end{align*}$$

Likewise, $\mathbb{E}_{\mathbb{P}}[h\tilde{R}(ih, h, X_{ih}, X_{(i+1)h}) \mathbf{1}_{\|X_{(i+1)h} - X_{ih}\| \le \eta}] \le Ch$, so that

$$\mathbb{E}_{\mathbb{P}}[\log\left(\frac{p_{a,b}(ih, h, X_{ih}, X_{(i+1)h})}{p_{\bar{a},0}(ih, h, X_{ih}, X_{(i+1)h})}\right)] = \mathbb{E}_{\mathbb{P}}[\log\left(\frac{q_{a,b}(ih, h, X_{ih}, X_{(i+1)h})}{q_{\bar{a},0}(ih, h, X_{ih}, X_{(i+1)h})}\right)] + O(h^{1/4}) \quad (35)$$

and therefore, according to (35),

$$\begin{align*}
\mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) &= \mathbb{E}_{\mathbb{P}} \left[ \log \left( \prod_{i=1}^{N} \frac{p_{a,b}(ih, h, X_{ih}, X_{(i+1)h})}{p_{\bar{a},0}(ih, h, X_{ih}, X_{(i+1)h})} \right) \right] \\
&= \sum_{i=1}^{N} \mathbb{E}_{\mathbb{P}} \left[ \log \left( \frac{q_a(ih, h, X_{ih}, X_{(i+1)h})}{q_{\bar{a}}(ih, h, X_{ih}, X_{(i+1)h})} \right) \right] + O(Nh^{1/4}).
\end{align*}$$

**Step 3. Control of the entropy of the proxy transitions.** We have

$$\log \left( \frac{q_a(ih, h, X_{ih}, X_{(i+1)h})}{q_{\bar{a}}(ih, h, X_{ih}, X_{(i+1)h})} \right) = -\frac{1}{2} \log \left( \frac{\det(a(ih, X_{ih}))}{\det(\bar{a}(ih, X_{ih}))} \right) \\ - \frac{1}{2} h^{-1} (X_{(i+1)h} - X_{ih})^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) (X_{(i+1)h} - X_{ih}).$$
Next, we have

$$
\begin{align*}
& \frac{1}{2} h^{-1} (X_{(i+1)h} - X_{ih})^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) (X_{(i+1)h} - X_{ih}) \\
&= -\frac{1}{2} h^{-1} \left( \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s \right)^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) \left( \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s \right) + \zeta_i^h,
\end{align*}
$$

where

$$
\begin{align*}
\zeta_i^h &= h^{-1} \left( \int_{ih}^{(i+1)h} b(s, X_s) ds \right)^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) \left( \int_{ih}^{(i+1)h} b(s, X_s) ds \right) \\
&\quad + h^{-1} \left( \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s \right)^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) \left( \int_{ih}^{(i+1)h} b(s, X_s) ds \right) \\
&\quad + h^{-1} \left( \int_{ih}^{(i+1)h} b(s, X_s) ds \right)^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) \left( \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s \right)
\end{align*}
$$

The boundedness of $b$ and $\sigma$ (and its inverse) readily yields $\max_{1\le i\le N} \mathbb{E}_{\mathbb{P}}[\|\zeta_i^h\|] = O(h^{1/2})$ by applying Cauchy-Schwarz inequality repeatedly. Also, writing $h^{-1/2} \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s = \sigma(ih, X_{ih}) + \xi_i^h$, where $\xi_i^h = h^{-1/2} \int_{ih}^{(i+1)h} (\sigma(s, X_s) - \sigma(ih, X_{ih})) dB_s$ satisfies $\max_{1\le i\le N} \mathbb{E}_{\mathbb{P}}[\|\xi_i^h\|] = O(h^{1/2})$ thanks to the smoothness of $(t, x) \mapsto \sigma(t, x)$ (at least Lipschitz in both variables) and standard moment estimates, we obtain

$$
\begin{align*}
& -\frac{1}{2}h^{-1} \left( \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s \right)^* (a(ih, X_{ih})^{-1} - \bar{a}(ih, X_{ih})^{-1}) \left( \int_{ih}^{(i+1)h} \sigma(s, X_s) dB_s \right) \\
&= \frac{1}{2} \operatorname{Tr} (\bar{a}(ih, X_{ih})^{-1} (a(ih, X_{ih}) - \bar{a}(ih, X_{ih}))) + \tilde{\xi}_i^h,
\end{align*}
$$

where now, in the same way as for the remainder term $\xi_i^h$, we readily have
$\max_{1\le i\le N} \mathbb{E}_{\mathbb{P}}[\|\tilde{\xi}_i^h\|] = O(h^{1/2})$. Putting together all our estimates, we obtain

$$
\log \left( \frac{q_a(ih, h, X_{ih}, X_{(i+1)h})}{q_{\bar{a}}(ih, h, X_{ih}, X_{(i+1)h})} \right) = -\frac{1}{2} \log \left( \frac{\det(a(ih, X_{ih}))}{\det(\bar{a}(ih, X_{ih}))} \right) \\
+ \frac{1}{2} \operatorname{Tr} (\bar{a}(ih, X_{ih})^{-1} (a(ih, X_{ih}) - \bar{a}(ih, X_{ih}))) + \tilde{\zeta}_i^h,
$$

for some stochastic remainder term $\tilde{\zeta}_i^h$ that satisfies the moment estimate
$\max_{1\le i\le N} \mathbb{E}_{\mathbb{P}}[\|\tilde{\zeta}_i^h\|] = O(h^{1/2})$. By definition of the relative entropy $S^I$, we thus obtain

$$
\mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) = \frac{1}{2} \sum_{i=1}^{N} \mathbb{E}_{\mathbb{P}}[S^I(a(ih, X_{ih}) | \bar{a}(ih, X_{ih}))] + O(Nh^{1/4}).
$$

Finally, using that $h \sim N^{-1}$, we obtain

$$
h \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) = \frac{1}{2} \sum_{i=1}^{N} h \mathbb{E}_{\mathbb{P}} [S^{\mathcal{I}}(a(ih, X_{ih}) | \bar{a}(ih, X_{ih}))] + O(h^{1/4}) \quad (36)
$$

and $\lim_{h \to 0} h \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) = \mathbb{E}_{\mathbb{P}} \left( \int_0^1 S^{\mathcal{I}}(a(s, X_s) | \bar{a}(s, X_s) ds) \right)$ by Riemann approximation,
which is the desired result. More specifically, one readily checks that

$$
\int_{ih}^{(i+1)h} S^{\mathcal{I}}(a(s, X_s) | \bar{a}(s, X_s) ds = h S^{\mathcal{I}}(a(ih, X_{ih}) | \bar{a}(ih, X_{ih}) + \rho_i^h,
$$

where $\max_{1\le i\le N} \mathbb{E}[\|\rho_i^h\|] = O(h^{1/2})$ using that both $S^{\mathcal{I}}$ and $a, \bar{a}$ are smooth functions.
**Proof of Proposition 2 ii)** Since $\mathbb{P}^h$ is Markov we rewrite the relative Entropy using the transition probabilities $\mathcal{H}(\mathbb{P}^h|\bar{\mathbb{P}}^h) = \mathcal{H}(\mathbb{P}_0|\rho_0) + \sum_{i=0}^{N-1} \mathbb{E}_{\bar{\mathbb{P}}_i^h}(\mathcal{H}(\mathbb{P}_{i \to i+1}^h|\bar{\mathbb{P}}_{i \to i+1}^h))$. It is therefore sufficient to focus on one transition and show $\mathcal{H}(\mathbb{P}_{i \to i+1}^h|\bar{\mathbb{P}}_{i \to i+1}^h) \ge \frac{1}{2} S^\mathcal{I}(a_i^h(X_i^h)|\bar{a})$. Let us simplify the notations and look at $\mathcal{H}(\mathbb{P}_x|\bar{\mathbb{P}}_x)$ where $\mathbb{P}_x$ is a probability measure over a set $Y = \mathbb{R}_{t_{i+1}}$ parameterized by $x \in X = \mathbb{R}_t$ and $\bar{\mathbb{P}}_x$ a probability measure over $Y$ with normal law $\Gamma_{x \to y}(0, \bar{a})$. This inequality is a consequence of the well-known fact that the probability that minimizes the entropy under constrained second moments is a Gaussian measure. In the following, we prove it using the dual formulation of entropy, recalled hereafter, see [35, lemma 9.4.4.]

$$ \left\{ \begin{array}{l} \mathcal{H}(\mathbb{P}_x|\bar{\mathbb{P}}_x) = \sup_{f \in C_b^0(Y)} \beta(f) \\ \beta(f) := \int f(y) d\mathbb{P}_x(y) - \int \exp(f(y)) - 1 d\bar{\mathbb{P}}_x(y). \end{array} \right. \qquad (37) $$

We plan on restricting the test functions $f$ to quadratic functions $y \to p_{x,C,A}(y) = C + A(y-x)^2$ which are not bounded functions. To do so, we remark that the truncation

$$ p_{x,C,A,R}(y) = \begin{cases} p_{x,C,A}(y) & \text{if } |y-x| \le R \\ C + A R^2 & \text{if } |y-x| \ge R, \end{cases} \qquad (38) $$

is a pointwise increasing (in the variable $R$) sequence of functions if $A \ge 0$ and it is also the case for $\exp(p_{x,C,A,R})$. It is a pointwise decreasing sequence of functions if $A \le 0$ and it is also the case for $\exp(p_{x,C,A,R})$. As a consequence, in both cases, one can apply the monotone convergence theorem to obtain

$$ \beta(p_{x,C,A}) = \lim_{R \to \infty} \beta(p_{x,C,A,R}) \le \mathcal{H}(\mathbb{P}_x|\bar{\mathbb{P}}_x). \qquad (39) $$

We can now directly optimize on the family of quadratic functions. plugging $f := C + A(y-x)^2$ in the integrals we get

$$ \beta(C + A \frac{(y-x)^2}{2}) = C + h \frac{D(x)}{2} A - \int \exp \left( C + A \frac{(y-x)^2}{2} \right) - 1 d\bar{\mathbb{P}}_x(y) $$

$(C, A) \to \beta(C + A(y-x)^2)$ is strictly concave. The optimality in $C$ gives

$$ 1 = (2\pi h \bar{a})^{-\frac{1}{2}} \exp(C) \int \exp \left( \left(A - \frac{1}{h\bar{a}}\right) \frac{(y-x)^2}{2} \right) dy. \qquad (40) $$

We identify above a Gaussian probability measure over $y$ with mean $x$ and standard deviation $\alpha$ defined by $\alpha^2(A) = (\frac{1}{h\bar{a}} - A)^{-1}$ (for small $h$ and $\bar{a} > 0$ this is always well defined). The normalizing constant gives: $(2\pi h\bar{a})^{-\frac{1}{2}} \exp(C) = (2\pi\alpha^2(A))^{-\frac{1}{2}}$. We can eliminate $C$ and the function to maximize in $A$ now is (the second integral vanishes) $A \to \frac{1}{2}\log(1 - h\bar{a}A) + h\frac{D(x)}{2}A$ the optimal $A$ satisfies $1 - h\bar{a}A = \frac{\bar{a}}{D}$ $hDA = \frac{D(x)}{\bar{a}} - 1$. Therefore we get the result using the inequality (39)

$$ \sup_{ \begin{subarray}{l} (C, A) \in (\mathbb{R} \times \mathbb{R}) \\ f := C + A(y-x)^2 \end{subarray} } \beta(f) = -\log \left( \frac{D(x)}{\bar{a}} \right) + \frac{D(x)}{\bar{a}} - 1. $$
## 5.2 Proof of theorem 10 [Convergence of Markov Chain to diffusion]

## 5.3 Convergence of Markov chains to diffusions

We consider an inhomogeneous (family in h) of transition probabilities $\mathbb{P}_{i \to i+1}^h(x, dy)$ on $\mathbb{R}^d$ for $i = 0, \dots, N-1$. We associate to it a discrete time Markov chain $X_0, \dots, X_{N-1}$ that is turned itself into a continuous time process $x(t)$ with $x(ih) = X_i$ and linear interpolation for $t \in [ih, (i+1)h)$ for $i = 0, \dots, N-1$ with law $\tilde{\mathbb{P}}^h \in \mathcal{P}(\Omega)$. The discrete sampling $x(0), x(h), \dots, x((N-1)h)$ is a Markov sequence with inhomogeneous transition

$$ \tilde{\mathbb{P}}^h(x((i+1)h) \in dy | x(ih) = x) = \mathbb{P}_{i \to i+1}^h(x, dy). $$

We define the piecewise constant in time drift and diffusion coefficients:

$$ \left\{ \begin{array}{l} \tilde{b}_t^h(x) := \frac{1}{h} \mathbb{E}_{\mathbb{P}_{[t-h-1] \to [t-h-1] + 1}}^h [(X_{i+1} - x)1_{\|X_{i+1} - x\| \le 1}], \\ \tilde{a}_t^h(x) := \frac{1}{h} \mathbb{E}_{\mathbb{P}_{[t-h-1] \to [t-h-1] + 1}}^h [(X_{i+1} - x)(X_{i+1} - x)^* 1_{\|X_{i+1} - x\| \le 1}]. \end{array} \right. \quad (41) $$

We also need a conditional moment of order larger than 4 ($\alpha > 0$) referred to abusively as “kurtosis” (kurtosis corresponds to $\alpha = 0$)

$$ c_i^h(x_i) := h^{-(2+\alpha)} \mathbb{E}_{\mathbb{P}_{i \to i+1}^h} (\|X_{i+1} - x_i\|^{4+2\alpha}). \quad (42) $$

Note that when $X_i^h$ is the discrete sampling of a continuous diffusion with bounded characteristics, $c_i^h(x_i)$ is controlled uniformly by the bound on the characteristics. (see Lemma 11).

**Theorem 10.** Assuming $(\tilde{b}_t^h, \tilde{a}_t^h)$ defined in (41) and $(c_i^h)$ defined in (42) satisfy

$$ \left\{ \begin{array}{l} i) \lim_{h \searrow 0} \sup_{|x| \le R, t \in [0,1]} \|\tilde{b}_t^h(x) - \tilde{b}_t^0(x)\| = 0, \\ ii) \lim_{h \searrow 0} \sup_{|x| \le R, t \in [0,1]} \|\tilde{a}_t^h(x) - \tilde{a}_t^0(x)\| = 0, \\ iii) \lim_{h \searrow 0} \max_{i=1, \dots, N-1} \sup_{\|x\| \le R} h c_i^h(x) = 0, \end{array} \right. \quad (43) $$

for some pair $(\tilde{b}_t^0(x), \tilde{a}_t^0(x))$ of integrable functions. Then, $\tilde{\mathbb{P}}^h$ narrowly converges to $\tilde{\mathbb{P}}^0$ a diffusion process in the weak sense with drift and diffusion coefficients $(\tilde{b}_t^0, \tilde{a}_t^0)$.

This is a generalisation of [26, Th. 11.2.3] to inhomogeneous in time Markov chain. We give here some remarks about space-time transformations in the inhomogeneous case that are not detailed in [26].

Suppose we have a inhomogeneous (family in h) of transition probabilities $\mathbb{P}_{i \to i+1}^h(x, dy)$ on $\mathbb{R}^d$ for $i = 0, \dots, N-1$. We associate to it a discrete-time Markov chain $X_0, \dots, X_{N-1}$ that is turned itself into a continuous time process $x(t)$ with $x(ih) = X_i$ and linear interpolation for $t \in [ih, (i+1)h)$ for $i = 0, \dots, N-1$. The discrete sampling $x(0), x(h), \dots, x((N-1)h)$ is a Markov sequence with inhomogeneous transition

$$ \mathbb{P}(x((i+1)h) \in dy | x(ih) = x) = \mathbb{P}_{i \to i+1}^h(x, dy). $$

We next consider an equivalent homogeneous Markov process/chain via a simple time/space transformation. Define the deterministic process over the integers $N(i) = i$ extended into a continuous time (deterministic) process $n(ih) = N(i)h = ih$ over the times of the form $ih$ and with (trivial) linear interpolation so that simply $n(t) = t$. Define the new process
$y(t) = (n(t), x(t))$. We claim that $y_t$ is a homogeneous Markov sequence on the times of the form ih with state space $[0, 1] \times \mathbb{R}^d$ (actually $\{kh, k = 0, \dots, N-1\} \times \mathbb{R}^d$ but we embed everything into $[0, 1] \times \mathbb{R}^d$). Indeed

$$
\begin{aligned}
& \mathbb{P}(y((i+1)h) \in dvdy \mid y(ih) = (u, x)) \\
& = \mathbb{P}(n((i+1)h) \in dv, x((i+1)h) \in dy \mid n(ih) = u, x(ih) = x)) \\
& = \delta_{u+h}(dv)\mathbb{P}_{[uh^{-1}] \to [uh^{-1}] + 1}^{h}(x, dy) =: \Pi^{h}((u, x), dvdy),
\end{aligned}
$$

and the last expression does not depend on $i$ hence the homogeneity. With the same notation as Stroock-Varadhan, this defines the family of (homogeneous) transitions over which we are going to apply the standard homogeneous Stroock-Varadhan result. It suffices to show that our condition (43) entails the three classical conditions of Stroock-Varadhan, i.e. conditions (2.4)-(2.5)-(2.6) p. 268.

Following the definition given at the end of p. 267 in Stroock-Varadhan, we have semi-
explicit expressions for the drift and diffusion coefficients, (with $x \in \mathbb{R}^d$ that becomes $(u, x) \in$
$[0, 1] \times \mathbb{R}^d$) namely

$$
\begin{aligned}
b_h((u,x)) &= \frac{1}{h} \int_{|(v,y)-(u,x)|\le 1} ((v,y)-(u,x))\Pi^h((u,x), dvdy) \\
&= \frac{1}{h} \begin{cases}
\int_{|y-x|\le 1} \int_{[0,1]} (v-u)\delta_{u+h}(dv)\mathbb{P}^h_{[uh^{-1}]\to [uh^{-1}] + 1}(x, dy) \\
\int_{|y-x|\le 1} (y-x)\int_{[0,1]} \delta_{u+h}(dv)\mathbb{P}^h_{[uh^{-1}]\to [uh^{-1}] + 1}(x, dy)
\end{cases} \\
&= \begin{cases}
\mathbb{P}^h_{[uh^{-1}]\to [uh^{-1}] + 1}(x, \mathcal{B}_{\mathbb{R}^d}(x, 1)) \\
\frac{1}{h} \int_{|y-x|\le 1} (y-x)\mathbb{P}^h_{[uh^{-1}]\to [uh^{-1}] + 1}(x, dy)
\end{cases}
\end{aligned}
$$

and conditions i) and iv) of (43) entail the convergence, for every $R > 0$:
$\lim_{h \to 0} \sup_{u \in [0,1], |x| \le R} \|b_h(u,x) - (p_u(x), \tilde{b}_u^0(x))\| = 0$ which is (2.4) of Stroock-Varadhan
with the state variable $(u,x) \in [0, 1 \times \mathbb{R}^d)$ instead of $x \in \mathbb{R}^d$. Likewise, with a slight abuse of
notation for the three components of the $2 \times 2$ symmetric matrix that defines the diffusion
matrix, we have

$$
\begin{aligned}
a_h((u,x)) &= \frac{1}{h} \left\{
\begin{aligned}
& \int_{|y-x| \le 1} \int_{[0,1]} (v-u)^2 \delta_{u+h}(dv) \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h}(x, dy) \\
& \int_{|y-x| \le 1} (y-x)^2 \int_{[0,1]} \delta_{u+h}(dv) \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h}(x, dy) \\
& \int_{|y-x| \le 1} (y-x) \int_{[0,1]} (v-u) \delta_{u+h}(dv) \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h}(x, dy)
\end{aligned}
\right. \\
&= \left\{
\begin{aligned}
& h \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h}(x, \mathcal{B}_{\mathbb{R}^d}(x, 1)) \\
& \frac{1}{h} \int_{|y-x| \le 1} (y-x)^2 \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h}(x, dy) \\
& \int_{|y-x| \le 1} (y-x) \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h}(x, dy)
\end{aligned}
\right.
\end{aligned}
$$
$$ \rightarrow_{h \to 0} \begin{cases} 0 \\ \frac{1}{h} \int_{|y-x| \le 1} (y-x)^2 \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1}^{h} (x, dy) \\ 0 \end{cases} $$

and conditions ii) of (43) entail the convergence $\lim_{h \to 0} \sup_{u \in [0,1], |x| \le R} \|a_h(u, x) - \begin{pmatrix} \tilde{a}^0(x) & 0 \\ 0 & 0 \end{pmatrix}\| = 0$ which is (2.5) of Stroock-Varadhan. Finally, we check (2.6) of Stroock-Varadhan, namely for $\varepsilon > 0$ given

$$ \lim_{h \to 0} \sup_{|x| \le R, u \in [0,1]} \Delta_h^\varepsilon((u, x), [0, 1] \times \mathbb{R}^d \setminus \mathcal{B}_{[0,T]} \times \mathbb{R}^d((u, x), \varepsilon)) = 0. \quad (44) $$

Setting $\mathcal{B}^c(u, x, \varepsilon) = [0, 1] \times \mathbb{R}^d \setminus \mathcal{B}_{[0,T]} \times \mathbb{R}^d((u, x), \varepsilon)$, we have

$$ \begin{aligned} \Delta_h^\varepsilon((u, x), \mathcal{B}^c(u, x, \varepsilon)) &= \frac{1}{h} \Pi^h((u, x), \mathcal{B}^c(u, x, \varepsilon)) \\ &= \frac{1}{h} \delta_{u+h}([\varepsilon, 1]) \mathbb{P}_{[uh^{-1}] \to [uh^{-1}]+1} (x, \mathbb{R}^d \setminus \mathcal{B}_{\mathbb{R}^d}(x, \varepsilon)) \\ &\le \frac{1}{h} \sup_{x \in \mathbb{R}^d} \max_{0 \le i \le N-1} \mathbb{P}_{i \to i+1} (x, \mathbb{R}^d \setminus \mathcal{B}_{\mathbb{R}^d}(x, \varepsilon)) \\ &\le \frac{1}{h} \sup_{x \in \mathbb{R}^d} \max_{0 \le i \le N-1} \mathbb{E}_{\mathbb{P}_{i,i+1}} [|X_{i+1} - X_i|^{4+2\alpha} |X_i = x|] \\ &\le h^\alpha \varepsilon^{-4-2\alpha} \sup_{x \in \mathbb{R}^d} \max_{0 \le i \le N-1} h c_i^h(x). \end{aligned} $$

Giving (44) using the definition of $c_i^h$ (42) and the Kurtosis bound constraint in $(\mathcal{V}^h)$ which enforces the uniform boundedness of $hc_i^h(x)$.

## 5.4 Proof of Theorem 3 [Well posedness for the discrete problem]

We give hereafter the proof of the existence of a solution to $(\mathcal{V}^h)$. This solution is unique due to the entropy, in contrast with the continuous problem.

The cost functional in $(\mathcal{V}^h)$ is (12):

$$ \mathcal{I}^h(\mathbb{P}^h) := h \sum_{i=0}^{N-1} \mathbb{E}_{\mathbb{P}_i^h} \left( F(b_i^h(X_i^h), a_i^h(X_i^h)) \right) + h \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) + \mathcal{D}(\mathbb{P}_0^h, \rho_0) + \mathcal{D}(\mathbb{P}_1^h, \rho_1) $$

with $(b_i^h, a_i^h)$, the discrete drift and quadratic variation increments defined in (9). The strict convexity and lower semi-continuity of $\mathcal{I}^h$ are not immediately seen but obtained, as often in optimal transport, using a linear change of variable in the conditional moments $(b_i^h(X_i^h), a_i^h(X_i^h))$ uncovering the composition of the convex perspective function associated
to F and a linear operator. Using (9), simplifying and abusing notations, for $\beta = 1, 2, ...$:

$$ 
\begin{aligned}
\mathbb{P}^h & \to \mathbb{E}_{\mathbb{P}_i^h} \left( F\left(\frac{1}{h} \mathbb{E}_{\mathbb{P}_{i \to i+1}^h} \left( (X_{i+1}^h - X_i^h)^\beta \right) \right) \right) \text{ is decomposed as} \\
\mathbb{P}^h & \to (\mathbb{P}_i^h(x_i), \frac{1}{h} \mathbb{E}_{\mathbb{P}_{i, i+1}^h}(x_i, \dots) \left( (X_{i+1}^h - x_i)^\beta \right)), \quad \forall x_i \\
& \to \mathbb{E}_{\mathbb{P}_i^h} \left( F\left(\frac{1}{h} \frac{\mathbb{E}_{\mathbb{P}_{i, i+1}^h}(X_i^h, \dots) \left( (X_{i+1}^h - X_i^h)^\beta \right)}{\mathbb{P}_i^h(X_i^h)}\right) \right).
\end{aligned}
\quad (45) $$

where $\mathbb{P}_{i, i+1}^h(x_i, \dots)$ is to be understood as the measure on $\mathbb{R}_{t_{i+1}}$ obtained by freezing the first variable $x_i$ in the joint $\mathbb{R}_{t_i} \times \mathbb{R}_{t_{i+1}}$ probability $\mathbb{P}_{i, i+1}^h$.

The convexity of additional Kurtosis constraint is also a consequence of (45). The optimization problem ($\mathcal{V}^h$) therefore has a unique minimizer (the relative entropy is strictly convex).

The Markovianity of the minimiser a consequence of the structure of $\mathbb{I}^h$ and the additivity properties of the relative entropy giving

$$ \mathcal{H}(\mathbb{P}^h | \mathbb{P}^h) \ge \sum_{i=0}^{N-1} \mathbb{E}_{\mathbb{P}_i^h} \left( \mathcal{H}(\mathbb{P}_{i \to i+1}^h | \mathbb{P}_{i \to i+1}^h) \right) \text{ with equality if } \mathbb{P}^h \text{ is Markov} \quad (46) $$

(see [11] lemma 3.4).

## 5.5 Proof of Lemma 5 [Diffusion Coefficients rescaling]

**Scaling.** We start with an approximation step by scaling the volatility. Consider a solution $X(t)$ of the continuous problem with a finite cost, then introduce $X^{\alpha, \varepsilon}(t) = \sqrt{\alpha}X(t) + \sqrt{\varepsilon}B(t)$, where $B(t)$ is a $d$-dimensional Brownian motion, independent of $X$. As a consequence, $a_t^{\alpha, \varepsilon} = \alpha a_t + \varepsilon \text{Id}$. Since $a \in [\lambda \text{Id}, \Lambda \text{Id}]$, it implies that $(\alpha \lambda + \varepsilon) \text{Id} \le a_{\alpha, \varepsilon} \le (\alpha \Lambda + \varepsilon) \text{Id}$. Now, because $\lambda < \Lambda$, and for $\delta$ small enough, a simple computation shows that taking $\alpha = 1 - \frac{2\varepsilon}{\lambda+\Lambda}$ and $\varepsilon = \delta \frac{\lambda+\Lambda}{\Lambda-\lambda}$ with $c_1, c_2 > 0$ there holds $\lambda + \delta \le a_{\alpha, \varepsilon} \le \Lambda - \delta$ and $\sup |b_{\alpha, \varepsilon}| = \sqrt{\alpha} \sup |b| \text{ note that } \alpha < 1$. Moreover, when $\delta \to 0$, $\alpha \to 1$, $\varepsilon \to 0$. Now, remark that the boundary conditions at final time 1 is lost. The distribution at time $i = 0, 1$ of $Y$ is $g_{t, \varepsilon} \star [T_{\sqrt{\alpha}}]_\#(\rho(i))$, where $g_{t, \varepsilon}$ is the gaussian kernel of variance $\varepsilon t \text{Id}$ and $T_{\sqrt{\alpha}}(x) = \sqrt{\alpha}x$ is a rescaling. The same computations apply to a discrete process $X^h$, and its characteristics $a_h, b_h$.

**Bounds on the characteristics of the discretized process.** When we consider a time continuous discretisation a continuous process satisfying the bounds on $a, b$ we might lose the bounds. However, performing first the above rescaling gives enough margin to have the discretized process satisfy the bounds. There holds

$$ \mathbb{E}[X^{\alpha, \varepsilon}(t_{i+1}) - X^{\alpha, \varepsilon}(t_i) | X(t_i)] = \mathbb{E} \left[ \int_{t_i}^{t_{i+1}} b_t^{\alpha, \varepsilon} dt \right], \quad (47) $$

and

$$ \mathbb{E}[(X^{\alpha, \varepsilon}(t_{i+1}) - X^{\alpha, \varepsilon}(t_i))^2 | X(t_i)] = \mathbb{E} \left[ \int_{t_i}^{t_{i+1}} 2(X^{\alpha, \varepsilon}(t) - X^{\alpha, \varepsilon}(t_i)b_t^{\alpha, \varepsilon}) + \text{tr}(a_t^{\alpha, \varepsilon})dt \right] \quad (48) $$

$$ = \mathbb{E} \left[ \int_{t_i}^{t_{i+1}} \text{tr}(a_t^{\alpha, \varepsilon}) dt \right] + C(\Lambda, B)(h^{3/2}). \quad (49) $$
Therefore, for h small one can choose $\delta_h$ to rescale the process such that the bound holds for the characteristics of the rescaled discrete process $a^{\alpha, \varepsilon, h}$, $b^{\alpha, \varepsilon, h}$.

**Convergence of F.** This follows simply by observing that

$$ \mathbb{E}_Y \int_t F(a_Y, b_Y) = \mathbb{E}_Y \int_t G(a_Y, b_Y) \\ = \mathbb{E}_{X,B} \int_t G(a_Y, b_Y) = \mathbb{E}_{X,B} \int_t G(\alpha a_X + \varepsilon, \sqrt{\alpha} b_X) \quad (50) $$

since G is uniformly Lipschitz in its domain, the estimate on the functional I follows.

**Kurtosis.** The Kurtosis bound (42) is preserved whenever $\alpha < 1$ and $\varepsilon$ is small enough depending on $\alpha$.

**Entropy, continuous case.** In the time continuous case, the entropy term can be put into F, as a Lipschitz function of a, thus the convergence is shown as in (50).

**Entropy in the discrete case.** Let us first denote by $\mathbb{P}_\varepsilon^h$ the law of $X + \sqrt{\varepsilon}B$, $B$ a standard Brownian motion. We first use the change of variable

$$ X = (x_1, \dots, x_N) \mapsto Z = (x_1, x_2 - x_1, \dots, x_N - x_{N-1}). $$

Note that, under this change of coordinates, $\mathbb{P}_\varepsilon^h$ is the convolution of $\mathbb{P}^h$ with a (diagonal) Gaussian kernel G defined on $\mathbb{R}^{dN}$ such that $\int_X \|X\|^2 G(X) dX = \varepsilon$. We first have that

$$ \mathcal{H}(\mathbb{P}_\varepsilon^h | \bar{\mathbb{P}}^h) \le \int_{\mathbb{R}^{dN}} \mathcal{H}(\mathbb{P}^h | [T_X]_\sharp \bar{\mathbb{P}}^h) G(X) dX, $$

where $T_X$ denotes the translation by $X \in \mathbb{R}^{dN}$. This inequality is obtained thanks to the convexity of the entropy in the first variable. Let us denote by $\rho$ the image measure of $\mathbb{P}^h$ and $\rho_0$ the image measure of $\bar{\mathbb{P}}^h$ on $\mathbb{R}^{dN}$,

$$ \mathcal{H}(\mathbb{P}^h | [T_X]_\sharp \bar{\mathbb{P}}^h) = \int_z \rho(z) \log (\rho(z) / \rho_0(z + X)) dz. $$

Using the fact that $\rho_0$ is Gaussian, we have $\rho_0(z + X) = \exp(-2N\langle X, z \rangle + N\|X\|^2) \rho_0(z)$. We get

$$ \mathcal{H}(\mathbb{P}^h | [T_X]_\sharp \bar{\mathbb{P}}^h) = \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) + \int_z \rho(z) (-2N\langle X, z \rangle + N\|X\|^2) dz. $$

Integrating over X with respect to a the centered variable gives $\int_X \int_Z \rho(z) (-2N\langle X, z \rangle dz G(X) dX = 0$. This implies

$$ \int_X \mathcal{H}(\mathbb{P}^h | [T_X]_\sharp \bar{\mathbb{P}}^h) G(X) dX = \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) + \int_X N \|X\|^2 G(X) dX. $$

Thus, we get $h\mathcal{H}(\mathbb{P}_\varepsilon^h | \bar{\mathbb{P}}^h) \le h\mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) + \varepsilon$. We now treat similarly the scaling in $\alpha$,

$$ \begin{aligned} \mathcal{H}([S_{\sqrt{\alpha}}]_\sharp \mathbb{P}^h | \bar{\mathbb{P}}^h) &= \int_z \rho(z) \log \left( \alpha^{-dN/2} \rho(z) / \rho_0(\sqrt{\alpha}z) \right) dz \\ &= -dN/2 \log(\alpha) - \int_z \rho(z) \log (\rho_0(\sqrt{\alpha}z)) dz + \int_z \rho(z) \log (\rho(z)) dz \\ &= -dN \log(\alpha)/2 + \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) - \int_z \rho(z) 2N(\alpha - 1) \|z\|^2 dz. \end{aligned} \quad (51) $$

We obtain

$$ h\mathcal{H}(\mathbb{P}^h | [S_{\sqrt{\alpha}}]_\sharp \bar{\mathbb{P}}^h) = -d\log(\alpha)/2 + h\mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h) - \int_z \rho(z) 2(\alpha - 1) \|z\|^2 dz. $$

The quantity $\int_z \|z\|^2 \rho(z) dz$ is equal to the sum of the second moment of $\rho(t=0)$ and the quadratic variation of $\mathbb{P}^h$ which is finite uniformly in $h$. This completes the proof of the Lemma.

<footnote n="1">¹Note that this choice of regularization corresponds to a fixed regularization with a Gaussian kernel on the time interval $[0, 1]$ of variance $\varepsilon$.</footnote>
**Remark 6.** Another possibility for this proof is to apply the rescaling and regularisation also to the reference measure. Then $\mathcal{H}(\mathbb{P}_\varepsilon^h | \bar{\mathbb{P}}_\varepsilon^h) \le \mathcal{H}(\mathbb{P}^h | \bar{\mathbb{P}}^h)$ is straightforward. We then have to deal a specific relative entropy with a reference diffusion which depends on $\varepsilon$ via $\bar{a}_\varepsilon = \alpha\bar{a} + \sigma (\alpha \text{ and } \sigma \text{ depend on } \varepsilon)$. Then (20) becomes $\mathcal{I}_\varepsilon^0(\mathbb{P}_\varepsilon^0) \le \liminf_{h \searrow 0} \mathcal{I}^h(\mathbb{P}^h) + O(\varepsilon)$ where $\mathcal{I}_\varepsilon^0 = \mathcal{I}^0 + R_\varepsilon$ with $R_\varepsilon = \mathbb{E}_{\bar{\mathbb{P}}_\varepsilon^0}(S^{\mathcal{I}}(\tilde{a}_\varepsilon^0 | \bar{a}_\varepsilon) - S^{\mathcal{I}}(\tilde{a}_\varepsilon^0 | \bar{a})) = O(\bar{a}_\varepsilon - \bar{a}) = O(\varepsilon)$.

## 5.6 Proofs of Lemma 6 [Regularization time continuous case]

The first step concerns space regularization and consists of adding a Gaussian variable of variance $\sigma$ to the process $X(t)$. The corresponding density at time 0 is given by $g_\sigma * \rho(0)$ where $\rho(0)$ is the marginal at time 0 of $X(t)$.

The second step involves extrapolating the process $X$ in time with a process that has a finite cost. Choose $\sigma_0 \in [\lambda, \Lambda]$ and define on $[-2\varepsilon, 0]$, the solution of the heat equation starting at time $t = -2\varepsilon$ with the initial condition $g_{\sigma - (2\varepsilon\sigma_0)} * \rho(0)$. Here we need to assume $2\varepsilon\sigma_0 < \sigma$ which is always satisfied for $\varepsilon$ small enough. At time $t = 0$, the density equals $g_\sigma * \rho(0)$. For the times $t > 1$, we also use a diffusion with coefficient $\sigma_0$ so that $X(t)$ can be extended to the time interval $[-2\varepsilon, 1 + 2\varepsilon]$ with a finite cost outside the interval $[0, 1]$ of order $O(\varepsilon)$.

Let $k: \mathbb{R} \to \mathbb{R}_+$ be a smooth nonnegative function with support in the unit ball, such that $\int k(y)dy = 1$. We denote the kernel $\eta_\varepsilon(t, x) = k(t/\varepsilon)/\varepsilon$ with support in the time variable is thus contained in the ball of radius $\varepsilon$.

We consider

$$ \mathbb{P}_{t, \varepsilon}^0 := (\mathbb{P}_t^0 * \eta_\varepsilon) \qquad b_{t, \varepsilon}^0 = \frac{(b_t^0 d\mathbb{P}_t^0) * \eta_\varepsilon}{d\mathbb{P}_{t, \varepsilon}^0} \qquad a_{t, \varepsilon}^0 = \frac{(a_t^0 d\mathbb{P}_t^0) * \eta_\varepsilon}{d\mathbb{P}_{t, \varepsilon}^0}. \qquad (52) $$

By convexity, the drift and diffusion coefficients satisfy the hard constraints of the lemma's hypotheses. Moreover, these coefficients (drift and diffusion) are smooth in time (due to the convolution with $\eta_\varepsilon$) and space (due to the initial regularization in the first step) so that by standard arguments, there exists a unique Markov process solving the corresponding diffusion equation which is well defined on the time interval $[-\varepsilon, 1 + \varepsilon]$. We denote such a process by $X^\varepsilon(t)$. By the control on the time extrapolation and the convexity of the functional $\mathcal{I}^0$, we have

$$ \int_{-\varepsilon}^{1+\varepsilon} \int_x G(b_{t,\varepsilon}(x), a_{t,\varepsilon}(x)) d\rho_{t,\varepsilon}(x) dt \le \mathcal{I}^0(\rho_t^0, b_t^0, a_t^0) + O(\varepsilon). $$

We rescale the time interval from $[-\varepsilon, 1 + \varepsilon]$ to $[0, 1]$ to control the soft constraint on the boundary terms. Due to the hypotheses (hard constraints on the drift and diffusion coefficients), $Y := X^\varepsilon((1 + 2\varepsilon)(t + \varepsilon))$ satisfies the hard constraints on the drift $|b_{t,\varepsilon}| \le B$ and volatility $\lambda \le a_{t,\varepsilon}$ Id $\le \Lambda$. Importantly, the distributions of $Y(0)$ and $Y(1)$ are (infinitesimal) mixture of Gaussian convolutions of $\rho(0)$ and $\rho(1)$ the marginals of the process $X(t)$, since the kernel has compact support in time in $[-\varepsilon, \varepsilon]$ and the extension in time is by heat diffusion. As a consequence and from the hypotheses on $\mathcal{D}$, one has for $i = 0, 1$, $|\mathcal{D}(\mathbb{P}_{t=i}^Y, \rho_i) - \mathcal{D}(\mathbb{P}_{t=i}^X, \rho_i)| \le O(\varepsilon) + O(\sigma)$ (where we recall that $\rho_i$ are the (soft) boundary values). Thus, we obtain the result $\mathcal{I}(\mathbb{P}^Y) \le \mathcal{I}(\mathbb{P}^X) + O(\varepsilon) + O(\sigma)$.

## 5.7 Proof of Lemma 7

For the specific entropy, Proposition 2 i) gives for the entropic part of $\mathcal{I}(\mathbb{P}_\varepsilon^h)$:

$$ \lim_{h \searrow 0} h \mathcal{H}(\mathbb{P}_\varepsilon^h | \bar{\mathbb{P}}^h) = S^\mathcal{H}(\mathbb{P}_\varepsilon^0 | \bar{\mathbb{P}}). \qquad (53) $$
For the integrand part, we first start with (already in the proof of lemma 5 iv))

$$
\begin{align*}
\mathbb{E}[X(t_{i+1}) - X(t_i) | X(t_i)] &= \mathbb{E}_{\mathbb{P}_{i \to i+1, \varepsilon}^h} \left[ \int_{t_i}^{t_{i+1}} b_{t, \varepsilon}(X_t) dt + \int_{t_i}^{t_{i+1}} \sqrt{a_{t, \varepsilon}(X_t)} dB_t \right] \\
&= \mathbb{E}_{\mathbb{P}_{i \to i+1, \varepsilon}^h} \left[ \int_{t_i}^{t_{i+1}} b_{t, \varepsilon}(X_t) dt \right],
\end{align*}
$$

$$
\begin{align*}
\mathbb{E}[(X(t_{i+1}) - X(t_i))^2 | X(t_i)] &= \mathbb{E}_{\mathbb{P}_{i \to i+1, \varepsilon}^h} \left[ \int_{t_i}^{t_{i+1}} 2(X(t) - X_i^h) b_{t, \varepsilon}(X(t)) + \text{tr}(a_{t, \varepsilon}(X_t)dt \right] \\
&= \mathbb{E}_{\mathbb{P}_{i \to i+1, \varepsilon}^h} \left[ \int_{t_i}^{t_{i+1}} \text{tr}(a_{t, \varepsilon}(X_t)dt \right] + O(h^{1+\alpha}),
\end{align*}
$$

where $\alpha > 0$ depends on the integrability of $b$. This formula implies that the discrete characteristic coefficients $b_i^h$ and $a_i^h$ for $\mathbb{P}_\varepsilon^h$ satisfy for $h$ small enough $\lambda Id < a < \Lambda Id$ and $|b| < B$. We also have:

$$
h \mathbb{E}_{\mathbb{P}_{i, \varepsilon}^{h}} [F(b_{i, \varepsilon}^{h}(X(t_{i})), a_{i, \varepsilon}^{h}(X(t_{i})))] =
$$

$$
h \mathbb{E}_{\mathbb{P}_{i, \varepsilon}^h} \left[ G \left( \frac{1}{h} \mathbb{E}_{\mathbb{P}_{i \to i+1, \varepsilon}^h} \left[ \int_{t_i}^{t_{i+1}} b_{t, \varepsilon}(X(t)) dt \right], \frac{1}{h} \mathbb{E}_{\mathbb{P}_{i \to i+1, \varepsilon}^h} \left[ \int_{t_i}^{t_{i+1}} \text{tr}(a_{t, \varepsilon}(X(t))) + O(h^\alpha) dt \right] \right) \right]. \quad (54)
$$

Using Jensen's inequality and the Lipschitz assumption on $G$, we get

$$
h \mathbb{E}_{\mathbb{P}_{i, \varepsilon}^{h}} [G(b_{i, \varepsilon}^{h}(X(t_i)), a_{i, \varepsilon}^{h}(X(t_i)))] \le \int_{t_i}^{t_{i+1}} \mathbb{E}_{\mathbb{P}_{i, \varepsilon}^{h}} G(b_{t, \varepsilon}(X(t), \text{tr}(a_{t, \varepsilon}(X(t))))) dt + O(h^{\alpha+1}).
$$

The control on the time boundary marginal terms follows from lemma 5 and 6.

5.8 Proof of Lemma 8 [Regularization time discrete case]

Tightness: We need tightness of the quantities $\mathbb{P}^h$, $m^h$ and $n^h$ defined by

$$
m_i^h = \frac{1}{h} \left( \int_{x_{i+1}} (x_{i+1} - x_i) \mathbb{P}_{i \to i+1}^h (x_{i+1} | x_i) \mathbb{P}_i^h (x_i) \right), \quad (55)
$$

$$
n_i^h = \frac{1}{h} \left( \int_{x_{i+1}} (x_{i+1} - x_i)^2 \mathbb{P}_{i \to i+1}^h (x_{i+1} | x_i) \mathbb{P}_i^h (x_i) \right). \quad (56)
$$

These quantities are discrete in time; we will define this quantity continuously to obtain
tightness in the space of time-space measures. By considering either the Markov chain or the
time-linearly interpolated process, we can consider $\mathbb{P}^h$ either as discrete probabilities or as
probabilities on the set of continuous paths. Let us define more explicitly the time-dependent
path of marginals $\rho^h(t) = \sum_{i=0}^N \mathbf{1}_{[t_i, t_{i+1})} \mathbb{P}_i^h$. To obtain tightness of the $\rho^h$, we use linear
interpolation of $\mathbb{P}^h$ denoted $\tilde{\mathbb{P}}^h$. We start with the tightness of $\mathbb{P}^h$. Recall that we have a
control on the $\alpha$-moments ($\alpha > 4$) of $\mathbb{P}^h$. In particular, the linearly interpolated process
$X^h(t)$ satisfies, for $s, t \in [0, T]$, $|s-t| < 1/N$, for some $c > 0$, $\mathbb{E}[\|X_s^h - X_t^h\|^\alpha] \le c|s-t|^{\alpha/2-1}$.
The Kolmogorov lemma implies that the sequence $\tilde{\mathbb{P}}^h$ only charges Hölder trajectories with
Hölder coefficient $1/2 - 2/\alpha - \varepsilon' > 0$ (positive since $\alpha > 4$) for all positive $\varepsilon'$ sufficiently
small. Hence the sequence $\tilde{\mathbb{P}}^h$ is tight and consequently, so is $\rho^h(t)$.

To apply the criterion in [26], we need to check if the "jump condition" is satisfied on $\mathbb{P}^h$.
It is guaranteed by the kurtosis bound. We now deal with quantities $m_i^h$ and $n_i^h$. We have
$$m^h = \sum_{i=0}^N \mathbf{1}_{[t_i, t_{i+1})} m_i^h, \quad n^h = \sum_{i=0}^N \mathbf{1}_{[t_i, t_{i+1})} n_i^h.$$
(defined on $[0, 1] \times \mathbb{R}^d$) $m^h$ and $n^h$ are tight. For that, we apply the De la Vallée Poussin lemma. We propose to bound the $1 < 1+\alpha < 2$ moment, using the inequality $ab \le \frac{1}{2}(a^2+b^2)$,

$$ \int_0^1 \int_{\mathbb{R}^d} |x|^{1+\alpha} n^h(t, x) dt \le h \left( \sum_{i=1}^N \frac{1}{2} \int |x_i|^{2+2\alpha} \mathbb{P}_i^h + \frac{1}{2h^2} \int |x_{i+1} - x_i|^4 \mathbb{P}^h \right). \quad (57) $$

The second term on the right-hand side is bounded due to the Kurtosis bound. It suffices to prove a uniform bound in time of $(2+2\alpha)$-moment to bound the first term, which is bounded by the 4-moment when $\alpha < 1$ which is true. Using a telescopic sum, there exists a positive constant $M > 0$ such that

$$ \int_0^1 \int_{\mathbb{R}^d} |x_i|^4 \mathbb{P}_{i,i+1}^h (x_i, x_{i+1}) dt \le M \sum_{i=1}^N (\int |x_0|^4 + \sum_{j=1}^i |x_{i+1} - x_i|^4 \mathbb{P}^h). \quad (58) $$

Our hypothesis bounds the first term on the right-hand side on $\rho_0$. The second term is bounded by $O(h)$ due to the kurtosis bound. The proof is similar for the quantity $m^h$.

**Scaling the coefficients and space regularization:** To prepare for the time extrapolation below, we transform $\mathbb{P}^h$ by using Lemma 5. Denoting $X^h(t)$ any interpolation of $\mathbb{P}^h$, we consider the transformation of the type (for well-chosen parameters $\alpha, \delta$, see Lemma 5)

$$ Y(t) = \sqrt{\alpha} X^h(t) + \sqrt{\delta} B(t) + B_0, \quad (59) $$

with $B(t)$ a standard Brownian motion and $B_0$ a centered Gaussian variable of variance $\varepsilon$ Id. Lemma 5 gives that $Y$ satisfies the estimates on the functional and the hard constraints on the drift and the volatility. The kurtosis bound is also trivially satisfied.

**Time extrapolation:** We want to again define the convolution in time and space of the probability $\mathbb{P}^h$ extended on the path space by linear interpolation of the curves. However, it is not defined outside the time interval $[0, 1]$, which is needed to regularize the process in time and control the density values at times 0, 1. We use the two first steps of the proof of Lemma 6 to construct an extrapolation on the time interval $[-2\varepsilon, 0]$ and $[1, 1+2\varepsilon]$ with finite cost $\mathcal{I}^0(\mathbb{P}) + O(\varepsilon) + O(\delta)$. We now discretize it in time as in the Gamma-Limsup proof so that the cost $\mathcal{I}^h$ on $[-2\varepsilon, 0]$ and $[1, 1+2\varepsilon]$ is of the order $O(\varepsilon) + O(h)$. This discretization gives us a time probability measure in discrete time still denoted by $\mathbb{P}^h$.

**Time regularization:** We use $\eta_\varepsilon(s, x) := k(s/\varepsilon)/\varepsilon$ as in Lemma 6. We now introduce a time-dependent curve with values in the convex set of plans, i.e. probability measures on the product space,

$$ \Pi^h(t, x, y) = \sum_i \mathbf{1}_{[t_i, t_{i+1})} \mathbb{P}_{i, i+1}^h(x, y), \quad (60) $$

which is a collection of time-dependent plans indexed by $t \in [-2\varepsilon, 1+2\varepsilon]$. To a given plan $\pi(x, y) \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d)$, one can define its conditional moments (exactly as in Formula (55)), in probabilistic notations $m^h(\pi) := \frac{1}{h} \int_y (y-x)\pi(x, dy)$, and $n^h(\pi) := \frac{1}{h} \int_y (y-x)^2\pi(x, dy)$. Now, we observe that $m^h(t) = m^h(\Pi(t))$, $n^h(t) = n^h(\Pi(t))$. In addition, the marginalization of $\Pi$ on the first variable gives $\rho^h(t)$. Note also that the following property holds by definition $\forall t$ $[p_2]_\Pi^h(t) = [p_1]_\Pi^h(t+h)$, where the $[p_i]_\Pi^h$, $i = 1, 2$ is the pushforward operator on the first and second variables. We now regularize in time this object $\Pi_\varepsilon^h = \eta_\varepsilon \star \Pi^h$. Importantly, the regularization on $\Pi^h$ induces a regularization on the quantities $\rho, m, n$: By linearity, we have that

$$ [\pi_1]_\Pi^h = \eta_\varepsilon \star \rho^h \quad m^h(\Pi_\varepsilon^h) = \eta_\varepsilon \star m^h \quad n^h(\Pi_\varepsilon^h) = \eta_\varepsilon \star n^h. $$

In addition, it defines a discrete Markov chain: By linearity, we have the property: $\forall t$ $[p_2]_\Pi^h(t) = [p_1]_\Pi^h(t+h)$. From $\Pi_\varepsilon$ we construct a probability on the path space: we
evaluate $\Pi_{\varepsilon}(ih)$ which is a plan and we note that $[p_2]_{\sharp}\Pi_{\varepsilon}(ih) = [p_1]_{\sharp}\Pi_{\varepsilon}((i+1)h)$, by the property mentioned above. By the standard gluing lemma [35], we construct a corresponding Markov chain using the joint probabilities $\Pi_{\varepsilon}(ih)$ that we denote $\mathbb{P}_{\delta, \varepsilon}^h$. By convexity of the integral part of the functional, we estimate the cost

$$ \mathcal{I}^h(\mathbb{P}_{\delta, \varepsilon}^h) \le \mathcal{I}^h(\mathbb{P}^h) + O(\varepsilon) + O(h) + O(\delta), \quad (61) $$

where the term $O(\varepsilon) + O(h)$ is due to the extension in time of $\mathbb{P}^h$ and its time discretization. On $\rho_{\delta, \varepsilon}^h, m_{\delta, \varepsilon}^h, n_{\delta, \varepsilon}^h$ which are smoothing of a tight sequence of measures, we can extract a subsequence that is bounded in $L_{t,x}^{\infty}$ and smooth in both time and space, uniformly converging on compact subsets.

**Time-rescaling:** Similarly to Lemma 6, we rescale time so that the quantities are all defined on the $[0, 1]$. The cost of this rescaling is $O(\varepsilon)$. Note that by rescaling, the discretization timestep is also rescaled, by abuse of notation, we still keep the same letter $h$. More importantly, the rescaling increases the volatility linearly in $\varepsilon$ and we need to check that the hard constraints on the volatility are preserved. This is indeed the case since Lemma 5 has been applied beforehand.

**Riemann integral:** Recall the functional $\mathcal{J}$ as having three arguments

$$ \mathcal{J}(\rho, m, n) = \int_{t,x} F(m/\rho, n/\rho) + S_{\bar{a}}^{\mathcal{I}}(n/\rho) d\rho(t, x)dt, $$

where $S_{\bar{a}}^{\mathcal{I}}$ denotes the specific entropy with reference diffusion $\bar{a}$ using $n, \rho$ as arguments. Note that we want to study $\mathcal{J}(\rho(\Pi(t)), m^h(\Pi(t)), n^h(\Pi(t)))$. Then, the discrete formulation $\mathcal{I}^h$ (as well as the specific entropy) can be viewed as a Riemann integral in time:

$$ \mathcal{I}^h(\mathbb{P}^h) \ge \sum_{i=0}^{N} h \int_x F(m^h(\Pi(ih))/\rho(\Pi(ih)), n^h(\Pi(ih))/\rho(\Pi(ih))) \\ + S_{\bar{a}}^{\mathcal{I}}(n^h(\Pi(ih))/\rho(\Pi(ih)))d\rho(\Pi(ih)), \quad (62) $$

where $\rho(\Pi(ih))$ is the marginal of $\Pi(ih)$ on the first variable. In addition, we will use below that this Riemann sum coincides with the functional $\mathcal{J}$ evaluated on $(\tilde{\rho}, \tilde{m}, \tilde{n})$. By regularization, $(m_{\delta, \varepsilon}^h(t), n_{\delta, \varepsilon}^h(t), \rho_{\delta, \varepsilon}^h(t))$ is smooth in time. It is also the case for the function integrated in space. To prove it, we use that the drift and volatility coefficients are bounded in $L^{\infty}$. We have

$$ \begin{aligned} & \left| \int_x F(m_{\delta, \varepsilon}^h(t)/\rho_{\delta, \varepsilon}(t), n_{\delta, \varepsilon}^h(t)/\rho_{\delta, \varepsilon}(t)) + S_{\bar{a}}^{\mathcal{I}}(m_{\delta, \varepsilon}^h(t)/\rho_{\delta, \varepsilon}(t), n_{\delta, \varepsilon}^h(t)/\rho_{\delta, \varepsilon}(t))d\rho_{\delta, \varepsilon}(t) \right. \\ & \qquad \left. - F(m_{\delta, \varepsilon}^h(t')/\rho_{\delta, \varepsilon}(t'), n_{\delta, \varepsilon}^h(t')/\rho_{\delta, \varepsilon}(t')) + S_{\bar{a}}^{\mathcal{I}}(m_{\delta, \varepsilon}^h(t')/\rho_{\delta, \varepsilon}(t'), n_{\delta, \varepsilon}^h(t')/\rho_{\delta, \varepsilon}(t'))d\rho_{\delta, \varepsilon}(t') \right| \\ & \le M \int_x |\rho_{\delta, \varepsilon}(t') - \rho_{\delta, \varepsilon}(t)| \le MC(\varepsilon)|t - t'|, \quad (63) \end{aligned} $$

where $M$ bounds the integrand (which is possible since all the terms are bounded independently of $h, \varepsilon, \delta$) and $C(\varepsilon)$ is a constant coming from the convolution kernel and depending on $\varepsilon$. Therefore, we get

$$ \mathcal{I}^h(\mathbb{P}_{\delta, \varepsilon}^h) \ge \mathcal{J}(\tilde{\rho}_{\delta, \varepsilon}^h, \tilde{m}_{\delta, \varepsilon}^h, \tilde{n}_{\delta, \varepsilon}^h) \ge \mathcal{J}(\rho_{\delta, \varepsilon}^h, m_{\delta, \varepsilon}^h, n_{\delta, \varepsilon}^h) + O(h), $$

which gives the inequality (19). Note that the $O(h)$ is due to the difference between the Riemann sum and the integral. By lower-semicontinuity of $\mathcal{J}$, one has $\mathcal{J}(\rho_{\delta, \varepsilon}, m_{\delta, \varepsilon}, n_{\varepsilon}) \le \lim_{h \to 0} \mathcal{J}(\rho_{\delta, \varepsilon}, m_{\delta, \varepsilon}, n_{\varepsilon})$. In particular, the two previous inequalities imply iii) of the lemma and $\mathcal{J}(\rho_{\delta, \varepsilon}, m_{\delta, \varepsilon}, n_{\varepsilon}) \le \lim_{h \to 0} \mathcal{I}^h(\mathbb{P}_{\delta, \varepsilon}^h)$.
**Conclusion:** The coefficients $m_{\delta,\varepsilon}^h, n_{\delta,\varepsilon}^h$, and $\rho_{\delta,\varepsilon}^h$ converge uniformly in space and time: It implies that $b_h(ih) = \frac{m_{\delta,\varepsilon}^h(ih)}{\rho_{\delta,\varepsilon}^h(ih)}$ and $a_h(ih) = \frac{n_{\delta,\varepsilon}^h(ih)}{\rho_{\delta,\varepsilon}^h(ih)}$ converge uniformly on every compact in time and space since $\rho_\varepsilon$ is bounded below by a positive constant on every compact set. Thus, we can apply Theorem 10 to obtain ii) of the lemma: the sequence $\mathbb{P}_\varepsilon^h$ uniformly converges to $\mathbb{P}_\varepsilon$ the probability on the path space of a diffusion process. By uniform convergence, the drift and diffusion coefficients can be identified by $b_{\delta,\varepsilon}(t,x) = \frac{m_{\delta,\varepsilon}(t,x)}{\rho_{\delta,\varepsilon}^h(t,x)}$ and $a_{\delta,\varepsilon}(t,x) = \frac{n_{\delta,\varepsilon}(t,x)}{\rho_{\delta,\varepsilon}^h(t,x)}$.
Finally, we get

$$ \mathcal{J}(\rho_{\delta,\varepsilon}, m_{\delta,\varepsilon}, n_{\delta,\varepsilon}) \leq \liminf_{h \to 0} \mathcal{I}^h(\mathbb{P}^h) + O(\varepsilon) + O(\delta). $$

## 5.9 Additional Lemmas

**Lemma 11.** The quantities $b_t, a_t$ are the characteristic coefficients of the semi-martingale $X$ with law $\mathbb{P} \in \mathcal{P}^1$. For every $p \ge 1$ and $q, q' > 1$, we have, in full generality:

$$ \mathbb{E}[|X_t - X_s|^p] \le C_p \left( (t-s)^{p(1-\frac{1}{q})} \left( \mathbb{E}\left[\int_0^1 |b_u(X_u)|^{\max(p,q)} du\right] \right)^{\min(p/q,1)} \right. \\ \left. \qquad + (t-s)^{\frac{p}{2}(1-\frac{1}{q'})} \left( \mathbb{E}\left[\int_0^1 |a_u(X_u)|^{\max(p/2,q')} du\right] \right)^{\min(p/2q',1)} \right), \quad (64) $$

where $C_p > 0$ only depends on $p$.

*Proof of the lemma.* In the following, the notation $C_p$ denotes a positive number that only depends on $p$ and that may vary at each occurrence. We start with the following observation: for any measurable random process $f_u$, we have, for any $q > 1$, by Hölder's inequality

$$ \mathbb{E}\left[\left|\int_s^t f_u du\right|^p\right] \le (t-s)^{p(1-\frac{1}{q})} \begin{cases} \mathbb{E}\left[\int_0^1 |f_u|^p du\right] & \text{if } p \ge q \\ \mathbb{E}\left[\int_0^1 |f_u|^q du\right]^{p/q} & \text{if } p < q \end{cases} \\ = (t-s)^{p(1-\frac{1}{q})} \left(\mathbb{E}\left[\int_0^1 |f_u|^{\max(p,q)} du\right]\right)^{\min(1,p/q)} \quad (65) $$

Next, from $X_t - X_s = \int_s^t b_u(X_u) du + \int_s^t a_u(X_u)^{1/2} dB_u$ and thus $\mathbb{E}[|X_t - X_s|^p] \le C_p(\mathbb{E}[\int_s^t b_u(X_u) du]^p + \mathbb{E}[\int_s^t a_u(X_u)^{1/2} dB_u]^p)$ since $p \ge 1$, we successively have $\mathbb{E}[\int_s^t b_u(X_u) du]^p \le (t-s)^{p(1-\frac{1}{q})} (\mathbb{E}[\int_0^1 |b_u(X_u)|^{\max(p,q)} du])^{\min(1,p/q)}$ and

$$ \mathbb{E}\left[\left|\int_s^t a_u(X_u)^{1/2} du\right|^p\right] \le C_p \mathbb{E}\left[\left|\int_s^t a_u(X_u) du\right|^{p/2}\right] \\ \le (t-s)^{\frac{p}{2}(1-\frac{1}{q'})} \left(\mathbb{E}\left[\int_0^1 |a_u(X_u)|^{\max(p/2,q')} du\right]\right)^{\min(1,p/2q')} \quad (66) $$

by (65) and the Burkholder-Davis-Gundy inequality. The conclusion follows. $\square$

A brief excerpt from chapter 9 of [35] and [36, 37]:

**Definition 1** (perspectives functions over measures a.k.a. entropies). Assuming that $F(\cdot) = \theta(|\cdot|) : \mathbb{R}^d \to [0, \infty]$ is a convex lower semi-continuous function with superlinear growth at infinity. The perspective function (aka the entropy) of $F$ evaluated at the measures $(M, \nu)$ is defined as

$$ (\nu, M) \to \mathcal{F}(M|\nu) = \begin{cases} \mathbb{E}_{\nu}(F(\frac{dM}{d\nu})) & \text{if } M \ll \nu. \\ +\infty & \text{else.} \end{cases} $$

The classic example is the Shanon relative entropy $F(\cdot) = -\log(\cdot|\bar{M})$. The “nice” properties of $\mathcal{F}$ are inherited from its dual formulation (37),
**Definition 2** (Weak Fokker-Planck Solutions). We will refer to this property as: "$(\nu_t, b_t, a_t) \in FP(\rho_0, \rho_1)$".

A curve $(\nu_t)_{t \in [0,1]}$ of probability measures is a weak solution of the Fokker-Planck with drift coefficients $(b_t)$ and diffusion coefficients $(a_t)$ if $(b_t, a_t) \in L_{t,x}^1(\nu_t)$ and for all $f \in C_b^{1,2}((0,1) \times \mathbb{R}^d)$

$$ \int_0^1 \int [\partial_t f + b_t \partial_x f + a_t \partial_{xx} f] d\nu_t(x) dt = 0 $$

In our case we add the initial-final marginal conditions: $\nu_{0,1} = \rho_{0,1}$. It is not restrictive to assume that $(\nu_t)$ is narrowly continuous. For more details and $d > 1$ see [37] definition 2.2 and remark 2.3.

**Proposition 12.** i) $\mathcal{F}$ in definition 1 is jointly convex and lower semi continuous in $(M, \nu)$.

ii) $(M^h, \nu^h)$ two sequences of Borel positive measure in $\mathbb{R}^d$, such that, $\nu^h$ weakly converges to $\nu^0$, $M^h$ is a.c. w.r.t. $\nu_h$ for all $h$ and $\sup_h \mathbb{E}_{\nu^h}(F(\frac{dM^h}{d\nu^h})) < \infty$. Then $M^0$ is a.c. w.r.t. $\nu^0$ and $\liminf_{h \searrow 0} \mathbb{E}_{\nu^h}(F(\frac{dM^h}{d\nu^h})) \ge \mathbb{E}_{\nu^0}(F(\frac{dM^0}{d\nu^0}))$.

iii) Let $\eta_\epsilon$ be a smooth regularisation kernel with $k$ bounded derivatives. Denote $M_\epsilon = M * \eta_\epsilon$, $\nu_\epsilon = \nu * \eta_\epsilon$. Then, $\mathcal{F}(M_\epsilon|\nu_\epsilon) \le \mathcal{F}(M|\nu)$

iv) If $(\nu_t)$ is a weak solution of $FP_{\rho_{0,1}}(b_t, a_t)$, $(b_t, a_t) \in \mathcal{L}^p(\nu_t)$ $(\nu_{t,\epsilon})$ is a weak solution of $FP_{\rho_{0,1,\epsilon}}(b_{t,\epsilon}, a_{t,\epsilon})$ where for all $t$ $(\nu_{t,\epsilon} := (\nu_t * \eta_\epsilon)$ and we mollify the "moments":

$$ b_{t,\epsilon} = \frac{(b_t d\nu_t) * \eta_\epsilon}{d\nu_{t,\epsilon}} \qquad a_{t,\epsilon} = \frac{(a_t d\nu_t) * \eta_\epsilon}{d\nu_{t,\epsilon}}. $$

Then, $(b_{t,\epsilon}, a_{t,\epsilon})$ are in $\mathcal{L}^p(\nu_{t,\epsilon})$ and well defined uniformly bounded $C^k$ densities: $\sup_{t \in [0,1]} \|(b_{t,\epsilon}, a_{t,\epsilon})\|_{C_b^k(B)} < +\infty$ for all bounded set $B \in \mathbb{R}^d$.

v) Applying iii) to the setting in iv) we get $\mathbb{E}_{\nu_{t,\epsilon}}(F(b_{t,\epsilon})) \le \mathbb{E}_{\nu_t}(F(b_t))$.

*Proof.* These results are a direct application of the dual form of the entropy given in lemma 9.4.3 and 9.4.4 [35], see also lemma A.1 [37]. □