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

HAL
open science

# The LQR-Schrödinger Bridge

Marc Lambert

► To cite this version:

Marc Lambert. The LQR-Schrödinger Bridge. 64th IEEE Conference on Decision and Control, Dec 2025, Rio de Janeiro, France. hal-05108441v2

HAL Id: hal-05108441

https://inria.hal.science/hal-05108441v2

Submitted on 22 Nov 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The LQR-Schrödinger Bridge

Marc Lambert

INRIA - Ecole Normale Supérieure - PSL Research university

DGA - French Procurement Agency

marc.lambert@inria.fr

## Abstract

We consider the Schrödinger bridge problem in discrete time, where the pathwise cost is replaced by a sum of quadratic functions, taking the form of a pathwise linear quadratic regulator (LQR) cost. This cost comprises potential terms, which act as attractors, and kinetic terms, which control the diffusion of the process. When the two boundary marginals are Gaussian, we show that the LQR-Schrödinger bridge problem can be solved in closed form. We follow the dynamic programming principle, interpreting the Kantorovich potentials as cost-to-go functions. Under the LQR-Gaussian assumption, these potentials can be propagated exactly through backward and forward passes, leading to a system of dual Riccati equations, well known in estimation and control. This system converges rapidly in practice. We then show that the optimal process is Markovian and compute its transition kernel in closed form as well as the Gaussian marginals. Through numerical experiments, we demonstrate that this approach can be used to construct complex, non-homogeneous Gaussian processes with acceleration and loops, given well-chosen attractive potentials. Moreover, this approach allows extending the Bures transport between Gaussian distributions to more complex geometries with negative curvature.

## 1 INTRODUCTION

In his 1931 foundational article Chetrite et al. [2021], Schrödinger addressed the problem of determining the most likely process given two observations of its marginal distributions: $p_0$ at the initial time and $p_K$ at a subsequent time. Following a maximum entropy approach, he considered as a prior, that the process is a Brownian motion starting from $p_0$ and then correct this hypothesis based on the observed marginal $p_K$. At that time, probability theory had not yet been formally established, yet Schrödinger successfully modeled this problem using a system of coupled equations derived purely from physical considerations. This system of equations, known as the Schrödinger system, has attracted renewed interest in recent decades Léonard [2014], Gentil et al. [2017], Peyré and Cuturi [2019], Chen et al. [2021] due to its connection with entropy-regularized optimal transport and the Sinkhorn algorithm Sinkhorn [1964], which can be used to solve the optimal transport problem efficiently Cuturi [2013].

The Schrödinger system is based on the dual formulation of optimal transport, where two Lagrangian functions $V_0$ and $V_K$, also called Kantorovich potentials Kantorovich [1958], encode the constraints associated with the two bound marginals $p_0$ and $p_K$. The potential $V_K$ can be interpreted as a final cost in control theory, which can be propagated backward via dynamic programming to compute the so-called cost-to-go $V_k$ at intermediate times. Without an initial constraint, a single backward pass is sufficient; however, if we impose $p_0$, then a system of forward-backward equations must be solved. When the reference measure is a Brownian motion, efficient solvers such as the Sinkhorn algorithm exist. The case where the reference is a linear stochastic process is more involved
and has received significant attention Levy et al. [1990], Beghi [1996], Chen et al. [2016], Bakolas [2016], Bunne et al. [2022]. In particular, in Chen et al. [2016], the Schrödinger bridge problem is shown to be equivalent to a stochastic control problem with a quadratic cost on the control where the solution reduces to a system of forward-backward Riccati equations.

We propose an alternative formulation of the above problem: instead of specifying the reference measure through a linear stochastic process, we define it via a joint distribution determined by a pathwise LQR cost Kalman [1960]. This cost incorporate both a quadratic penalty on the control and a quadratic penalty that enforces proximity to a reference path. This setting can be useful not only to steer the covariance of the state to a final value Hotz and Skelton [1985] but also to guide it along a path or a tube Okamoto and Tsiotras [2019]. We show that the resulting problem reduces to a system of Riccati equations bearing a strong resemblance to the Kalman equations for estimation and control. These equations differ from previous works: while the forward equations in Chen et al. [2016], Bakolas [2016] describe the evolution of the state covariance matrix for the closed-loop system, our forward equations govern the evolution of the forward Kantorovich potential. Moreover, the pathwise cost introduces additional coupling in the system. Our results are derived in the discrete-time case and restricted to the setting where the velocity is directly controlled, so that the control takes the form $u_k \propto x_{k+1} - x_k$. As in the least action principle, the trajectory will then be fully determined by the cost function. This setting is also common in the mean field game approach of Schrödinger bridge Benamou et al. [2019]. Extending the framework to control through a channel for a pathwise LQR cost is left for future work.

More precisely, we define our reference measure on the path $x_0, \dots, x_K$ as a Gibbs distribution $\exp(-\frac{1}{\epsilon} \sum_{k=0}^{K-1} \ell_k(x_k, x_{k+1}))$ where $\epsilon$ is a temperature parameter and $\ell_k(x_k, x_{k+1})$ is the LQR cost transition which takes the form:

## The LQR cost

$$ \frac{1}{2}(x_k - x_k^*)^T Q_k (x_k - x_k^*) + \frac{1}{2}(x_{k+1} - x_k)^T R_k (x_{k+1} - x_k). $$

We can reinterpret this cost in the context of stochastic processes as follows: the left-hand term encodes potential energy through the penalization matrix $Q_k$ which attracts the process around the control points $x_k^*$. The right-hand term encodes kinetic energy through the penalization matrix $R_k$, which limits the diffusion of the process. While the kinetic term is consistent with the dynamic formulation of optimal transport Benamou [2000], the potential term is more unusual and will play an important role in our model. When considering Gaussian marginals and this quadratic LQR cost, the Kantorovich potential becomes quadratic and can be represented by symmetric positive definite matrices $P_0^\oplus, P_K^\ominus$ and drift terms $\alpha_0^\oplus, \alpha_K^\ominus$. We can then show that the forward and backward passes of the Schrödinger system can be computed in closed form, yielding a system of discrete algebraic Riccati equations Lancaster and Rodman [2002] on matrices:

## The forward-backward discrete Riccati equations

$$ P_k^\ominus = Q_k/\epsilon + P_{k+1}^\ominus - P_{k+1}^\ominus (R_k/\epsilon + P_{k+1}^\ominus)^{-1} P_{k+1}^\ominus \quad (1) $$

$$ P_{k+1}^{\oplus -1} = \epsilon R_k^{-1} + P_k^{\oplus -1} - P_k^{\oplus -1} (\epsilon Q_k^{-1} + P_k^{\oplus -1}) P_k^{\oplus -1}. \quad (2) $$

This closely resembles the so-called dual Riccati equations for control (backward) and estimation (forward). Using the resulting Kantorovich potentials, we obtain a non-homogeneous Gaussian Markov process governed by the time-varying weighting matrices $Q_k$ and $R_k$. This non-homogeneity
offers a versatile modeling framework. In particular, we recover the Bures transport between Gaussians when no potential is active, while activating potentials extends it to more complex geometries problems, as illustrated in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: As we increase the potential force, shown in green, the Gaussians’ Bures geodesics, shown in blue, are deformed and ultimately reduce to a single point where the information is maximal.</figcaption></figure>

To compute this solution, we proceed in several steps. First, in Section 2 we recall the dual Kantorovich formulation of the Schrödinger problem Léonard [2001], Léonard [2014] and explain its relation to the optimal transport problem. Next, in Section 3, we consider a structured reference measure with a pairwise loss. In this case, the optimal solution is Markovian, and the non-homogeneous transition kernel can be computed recursively in discrete time Jamison [1975], Beghi [1996]. We re-derive this result using a simple re-normalization procedure. Our main contribution is given in Section 4 where we consider the case in which the pairwise loss function has a LQR quadratic form and assume also the marginals are Gaussian. We show the forward and backward passes are given by the Riccati equations on covariance and equations on drifts. In Section 5, we illustrate the behavior of these optimal Gaussian processes in two dimensions for different types of LQR costs where we activate the potentials $Q_k$ at certain key times $k$ to control the path and shape of the Gaussian process.

## Notations

We use the convention used in statistical physics and probabilistic graphical models Wainwright and Jordan [2008]. To simplify the exposition, we assume that all probability measures admit density functions on the $d$-dimensional Euclidean space $\mathbb{R}^d$. This suffices, since we ultimately derive an algorithm for Gaussian distributions, though the results of Sections 2 and 3 extend to more general measures, including discrete ones. Moreover, following the graphical model approach, we construct conditional probabilities from joint probabilities in discrete time, i.e., multivariate density functions, which significantly simplifies the theory of Markov processes. Mathematically speaking, this means we assume that all measures $dp$ considered here are absolutely continuous with respect to the Lebesgue measure $dx$ and admit a density function $p(x)$ such that $dp := p(x)dx$. $p(x)$ represents the probability distribution $\mathbb{P}$ of the random variable $X$: $\mathbb{P}(X \in A) := \int_A p(x)dx$ and called without differentiation either a density function or a probability distribution. To simplify notation, we use the same convention as in probabilistic graphical models, and write $p(x)$ and $p(y)$ for distinct probability distributions, using dummy variables to differentiate them. We define a joint probability distribution as $p(x_0, \dots, x_K)$ and assume that the underlying graphical model is a linear chain, such that the index $k$ can be interpreted as time. Given this assumption, we can define the conditional probabilities from this joint probability distribution as follows: $p(x_0, \dots, x_K) = p(x_0) \prod_{k=1}^K p(x_k|x_0, \dots, x_{k-1})$ where the dependencies may simplify if the linear chain is a Markov chain. When $p$ is a normalized probability measure and $f$ not, we use the notation $p(x) \propto f(x)$ to indicate that the equality holds up to a normalization constant such that $p(x) = \frac{1}{Z}f(x)$, where
$$Z = \int f(x)dx.$$

The (negative) entropy is defined by

$$H(dp\|dx) := \int \log \frac{dp}{dx} dp = \int p(x) \log p(x) dx,$$

where $dp$ is the Radon-Nikodym derivative which corresponds to the density function $p$ with respect to the Lebesgue measure $dx$ here. For brevity, we write $H(dp\|dx) := H(p)$, omitting the reference measure in this case. When considering the relative entropy with respect to a measure $dq = q(x)dx$, which is not simply the Lebesgue measure, we use the Kullback-Leibler (KL) notation:

$$KL(p\|q) := H(dp\|dq) = \int \log \frac{p(x)}{q(x)} p(x) dx.$$

Finally, given a square-integrable function $V \in L^2(\mathbb{R}^d)$ and a probability density $p(x)$, we use $\int V(x)p(x)dx := \mathbb{E}_p[V(x)] := \int Vdp := \langle V, p \rangle$ interchangeably, with the last one emphasizing the duality product between square-integrable functions and measures. Extending these notations to the joint probability $p(x_0, x_K)$, we see that a marginal constraint $\int p(x_0, x_K)dx_0 = p_K(x_K)$ implies that $\int V(x_K)(\int p(x_0, x_K) - p_K(x_K)dx_0)dx_K = 0$ for any square-integrable functions $V$, or equivalently $\langle V, p - p_K \rangle = 0$. We can then see a marginal constraint as a linear constraint in a Hilbert space, which will be useful for constructing the dual Kantorovich problem from a Lagrangian functional. Finally, from the above definitions, we immediately deduce the fundamental relation $KL(p\|\exp(-V)) = H(p) + \langle V, p \rangle$ which is related to the duality between free energy and entropy.

# 2 Recall on Schrödinger bridge

## 2.1 Static Schrödinger bridge

Given a cost function $c: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$, the static Schrödinger bridge is given by the entropic-regularized optimal transport problem from a measure $p_0$ to a measure $p_K$:

$$\min_{p \in \mathcal{P}(p_0, p_K)} \mathbb{E}_p[c(x_0, x_K)] + \varepsilon H(p) = \min_{p \in \mathcal{P}(p_0, p_K)} \varepsilon \mathrm{KL} \left( p(x_0, x_K) \middle\| \exp(-c(x_0, x_K)/\varepsilon) \right), \quad (3)$$

where $\mathcal{P}(p_0, p_K)$ is the set of joint probability distributions (or couplings) on $(x_0, x_K)$ with marginal constraints $(p_0, p_K)$. $\varepsilon$ is a temperature parameter, as $\varepsilon \to 0$, we recover the classical optimal transport problem.

## 2.2 Dynamic Schrödinger bridge

The above problem can be generalized to the dynamic Schrödinger bridge if we replace the cost $c(x_0, x_K)$ by a loss along a path $\ell(x_0, x_1, \dots, x_{K-1}, x_K)$ introducing $K-1$ intermediate variables $x_1, \dots, x_{K-1}$:

$$\begin{align}
& \min_{p \in \mathcal{P}(p_0, p_K)} \mathbb{E}_p[\ell(x_0, \dots, x_K)] + \varepsilon H(p) \tag{4} \\
&= \min_{p \in \mathcal{P}(p_0, p_K)} \varepsilon \mathrm{KL} \left( p(x_0, \dots, x_K) \middle\| \exp(-\ell(x_0, \dots, x_K)/\varepsilon) \right) \nonumber \\
&= \min_{p \in \mathcal{P}(p_0, p_K)} \varepsilon \mathrm{KL} \left( p \middle\| \exp(-\ell/\varepsilon) \right), \nonumber
\end{align}$$
where $\mathcal{P}$ is now the set of joint probabilities on the path $x_0, \dots, x_K$ with two boundary marginal constraints $(p_0, p_K)$. The right term $\exp(-\ell/\epsilon)$ is called the reference measure and is not normalized. This lack of normalization plays a fundamental role in the duality between free energy and entropy, as we discuss in the next paragraph. If the reference measure is a Brownian motion, the static and dynamic Schrödinger bridges give the same optimal coupling Peyré and Cuturi [2019].

## 2.3 The Kantorovich dual formulation

We now introduce the Kantorovich dual formulation of the Schrödinger bridge problem, for the derivation in the case of general measures (not just density measures), please refer to Léonard [2014], Léonard [2001]. The marginal constraints can be rewritten, for all $V_0, V_K \in L^2(\mathbb{R}^d)$, as:

$$ \int V_0 dp = \int V_0 dp_0; \quad \int V_K dp = \int V_K dp_K, $$

where $V_0$ and $V_K$ will play the role of the Lagrange parameters satisfying the linear constraint $\langle V_0, p - p_0 \rangle = 0$ and $\langle V_K, p - p_K \rangle = 0$. These Lagrange parameters are closely related to the so-called dual Kantorovich potentials which play an important role in the theory of optimal transport Brenier [1991], Benamou [2000]. Since the KL is convex, there is no duality gap Léonard [2001] and we can switch the min and the max such that the Lagrangian takes the form (where $C_{0,K} := \int V_0 dp_0 + \int V_K dp_K$):

$$ \begin{aligned} & \min_{p \in \mathcal{P}(p_0, p_K)} \text{KL} \left( p \mid \exp(-\ell/\epsilon) \right) \\ &= \min_p \max_{V_0, V_K} \text{KL} \left( p \mid \exp(-\ell/\epsilon) \right) + \langle V_0, p - p_0 \rangle + \langle V_K, p - p_K \rangle \\ &= \max_{V_0, V_K} \min_p \text{KL} \left( p \mid \exp(-V_0) \exp(-\ell/\epsilon) \exp(-V_K) \right) - C_{0,K}. \end{aligned} $$

The optimal probability is given by $p^* = \frac{1}{Z^*} \exp(-V_0) \exp(-\ell/\epsilon) \exp(-V_K)$ and the optimal value is:

$$ \begin{aligned} & \max_{V_0, V_K} - \log Z^* - C_{0,K} \\ &= \max_{V_0, V_K} - \log \int \exp(-V_0) \exp(-\ell/\epsilon) \exp(-V_K) - C_{0,K}, \end{aligned} $$

which is the dual free-energy problem Léonard [2001]. Finally we rewrite the optimal solution as follows:

$$ p^*(x_0, \dots, x_K) \propto \phi_0(x_0) r(x_0, \dots, x_K) \phi_K(x_K), \quad (5) $$

where we have introduced the Gibbs potentials $\phi_0$ and $\phi_K$ supposed normalized and the unnormalized reference measure $r$ defined as follows:

$$ \begin{aligned} \phi_0(x_0) &= \frac{1}{Z_0} \exp(-V_0(x_0)) \\ \phi_K(x_K) &= \frac{1}{Z_K} \exp(-V_K(x_K)) \\ r(x_0, \dots, x_K) &= \exp(-\ell(x_0, \dots, x_K)/\epsilon). \end{aligned} $$
## 2.4 The forward-backward Schrödinger system

Reinserting the marginal constraints in the optimal solution (5), we find that the Gibbs potentials $\phi_0, \phi_K$ must satisfy the Schrödinger system (given here up to a normalization constant):

$$p_0(x_0) \propto \phi_0(x_0) \int \phi_K(x_K) r(x_0, \dots, x_K) dx_1 \dots dx_K, \qquad (6)$$

$$p_K(x_K) \propto \phi_K(x_K) \int \phi_0(x_0) r(x_0, \dots, x_K) dx_0 \dots dx_{K-1}. \qquad (7)$$

When there are no intermediate points and only the two bounds $(x_0, x_K)$, we retrieve the iterative proportional fitting (IPF) algorithm Csiszár [1975]. When, in addition, the two marginal measures are discrete, we recover the Sinkhorn algorithm Sinkhorn [1964], which can be computed efficiently by Cuturi [2013]. However, when we have a pathwise reference measure with intermediate points, the problem is more involved and requires complex computation of integrals. In this case, we need to propagate the potential $\phi_K(x_K)$ backward in equation (6) to update $\phi_0(x_0)$. Indeed, if $\phi_0^\ominus(x_0) \propto \int \phi_K(x_K) r(x_0, \dots, x_K) dx_1 \dots dx_K$ is the backward propagated potential then we can compute the initial potential as follows $\phi_0 = p_0/\phi_0^\ominus$. The updated potential $\phi_0$ can then be propagated forward in equation (7) to update the terminal potential in the same way $\phi_K = p_K/\phi_K^\oplus$. The process is repeated until, hopefully, stabilization of the potentials. Guarantees of convergence may be hard to obtain, but we will see that in the LQR case, fast convergence is observed. This complex forward-backward process simplifies significantly if the reference measure factorizes into a pairwise structure (like Markovian measure). In such cases, the propagation can be computed recursively, allowing us to determine intermediate potentials $\phi_k^\oplus$ and $\phi_k^\ominus$. Furthermore, the optimal solution is necessarily a (non-homogeneous) Markov chain and we can explicitly compute all its transition kernels. This will be the focus of the next section.

# 3 The Markovian-Schrödinger bridge

## 3.1 The pairwise reference measure

We now refine the setting by considering a loss function that splits into pairwise terms as follows:

$$\ell(x_0, \dots, x_K) = \sum_{k=0}^{K-1} \ell_k(x_k, x_{k+1}). \qquad (8)$$

**Remark 1 (Connexion with least action principle)** This structured loss can be seen as the discrete counterpart of the action integral $S(x_0) = \int_0^T L(x_t, \dot{x}_t, t)dt$ where L is the Lagrangian encoding the potential and kinetic energy of a physical system. The action integral S was already proposed as a cost function in optimal transport in the continuous-time setting; leading to the rich theory of optimal transport on manifolds Villani [2008]. From equation (4), we see also that the Schrödinger bridge problem can be interpreted as a generalization of the principle of least action, with an additive entropic term and bounded constraints.

When considering the loss (8), the unnormalized reference measure factorizes as a pairwise Markov random field Wainwright and Jordan [2008]. Let's define $r(x_{k+1}, x_k) = \exp(-\ell(x_k, x_{k+1})/\varepsilon)$, we have the relation $r(x_0, \dots, x_K) = \exp(-\ell(x_0, \dots, x_K)/\varepsilon) = \prod_{k=0}^{K-1} r(x_{k+1}, x_k)$ from which we
deduce from equation (5):

$$p^*(x_0, \dots, x_K) \propto \phi_0(x_0) \prod_{k=0}^{K-1} r(x_{k+1}, x_k) \phi_K(x_K),$$

and the optimal solution is itself a pairwise Markov random field with bounded constraint.

From this pairwise structure, we can construct a probability distribution that is a non-homogeneous forward Markov chain $p^*(x_0, \dots, x_K) = p_0(x_0) \prod_{k=0}^{K-1} p^*(x_{k+1} | x_k)$ where the transition kernel $p^*(x_{k+1} | x_k)$ can be constructed from the loss $\ell_k$ as we will show shortly.

## 3.2 The optimal Kernel transition

To compute the optimal transition kernel $p^*(x_{k+1} | x_k)$, we follow the Bellman dynamic programming approach Bellman [2010]. Starting from the final potential $\phi_K(x_K)$, we first compute the last transition $p^*(x_K | x_{K-1})$ and then deduce the form of the next propagated potential $\phi_{K-1}^{\ominus}(x_{K-1})$. We then proceed by recursion to deduce all the subsequent transitions.

Starting from $x_K$, using the pairwise structure of $r$ and the Markov property of $p^*$, we have:

$$\begin{aligned}
p^*(x_K | x_{K-1}) &= \phi_0(x_0) \frac{r(x_0, \dots, x_{K-1})}{p^*(x_0, \dots, x_{K-1})} r(x_K, x_{K-1}) \phi_K(x_K) \\
&= Cr(x_K, x_{K-1}) \phi_K(x_K),
\end{aligned}$$

where we have regrouped all the terms which do not depend on $x_K$ in $C$. Since $p^*(x_K | x_{K-1})$ is a probability distribution, we must have $\int Cr(x_K, x_{K-1}) \phi_K(x_K) dx_K = 1$ and replacing $C$ by its expression we deduce:

$$\begin{aligned}
p^*(x_0, \dots, x_{K-1}) &= \phi_0(x_0) r(x_0, \dots, x_{K-1}) \int r(x_K, x_{K-1}) \phi_K(x_K) dx_K \\
&= \phi_0(x_0) r(x_0, \dots, x_{K-1}) \phi_{K-1}^{\ominus}(x_{K-1}),
\end{aligned}$$

where we have introduced the backward-propagated Gibbs potential

$$\phi_{K-1}^{\ominus}(x_{K-1}) := \int r(x_K, x_{K-1}) \phi_K(x_K) dx_K = \frac{1}{C}.$$

We can then proceed by recursion to find all the conditionals:

$$p^*(x_{k+1} | x_k) = \frac{1}{\phi_k^{\ominus}(x_k)} r(x_{k+1}, x_k) \phi_{k+1}^{\ominus}(x_{k+1}) \quad (9)$$

$$\phi_k^{\ominus}(x_k) = \int r(x_{k+1}, x_k) \phi_{k+1}^{\ominus}(x_{k+1}) dx_{k+1}, \quad (10)$$

where $r(x_{k+1}, x_k) = \exp(-\ell(x_{k+1}, x_k))$.

**Remark 2 (Connexion with Doob's transform & RL)** Through this re-normalization procedure in discrete time, we recover a Doob's recursion Jamison [1975], Beghi [1996]. Doob's transform is used to condition a reference process to a terminal constraint and is obtained by an importance reweighting of the reference process, where, from equation (9), the importance weight is $\frac{\phi_{k+1}^{\ominus}(x_{k+1})}{\phi_k^{\ominus}(x_k)}$. The recursion (9)-(10) also emerges as a special case of the variational dynamic programming principle introduced in Lambert et al. [2024], specifically when the velocity is controllable, leading to further connections with maximum entropy control Ziebart et al. [2010] and reinforcement learning Haarnoja et al. [2018].
To completely solved the Schrödinger bridge system we need to write also the recursion for the forward pass (7). Using the pairwiswe structure of $r$ it is straightforward, and we obtain:

$$ \phi_k^{\oplus}(x_{k+1}) = \int r(x_k, x_{k+1}) \phi_k^{\oplus}(x_k) dx_k. \qquad (11) $$

We have obtained the closed-form expressions for the optimal transition kernel $p^*(x_{k+1}|x_k)$ and the optimal marginals, $p^*(x_k)$, but the computation of the marginalizations in equations (10)-(11) remains intractable in general. We now refine the setting further to find closed-form formulas, considering a quadratic loss and Gaussian marginals. In this case, the Gibbs potentials $\phi_0, \phi_K$ are also Gaussian and preserve their Gaussian structure in both the backward and forward passes. This remarkable property of stability of the quadratic-Gaussian model is at the foundation of the Kalman theory for estimation and control Kalman [1960] and is the subject of the next section, where we will consider an LQR loss.

# 4 The LQR-Schrödinger bridge

## 4.1 Recall on results on Gaussian

Before delving into the details of the subject, we recall some notations and results on operations with Gaussians in dimension $d$, which will be useful in the sequel. A Gaussian probability distribution with mean $\mu \in \mathbb{R}^d$ and covariance matrix $\Sigma \in \mathbb{R}^{d \times d}$ is denoted by $q = \mathcal{N}(\mu, \Sigma)$. The associated Gaussian density function is $q(x) = \mathcal{N}(x|\mu, \Sigma)$. Given these notations, we have the following properties:

1. The product of two Gaussian densities $\mathcal{N}(\mu_1, \Sigma_1)$ and $\mathcal{N}(\mu_2, \Sigma_2)$ is, up to a normalization constant $Z$, a Gaussian density, defined by:

$$ \mathcal{N}(\mu^*, \Sigma^*) = \frac{1}{Z} \mathcal{N}(\mu_1, \Sigma_1) \mathcal{N}(\mu_2, \Sigma_2) \qquad (12) $$

$$ \mu^* = \Sigma^* (\Sigma_1^{-1} \mu_1 + \Sigma_2^{-1} \mu_2); \quad \Sigma^* = (\Sigma_1^{-1} + \Sigma_2^{-1})^{-1}; \\ Z = \mathcal{N}(\mu_1|\mu_2, \Sigma_1 + \Sigma_2). \qquad (13) $$

2. The convolution of two Gaussian densities is a Gaussian, defined by:

$$ \int \mathcal{N}(x|Az + b, \Sigma_x) \mathcal{N}(z|\mu_z, \Sigma_z) dz = \mathcal{N}(A\mu_z + b, \Sigma_x + A\Sigma_z A^T), $$

where $A \in \mathbb{R}^{d \times d}$ is a matrix and $b \in \mathbb{R}^d$ is a vector.

We will use these results to propagate the Gaussian potentials through the LQR quadratic cost.

## 4.2 The LQR cost with Gaussian marginals

Inspired by LQR control, we consider a pairwise quadratic loss that mimics the Lagrangian, with a potential energy term and a kinetic energy term:

$$ \ell_k(x_k, x_{k+1}) := \frac{1}{2}(x_k - x_k^*)^T Q_k (x_k - x_k^*) + \frac{1}{2}(x_{k+1} - x_k)^T R_k (x_{k+1} - x_k). $$

The potential energy depends on a weighting symmetric positive definite matrix $Q_k \in d \times d$, which operate as a path-following cost that enforces the process to follow a specific trajectory $x_0^*, \dots, x_K^*$.
The kinetic energy depends on a weighting symmetric positive definite matrix $R_k \in d \times d$, which penalizes the diffusion of the process and prevents it from making too big jumps.

The pairwise refrence measure is then:

$$r(x_k, x_{k+1}) = \exp(-\ell_k(x_k, x_{k+1})/\epsilon) \propto \mathcal{N}(x_k|x_k^*, \epsilon Q_k^{-1})\mathcal{N}(x_{k+1}|x_k^*, \epsilon R_k^{-1}).$$

Moreover, we suppose that the two bounds marginals are Gaussians distributions given by:

$$p_0 = \mathcal{N}(\mu_0, \Sigma_0); \quad p_K = \mathcal{N}(\mu_K, \Sigma_K).$$

Since $p_0, p_K$ are Gaussians, it is sufficient to choose a Gaussian model for the Gibbs potentials $\phi_0, \phi_K$ to respect the marginals constraints equations (7) and (6).

We thus consider the following Gaussian structure for the unknown Gibbs potentials:

$$\phi_0 = \mathcal{N}(\alpha_0, P_0^{-1}); \quad \phi_K = \mathcal{N}(\alpha_K, P_K^{-1}),$$

associated to the Kantorovich potentials:

$$V_0(x_0) = (x_0 - \alpha_0)^T P_0 (x_0 - \alpha_0)$$

$$V_K(x_K) = (x_K - \alpha_K)^T P_K (x_K - \alpha_K).$$

Finally, at subsequent times $k$, we note the forward propagated potential $\phi_k^{\oplus} = \mathcal{N}(\alpha_k^{\oplus}, P_k^{\oplus^{-1}})$, and the backward propagated potential $\phi_k^{\ominus} = \mathcal{N}(\alpha_k^{\ominus}, P_k^{\ominus^{-1}})$.

## 4.3 The forward-backward Riccati solutions

### 4.3.1 backward potentials

The recursive update for the backward propagated potentials (10) simplifies as:

$$\phi_k^{\ominus}(x_k) = \int r(x_k, x_{k+1}) \phi_{k+1}^{\ominus}(x_{k+1}) dx_{k+1} \\ \propto \mathcal{N}(x_k|x_k^*, \epsilon Q_k^{-1}) \int \mathcal{N}(x_{k+1}|x_k, \epsilon R_k^{-1}) \mathcal{N}(x_{k+1}|\alpha_{k+1}^{\ominus}, P_{k+1}^{\ominus^{-1}}) dx_{k+1} \\ \propto \mathcal{N}(x_k|x_k^*, \epsilon Q_k^{-1}) \mathcal{N}(x_k|\alpha_{k+1}^{\ominus}, \epsilon R_k^{-1} + P_{k+1}^{\ominus^{-1}}),$$

where we have used equation (13). Using equation (12) to find the normalized solution, we finally obtain:

$$\phi_k^{\ominus}(x_k) = \mathcal{N}(x_k|\alpha_k^{\ominus}, P_k^{\ominus^{-1}}),$$

where we recover equation (1) for the covariance matrix:

$$\begin{aligned} P_k^{\ominus} &= Q_k/\epsilon + (\epsilon R_k^{-1} + P_{k+1}^{\ominus^{-1}})^{-1} \\ &= Q_k/\epsilon + P_{k+1}^{\ominus} - P_{k+1}^{\ominus}(R_k/\epsilon + P_{k+1}^{\ominus})^{-1} P_{k+1}^{\ominus}, \end{aligned}$$

which comes from the Woodbury formula Hager [1989]. This update correspond to the LQR backward Riccati equation Kalman [1960] in the special case of an identity transition matrix $A_k = I$. And for the mean we obtain:

$$\begin{aligned} \alpha_k^{\ominus} &= P_k^{\ominus^{-1}} (Q_k/\epsilon x_k^* + (\epsilon R_k^{-1} + P_{k+1}^{\ominus^{-1}})^{-1} \alpha_{k+1}^{\ominus}) \\ &= \alpha_{k+1}^{\ominus} + P_k^{\ominus^{-1}} Q_k/\epsilon (x_k^* - \alpha_{k+1}^{\ominus}), \end{aligned}$$

which reduce to $\alpha_k^{\ominus} = \alpha_{k+1}^{\ominus}$ if $Q_k = 0$ (no bias term).
### 4.3.2 forward potentials

The recursive update for the forward propagated potentials (11) simplifies as:

$$ \phi_{k+1}^{\oplus}(x_{k+1}) = \int r(x_k, x_{k+1}) \phi_k^{\oplus}(x_k) dx_k \\ \propto \int \mathcal{N}(x_k|x_k^*, \epsilon Q_k^{-1}) \mathcal{N}(x_k|\alpha_k^{\oplus}, P_k^{\oplus-1}) \mathcal{N}(x_{k+1}|x_k, \epsilon R_k^{-1}) dx_k \\ \propto \int \mathcal{N}(x_k|\alpha_{k+1}^{\oplus}, (Q_k/\epsilon + P_k^{\oplus})^{-1}) \mathcal{N}(x_{k+1}|x_k, \epsilon R_k^{-1}) dx_k, $$

where we have used equation (12) such that $\alpha_{k+1}^{\oplus} = (Q_k/\epsilon + P_k^{\oplus})^{-1}(Q_k/\epsilon x_k^* + P_k^{\oplus} \alpha_k^{\oplus})$. Using equation (12), we finally obtain:

$$ \phi_k^{\oplus}(x_{k+1}) = \mathcal{N}(x_{k+1}|\alpha_{k+1}^{\oplus}, P_{k+1}^{\oplus-1}), $$

where we recover equation (2) for the covariance matrix:

$$ P_{k+1}^{\oplus-1} = \epsilon R_k^{-1} + (Q_k/\epsilon + P_k^{\oplus})^{-1} \\ = \epsilon R_k^{-1} + P_k^{\oplus-1} - P_k^{\oplus-1}(\epsilon Q_k^{-1} + P_k^{\oplus-1})P_k^{\oplus-1}, $$

which comes again from the Woodbury formula Hager [1989]. Note that if $Q_k = 0$ we may not use this formula.

**Remark 3 (Duality of Riccati equations)** Remarkably, we find again a Riccati equation on $P^{-1}$ but where the roles of R and Q are interchanged. This give a new interpretation on the so called Kalman duality where the inverse of LQR costs matrixes R and Q can be interpreted as the diffusion matrix of the model noise and “observation” noise respectively. Indeed $Q_k^{-1}$ can be interpreted as our uncertainty (in terms of covariance) on the prior information. If $Q_k^{-1} = 0$, then we can predict perfectly where will be the trajectory at $x_k^*$, otherwise we can only estimate a neighborhood around $x_k^*$.

Finally for the drift term we obtain the recursion:

$$ \alpha_{k+1}^{\oplus} = (Q_k/\epsilon + P_k^{\oplus})^{-1}Q_k/\epsilon x_k^* + (Q_k/\epsilon + P_k^{\oplus})^{-1}P_k^{\oplus} \alpha_k^{\oplus}, $$

which reduce to $\alpha_{k+1}^{\oplus} = \alpha_k^{\oplus}$ if $Q_k = 0$ (no bias term).

### 4.3.3 update of potentials

At the end of the first backward pass, we update the initial potential $\phi_0$ as follows:

$$ \phi_0(x_0) = \frac{p_0(x_0)}{\phi_0^{\ominus}(x_0)} = \frac{\mathcal{N}(\mu_0, \Sigma_0)}{\mathcal{N}(\alpha_0^{\ominus}, P_0^{\ominus-1})}. $$

If $P_0^{\ominus} \prec \Sigma_0^{-1}$ we obtain a well defined Gaussian given by:

$$ \phi_0(x_0) = \mathcal{N}(\alpha_0^*, P_0^{*-1}) \\ P_0^* = \Sigma_0^{-1} - P_0^{\ominus} \\ \mu_0^* = (\Sigma_0^{-1} - P_0^{\ominus})^{-1}(\Sigma_0^{-1} \mu_0 - P_0^{\ominus} \alpha_0^{\ominus}). $$

We then set $\phi_0^{\oplus}(x_0) = \phi_0(x_0)$ and proceed to a forward pass to propagate the updated potential. At the end, we obtain $\phi_K^{\oplus}(x_K)$ which can be used to update the terminal potential $\phi_K(x_K) = \frac{p_K(x_K)}{\phi_K^{\oplus}(x_K)}$. We iterate several backward and forward passes until convergence. After convergence, we compute the optimal transition kernel using an ultimate backward pass.
### 4.3.4 optimal transition kernel

The recursive update for the optimal transition kernel (9) simplifies as:

$$ 
\begin{aligned}
p^*(x_{k+1}|x_k) &= \frac{1}{\phi_k(x_k)} r(x_k, x_{k+1}) \phi_{k+1}(x_{k+1}) \\
&\propto \mathcal{N}(x_{k+1}|x_k, \varepsilon R_k^{-1}) \mathcal{N}(x_{k+1}|\alpha_{k+1}, P_{k+1}^{\ominus^{-1}}) \\
&= \mathcal{N}(x_{k+1}|m_k, (R_k/\varepsilon + P_{k+1}^{\ominus})^{-1}),
\end{aligned}
$$

where we have used equation (12) to find the normalized solution and $m_k = (R_k/\varepsilon+P_{k+1}^{\ominus})^{-1}(R_k/\varepsilon.x_k+P_{k+1}^{\ominus}\alpha_{k+1})$. Finally we obtain:

$$ p^*(x_{k+1}|x_k) = \mathcal{N}(x_{k+1}|x_k + \beta_k + K_k x_k, S_k^{-1}), $$

with $S_k = R_k/\varepsilon + P_{k+1}^{\ominus}$, $K_k = S_k^{-1}R_k/\varepsilon$,

and $\beta_k = S_k^{-1}P_{k+1}^{\ominus}\alpha_{k+1}$.

$K_k$ correspond to the LQR gain when the transition matrix is the identity and $\beta_k$ is a drift term.

### 4.3.5 optimal transport geodesics

Once we have the optimal transition kernel it is straightforward to compute the stochastic trajectories. Starting from $x_0 \sim \mathcal{N}(\mu_0, \Sigma_0)$ we have the following random walk in discrete time:

$$ x_{k+1} = x_k + \beta_k + K_k x_k + \nu_k; \quad \nu_k \sim \mathcal{N}(0, S_k^{-1}). $$

Computing the mean and the covariance of $x_{k+1}$ we easily obtain the Gaussian marginals which are given by the recursion:

$$ \mathcal{N}(\mu_{k+1}, \Sigma_{k+1}) = \mathcal{N}(\beta_k + (\mathbb{I} + K_k)\mu_k, S_k^{-1} + (\mathbb{I} + K_k)^T \Sigma_k (\mathbb{I} + K_k)). $$

**Remark 4 (Connexion with Bures geodesics)** When $Q_k = 0$, the dynamic Schrödinger bridge (4) is equivalent to the static Schrödinger bridge (3) and if $\varepsilon \to 0$, we recover the classical optimal transport problem. That means these sequence of marginals define a geodesics between $\mathcal{N}(\mu_0, \Sigma_0)$ and $\mathcal{N}(\mu_K, \Sigma_K)$ which matches the Bures-geodesic in this case. When we increase the temperature $\varepsilon$ we obtain the relaxed Bures geodesic Bunne et al. [2022]. More interesting, when $Q_k \succ 0$, we obtain a richer family of geodesics living on a manifold curved by the potential $Q_k$.

# 5 Numerical Results

We compute now the solution to the LQR-Schrodinger bridge problem in two dimensions ($d = 2$) for several types of LQR costs and temperature parameter $\varepsilon$.

## 5.1 Initialization

We provide below the setting used in our simulation and recall the parameters that need to be initialized:

* The Gaussian marginals mean and covariance $\mathcal{N}(\mu_0, \Sigma_0)$ and $\mathcal{N}(\mu_K, \Sigma_K)$ are shown as red ellipsoids in the plotted figure. The potentials $\phi_0$ and $\phi_K$ are initialized as $\mathcal{N}(\mu_0, \Sigma_0^{-1})$ and $\mathcal{N}(\mu_K, \Sigma_K^{-1})$.
* The number of points $K + 1$ of the path $x_0, \dots, x_K$ is fixed between 15 and 100 depending on whether it helps visualization. Note that we use an unscaled discretization step (i.e., $dt = 1$), so the value of the loss function depends on the number of points.

* For the temperature parameter $\varepsilon > 0$ we consider two regimes: the cold regime with $\varepsilon = 0.001$ and the warm regime with $\varepsilon = 1$. The number of passes is adapted based on the temperature regime. It is well known that the Sinkhorn algorithm converges slower at low temperatures, and we observe the same behavior in the case of LQR. At the warm temperature or when $K = 100$, we observe convergence in just one pass. At the low temperature and for a small number of points $K = 15$, we use up to 5 passes.

* The diffusion matrix $R_k$ in the LQR cost is assumed fixed to a constant diagonal value $R_k = rI$ with $r > 0$.

* The potential matrix is, in most cases, also constant over time and fixed to a constant diagonal value $Q_k = qI$ with $q > 0$. The control points $x_k^*$ are chosen to steer the distribution toward arbitrary locations. We show in green the 2D ellipsoid, which corresponds to the potential matrix $Q_k$ centered at $x_k^*$. When no green ellipsoid appears in the plotted figure, there are no control points and $Q_k = 0$.

## 5.2 Optimal transport with the LQR geometry

We first illustrate in Figure 2 how the LQR cost can be used to modify the curvature of the optimal transport problem by adjusting the matrices $R$ and $Q$. We consider the optimal transport between two Diracs in the warm regime and between two Gaussians in the cold regime. In the first case, we obtain a Brownian bridge where the paths are stochastic processes (see left plot in Figure 2). If we decrease the temperature to 0, these processes shrink into deterministic straight lines, corresponding to the Euclidean geodesic. In the second case, we replace the Diracs by two Gaussians and consider a low temperature. The optimal transport is given by the Bures transport map, and the samples should form nearly straight lines at this temperature. However, when an attractive potential is added to the path, the resulting geodesics differ from the Bures transport: the geodesics follow a negative curvature (see right plot in Figure 2).

<figure><img src="image_4.png" /><figcaption>Figure 2: The figure on the left shows a standard Brownian bridge between two Dirac marginals, which is obtained when we consider a warm temperature $\varepsilon = 1$. There is no potential in the LQR cost, such that $q = 0$ and $r = 100$. The figure on the right shows an optimal transport between two Gaussians, obtained with a cold temperature $\varepsilon = 0.001$. We have added a constant potential in the middle of the path, such that $q = 0.3$ and $r = 10$. The covariances of the Gaussian process are shown in blue at $3\sigma$ here, whereas the sampled trajectories are shown in black.</figcaption></figure>
## 5.3 Application in stochastic path following and Gaussian steering

Figures 3 and 4 illustrate how the LQR framework can be used to construct complex stochastic processes by modifying the positions of the waypoints or the eigenvalues of the potential covariance matrix.

<figure><img src="image_5.png" /><figcaption>Figure 3: Temperature parameter $\varepsilon = 0.001$. From left to right: **Wave** path following with $r = 1, q = 10, K = 15$, where we also use $K$ waypoints (green ellipsoid has been rescaled for better visualization in this case); **Zig-zag** case with $r = 10, q = 0.1, K = 100$, where we use 2 waypoints to steer the Gaussian process; **Scoubidou** case with $r = 10, q = 0.2, K = 100$, where we use 3 waypoints; **Twister** obstacle with $r = 10, q = 0.2, K = 200$, where we use one waypoint with a non-isotropic potential matrix, such that the Gaussian covariances of the marginals are twisted near the potential. The blue ellipsoids show the Gaussian marginals and the green one the potential forces. All the covariances are shown at $1\sigma$ here.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Same conditions as in Figure 3, but with the temperature parameter $\varepsilon = 1$. The process is no longer deterministic in the **Wave** path following and is much more diffusive in the other experiments.</figcaption></figure>

This allows the creation of processes that follow a specific path (e.g., the “Wave” case), avoid obstacles (e.g., the “Zig-zag” case), form loops (e.g., the “Scoubidou” case), or accelerate at certain points (e.g., the “Twister” case). In the cold regime (Figure 3), the process exhibits little diffusion, and the covariance of the Gaussian decreases due to the potential’s influence, in particular in the “Wave” experiment where the process tends to become deterministic due to the high potential value. In the warm regime, the Gaussian process becomes more diffusive and less constrained by the waypoints, as shown in Figure 4. In all these cases, the potential stabilizes after just one pass.

The sources of the code are available on Github on the following repository:

https://github.com/marc-h-lambert/LQR-SchrodingerBridge.

# 6 CONCLUSION

We have integrated ideas from LQR control into the optimal transport framework, leading to a highly effective approach with some surprising implications. Initially, we derive exact closed-form updates
for the propagation of the potentials, which take the form of backward and forward Riccati equations
that are dual to each other. These equations are inherently discrete-time and can be implemented
without approximation, enabling fast convergence of the potentials in practice. Additionally, we
compute the optimal transition kernel in closed form, which allows us to generate the associated
stochastic process and its marginals. Finally, we observe that a rich geometry emerges when a
potential is incorporated into the cost function: by strategically placing attractors, we can control
Gaussian distributions, deform them, or reduce their covariance. A natural next step, particularly
for control applications, would be to integrate prior information about the process, such as a known
passive dynamics, and to restrict the optimization to the control policy $p(u|x)$, linking this work to
our earlier research in variational dynamic programming.

# 7 ACKNOWLEDGMENTS

The author gratefully acknowledges C. Léonard for insightful discussions on the dual Kantorovich formulation of the Schrödinger bridge problem, Francis Bach for constructive feedback on the draft that helped emphasize the main ideas, and Amir Taghvaei for the invitation to present this work at the CDC conference. The author also thanks the anonymous reviewer for suggesting relevant references. This work was funded by the French Defence Procurement Agency (DGA) and by the French government under the management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

# References

Efstathios Bakolas. Optimal covariance control for stochastic linear systems subject to integral quadratic state constraints. In *2016 American Control Conference (ACC)*, pages 7231–7236, 2016.

A. Beghi. On the relative entropy of discrete-time markov processes with given end-point densities. *IEEE Transactions on Information Theory*, 42(5):1529–1535, 1996.

Richard Bellman. *Dynamic Programming*. Princeton University Press, 2010. ISBN 0691146683.

Jean-David Benamou, Guillaume Carlier, Simone Di Marino, and Luca Nenna. An entropy minimization approach to second-order variational mean-field games. *Mathematical Models and Methods in Applied Sciences*, 29(08):1553–1583, 2019.

Y Benamou, J.-D.; Brenier. A computational fluid mechanics solution to the monge-kantorovich mass transfer problem. *Numer. Math.*, 84:375–393, 2000.

Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. *Communications on Pure and Applied Mathematics*, 44:375–417, 1991.

Charlotte Bunne, Ya-Ping Hsieh, Marco Cuturi, and Andreas Krause. The schrödinger bridge between gaussian measures has a closed form. In *International Conference on Artificial Intelligence and Statistics*, 2022.

Yongxin Chen, Tryphon T. Georgiou, and Michele Pavon. Optimal steering of a linear stochastic system to a final probability distribution, part i. *IEEE Transactions on Automatic Control*, 61: 1158–1169, 2016.
Yongxin Chen, Tryphon T. Georgiou, and Michele Pavon. Stochastic control liaisons: Richard sinkhorn meets gaspard monge on a schrödinger bridge. *SIAM Review*, 63(2):249–313, 2021.

Raphaël Chetrite, Paolo Muratore-Ginanneschi, and Kay Schwieger. E. schrödinger’s 1931 paper “on the reversal of the laws of nature” [“Über die umkehrung der naturgesetze”, sitzungsberichte der preussischen akademie der wissenschaften, physikalisch-mathematische klasse, 8 n9 144–153]. *The European Physical Journal H*, 46, 2021. ISSN 2102-6467.

I. Csiszár. I-divergence geometry of probability distributions as minimization problems. *Annals of Probability*, 3:146–158, 1975.

Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In *Advances in Neural Information Processing Systems*, volume 26, 2013.

Ivan Gentil, Christian Léonard, and Luigia Ripani. About the analogy between optimal transport and minimal entropy. *Annales de la Faculté des sciences de Toulouse : Mathématiques*, 26(3): 569–600, 2017.

Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine. Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor. *35th International Conference on Machine Learning*, 2018.

William W. Hager. Updating the inverse of a matrix. *SIAM Rev.*, 31:221–239, 1989.

Anthony F. Hotz and Robert E. Skelton. A covariance control theory. In *1985 24th IEEE Conference on Decision and Control*, pages 552–557, 1985.

Benton Jamison. The markov processes of schrödinger. *Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete*, 32(4):323–331, 1975.

R.E. Kalman. Contributions to the theory of optimal control. *Matematica*, 5:102–119, 1960.

L. V. Kantorovich. On the translocation of masses. *Management Science*, 5:1–4, 1958.

Marc Lambert, Francis Bach, and Silvere Bonnabel. Variational dynamic programming for stochastic optimal control. In *2024 IEEE 63rd Conference on Decision and Control (CDC)*, pages 1863–1870, 2024.

Peter Lancaster and Leiba Rodman. *Algebraic Riccati Equations*. 01 2002. ISBN 978-0-19-853795-3.

C. Léonard. Minimization of energy functionals applied to some inverse problems. *Appl. Math. Optim.*, 44(3):273–297, 2001.

B.C. Levy, R. Frezza, and A.J. Krener. Modeling and estimation of discrete-time gaussian reciprocal processes. *IEEE Transactions on Automatic Control*, 35(9):1013–1023, 1990.

Christian Léonard. A survey of the schrödinger problem and some of its connections with optimal transport. *Discrete and Continuous Dynamical Systems*, 34:1533–1574, 2014. ISSN 1078-0947.

Kazuhide Okamoto and Panagiotis Tsiotras. Optimal stochastic vehicle path planning using covariance steering. *IEEE Robotics and Automation Letters*, 4(3):2276–2281, 2019.

Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. *Foundations and Trends® in Machine Learning*, 11(5-6):355–607, 2019.
Richard Sinkhorn. A Relationship Between Arbitrary Positive Matrices and Doubly Stochastic Matrices. *The Annals of Mathematical Statistics*, 35:876 – 879, 1964.

C. Villani. *Optimal Transport: Old and New*. Springer Berlin Heidelberg, 2008.

Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. *Foundations and Trends in Machine Learning*, 1(1–2):1–305, 2008.

Brian D. Ziebart, J. Andrew Bagnell, and Anind K. Dey. Modeling interaction via the principle of maximum causal entropy. *International Conference on Machine Learning*, pages 1255–1262, 2010.