Dynamic Optimal Transport

Optimal transport becomes especially powerful once distances between measures are seen as actions of moving mass. This chapter first develops the dynamic language: continuity equations describe admissible measure evolutions, while the Benamou--Brenier formula identifies $\Wass_2$ with a least-action principle. The path-space Schrodinger problem then provides its stochastic, entropy-regularized counterpart. After extending dynamic actions to other transport geometries, the chapter closes with variational mean field games, where congestion and a terminal cost turn the Benamou--Brenier action into a population-planning problem. These ideas prepare the gradient-flow and generative-model chapters that follow.

from pathlib import Path
import sys

from IPython.display import Image as DisplayImage
from IPython.display import display

here = Path.cwd()
myst_dir = None
for candidate in [here, here.parent, here / "myst", here.parent / "myst", here.parent.parent / "myst"]:
    if (candidate / "ot4ml_web.py").exists():
        myst_dir = candidate.resolve()
        sys.path.insert(0, str(myst_dir))
        break

if myst_dir is None:
    raise RuntimeError("Could not locate myst/ot4ml_web.py")

repo_root = myst_dir.parent
thumbnails = repo_root / "notebooks-figures" / "thumbnails"

def show_book_figure(name, width=760):
    display(DisplayImage(filename=str(thumbnails / f"{name}.png"), width=width))

Evolutions Over the Space of Measures¶

We start with the continuity equation because it is the common language for particles, densities and weak measure evolutions. It also makes precise which velocity fields actually move mass.

Lagrangian and Eulerian Descriptions¶

Consider an evolution $t\mapsto\alpha_t\in\mathcal P(\RR^d)$. It can be described in a Lagrangian way as the advection of particles along a time-dependent vector field $v_t(x)$:

Writing $T_t$ for the associated flow map, so that $T_t(x(0))=x(t)$, the advected measure is

For empirical measures, $\alpha_t=n^{-1}\sum_{i=1}^n\delta_{x_i(t)}$, each particle solves (1).

In the Eulerian description, the same motion is written directly on the evolving measure:

This PDE is often called the advection equation, the continuity equation, or Liouville’s equation when it acts on phase space. It is a classical PDE only when $\alpha_t$ has a smooth density. For general measures, and in particular for empirical measures, it is understood in the distributional sense: for every $\varphi\in C_c^1((0,1)\times\RR^d)$,

This weak equation is obtained from (3) by integration by parts. For smooth positive densities, the classical and weak formulations are equivalent; for particle clouds, the weak form remains meaningful.

The interior identity does not record endpoint traces or boundary conditions. The boundary-aware version supplies the admissible class used throughout the chapter.

From Measure Evolutions to Vector Fields¶

For a given evolution $(\alpha_t)_t$, there are typically infinitely many velocity fields $v_t$ satisfying

This non-uniqueness comes from the kernel of the weighted divergence. The linear space of vector fields that leave a measure $\alpha$ invariant is

It is usually non-trivial: if $\alpha$ is an isotropic Gaussian, $\mathcal H_\alpha$ contains rotational vector fields generated by anti-symmetric matrices.

Dacorogna--Moser Inversion¶

Reconstructing particles from an observed density evolution is therefore ill-posed. For a smooth positive density $\alpha_t=\rho_t\,\d x$, a simple choice, introduced by Dacorogna and Moser Dacorogna & Moser, 1990, imposes that the flux $\rho_t v_t$ is a gradient field. With a fixed convention for the inverse Laplacian,

with suitable boundary conditions, for instance vanishing at infinity. This formula is useful conceptually but delicate when $\rho_t$ vanishes, and it does not generally produce a gradient velocity field.

The classical Dacorogna--Moser construction uses the linear density path. If $\alpha_i=\rho_i\,\d x$ are smooth positive densities with the same total mass on a bounded connected domain $\Omega$, set

Choose a time-independent flux $w$ satisfying

for instance $w=-\nabla\phi$ with $\Delta\phi=\rho_1-\rho_0$ and Neumann boundary condition. Then

satisfies $\partial_t\rho_t+\operatorname{div}(\rho_t v_t)=0$. The flow $\partial_t T_t=v_t\circ T_t$, $T_0=\operatorname{Id}$, therefore transports $\rho_0\d x$ onto $\rho_t\d x$, and $T_1$ solves the prescribed-Jacobian problem $\rho_1(T_1(x))\det(\nabla T_1(x))=\rho_0(x)$.

Least-Square Inversion and Gradient Structure¶

A more robust choice, used implicitly in flow matching, optimal transport and Wasserstein gradient flows, is to select among all admissible velocities the one with smallest kinetic energy:

In general this inversion is still computationally demanding, but special choices of $(\alpha_t)_t$ lead to simpler formulas; this is the mechanism exploited later by flow matching in Section Generative Models via Flow Matching.

Benamou--Brenier Dynamic Formulation of OT¶

The dynamic formulation identifies $\Wass_2$ with the kinetic energy of the cheapest continuity-equation path. It is the point where OT becomes a least-action principle.

Benamou--Brenier Formulation¶

Instead of assuming that a whole curve $(\alpha_t)_{t\in[0,1]}$ is prescribed, one fixes only the endpoints $\alpha_0$ and $\alpha_1$ and minimizes the least-square energy (17). The theorem of Benamou and Brenier states that this geodesic energy is exactly the squared Wasserstein distance Benamou & Brenier, 2000.

Convex Momentum-Based Reformulation¶

Although (22) is not jointly convex in $(\alpha_t,v_t)$, it becomes convex after replacing velocities by momenta. Given $v\in L^2(\alpha;\RR^d)$, define the momentum

which is a finite $\RR^d$-valued measure. The nonlinear relation $\omega=\alpha v$ is eliminated by the quadratic perspective

This lower-semicontinuous convex function is positively 1-homogeneous: $J(\eta a,\eta m)=\eta J(a,m)$ for $\eta\geq0$. If $\lambda$ is any positive measure dominating both $\alpha$ and the total variation $|\omega|$, set

The value is independent of the dominating measure: both Radon--Nikodym densities change by the same factor, and the 1-homogeneity of $J$ cancels the change of reference measure. This is the integral functional associated with a convex normal integrand in the measure-valued relaxation of dynamic OT Ambrosio et al., 2006; see also the perspective construction in Rockafellar, 2015. Moreover,

The Benamou--Brenier problem therefore has the convex measure formulation

In the absolutely continuous case $\alpha_t=\rho_t\,\d x$ and $\omega_t=m_t\,\d x$, this reduces to the familiar integral of $J(\rho_t,m_t)=\norm{m_t}^2/\rho_t$, with the zero-density conventions already encoded in (27). This convex reformulation enables geodesic interpolation by convex optimization after discretization.

Dual Hamilton--Jacobi Formulation¶

The momentum formulation also has a useful dual. It turns the least-action problem into a Hamilton--Jacobi subsolution inequality for a scalar potential, with equality on the part of space-time actually visited by the optimal curve. With the no-$1/2$ convention of (30), the constants are as follows.

This also recovers the static Kantorovich inequality from a dynamic principle. If $\gamma$ is any smooth curve with $\gamma(0)=x$ and $\gamma(1)=y$, then

After integration and minimization over curves,

Thus $(-\phi_0,\phi_1)$ is a feasible static Kantorovich dual pair for the quadratic cost. At optimality the inequality is saturated on the endpoint pairs connected by the primal characteristics.

Figure Div displays these primal--dual relations for a one-dimensional mixture transport, including the Hamilton--Jacobi contact identity along the active mass.

show_book_figure("dynamic-benamou-brenier-duality")

One-dimensional Benamou--Brenier primal and dual solutions. The endpoints are Gaussian mixtures and the solution is computed from monotone quantile interpolation. The panels show the primal density, the momentum $m_t=\rho_t v_t$, and the dual Hamilton--Jacobi potential. Along the active transported mass, the notebook checks $m_t=\rho_t\partial_x\phi_t/2$ and $\partial_t\phi_t+|\partial_x\phi_t|^2/4=0$.

Proximal Splitting¶

The convex momentum formulation also explains the original Benamou--Brenier solver. After discretization, the ALG2 scheme can be read as a Douglas--Rachford splitting, equivalently ADMM on the Fenchel--Rockafellar dual Papadakis et al., 2014. Suppressing discretization indices, write $U=(\rho,m)$, let $\mathcal F(U)$ be the integral of the perspective action, and let $\mathcal G=\iota_{\mathcal C}$ be the indicator of the affine continuity constraint with prescribed endpoints. The problem is $\min_U \mathcal F(U)+\mathcal G(U)$.

The two proximal operators separate the nonlinear and linear parts: the prox of $\mathcal F$ is local in $(t,x)$ and amounts to the perspective proximal operator, whereas the prox of $\mathcal G$ is the orthogonal projection onto the divergence equation and endpoint constraints. Douglas--Rachford alternates these two simple operations.

Figure Div complements the Eulerian optimization viewpoint with the Lagrangian picture: matched particles travel along the straight characteristics of the minimizing curve.

show_book_figure("dynamic-benamou-brenier-geodesic")

Benamou--Brenier geodesic between two sampled silhouettes. A discrete quadratic OT plan between finely subsampled cat and two-disks point clouds induces the McCann interpolation $Z_t=(1-t)X+tY$, which is the Lagrangian realization of the least-action solution. The left panel renders local color images of the smaller-bandwidth kernel-smoothed densities with enough padding to include the full silhouettes. The right panel overlays shortened velocity arrows centered at evenly subsampled midpoint particles $Z_{1/2}$; each displayed arrow runs in data coordinates from a source-side tail to a target-side head along the matched characteristic direction $Y-X$, but is not drawn as the full endpoint segment from $X$ to $Y$.

The interactive demo keeps the same Lagrangian picture: particles are matched once, then move along straight characteristics. The time and velocity scale controls separate the path $\alpha_t$ from the underlying displacement field.

Path-Space Formulation¶

Let $\Ss=C([0,1];\RR^d)$ be the space of continuous paths endowed with the uniform topology. For $t\in[0,1]$ define the evaluation map

The Benamou--Brenier cost admits the equivalent formulation

The inner energy is understood as $+\infty$ outside the absolutely continuous paths. If $\alpha_0$ has a density, the minimizer $M^*$ is unique. Its time marginals reproduce the optimal curve: $\alpha_t=(e_t)_\sharp M^*$ for all $t$. Furthermore, for a.e. $t$, the conditional law of the path velocity is deterministic:

where $v_t^*$ is the optimal velocity field in the Benamou--Brenier formulation. Hence $M^*$ concentrates on straight-line geodesics and, for a.e. $t$, assigns exactly one direction at $\alpha_t$-a.e. spatial point.

Path-Space Schrödinger Problem¶

The path-space formulation above fills each endpoint pair by a deterministic least-action path. Schrödinger’s reciprocal problem replaces this path by the conditional trajectory of a noisy reference dynamics. Optimizing the conditional path laws leaves an entropic problem over endpoint couplings, which connects dynamic transport with the Sinkhorn problem of Chapter Paragraph.

From Least-Action Paths to Random Bridges¶

Let $\X$ be a Polish state space equipped with a compatible complete bounded metric and equip $\Om=C([0,1];\X)$ with the uniform metric. Then $\Om$ is Polish, the evaluations $e_t(\omega)=\omega_t$ are continuous, and regular conditional path laws exist.

For quadratic Euclidean transport, $\mathcal{A}(\omega)=\int_0^1\norm{\dot\omega_t}^2\d t$ on absolutely continuous paths and $+\infty$ otherwise. The induced endpoint cost is

It is lower semianalytic, hence universally measurable. In the quadratic Euclidean case, $c_{\mathcal{A}}(x,y)=\norm{x-y}^2$.

Entropic Path-Space Problem¶

Let $\mathsf R^\epsilon\in\Pp(\Om)$ be a Brownian, Langevin, or other reference path law at noise level $\epsilon$.

This is Schrödinger’s entropy projection of a prior dynamics onto prescribed endpoint marginals Schrödinger, 1931Léonard, 2014. For the reference $\d X_t=\sqrt\epsilon\,\d B_t$ started from $\alpha$, Girsanov’s formula gives, under its usual integrability assumptions,

when $M$ has controlled drift $u_t$. Thus the bridge is the least energetic drift change steering the Brownian prior from $\alpha$ to $\beta$.

A zero-noise conclusion requires a path-space large-deviation theorem, not only this definition. If $\mathsf R^\epsilon$ satisfies an LDP with speed $1/\epsilon$ and good action $\mathcal{A}$, and exponential tightness plus the endpoint constraints yield constrained $\Gamma$-convergence, then minima and suitably precompact minimizers converge to the unregularized path problem Léonard, 2012. For $\d X_t=\sqrt\epsilon\,\d B_t$, the rate action is $\frac12\int_0^1\norm{\dot\omega_t}^2\d t$, so the limiting endpoint cost is $\norm{x-y}^2/2$.

Brownian Bridges and Sinkhorn Couplings¶

For $\d X_t=\sqrt\epsilon\,\d B_t$, the unit-time Brownian transition density is

Thus $\epsilon\KL$ produces the cost $\norm{x-y}^2/2$; the usual quadratic Sinkhorn convention $e^{-\norm{x-y}^2/\epsilon}$ corresponds to Brownian noise variance $\epsilon/2$. More generally, suppose that for reference probability measures $\bar\alpha,\bar\beta$ the endpoint prior is, with $0<Z_\epsilon<+\infty$,

with $\alpha\ll\bar\alpha$ and $\beta\ll\bar\beta$. For every feasible coupling, the chain rule gives

The last three terms are fixed by the marginals, so (45) is the continuous Sinkhorn problem up to an additive constant. If Brownian motion starts from $\alpha$ and $\beta=b\,\d y$, the same identity contains the fixed one-body term $\epsilon\int\log b\,\d\beta$. Thus the needed domination is $\beta\ll\d y$ (and $\alpha\ll(\mathsf R^\epsilon)_0$ for another initial reference), not an absolute-continuity relation between $\alpha$ and $\beta$.

show_book_figure("sinkhorn-path-space-bridges")

Schematic endpoint couplings lifted to Brownian bridges. The discrete picture is literal for a reciprocal reference obtained by mixing endpoint-conditioned Brownian bridges with a discrete endpoint prior; an ordinary unconditioned Brownian reference cannot have an atomic terminal law at finite entropy.

Generalized Dynamic Wasserstein Distances¶

The quadratic Benamou--Brenier formula is only one instance of a broader fixed-mass dynamic language. The goal of this section is to define a large family of geodesic-like geometries on spaces of probability measures by modifying the action minimized in the Benamou--Brenier formula. An arbitrary action first defines only a path value; metric properties require the symmetry, homogeneity, nondegeneracy, and closure assumptions isolated below. All descent constructions are postponed to Generalized Dynamic Wasserstein Flows, where these distances are used to generate gradient-flow PDE models.

Path Actions¶

The common construction replaces the quadratic kinetic energy in the Benamou--Brenier formula by an instantaneous action while retaining the continuity equation and endpoint constraints.

In the mass-preserving Euclidean setting, the basic input is an instantaneous action $\mathbb A(\alpha,w)$, where $\alpha$ is the current measure and $w$ is an admissible velocity representative.

Equivalently, one may quotient by velocity fields that induce the same first-order variation of the measure. This value need not be symmetric or separate measures, and its square root need not satisfy the triangle inequality. For an $r$-homogeneous action satisfying Proposition Proposition: Homogeneous Dynamic Actions Define Distances, its $r$-th root is a distance. Some standard distances, such as $\Wass_p$, are first written with a $p$-homogeneous action and then squared by taking a constant-speed parametrization; this normalization is made explicit below. Different choices of $\mathbb A$ change the resulting geometry; Generalized Dynamic Wasserstein Flows later reuses these choices when dynamics are introduced.

Quadratic, or Riemannian, Tangent Actions¶

A particularly transparent case occurs when $w\mapsto\mathbb A(\alpha,w)$ is quadratic. For simplicity, take admissible velocities in $L^2(\alpha;\RR^d)$; in some applications this Hilbert space is replaced by a closed subspace encoding additional constraints. Suppose the polarization of $\mathbb A$ is represented by a positive self-adjoint operator $Q_\alpha:L^2(\alpha;\RR^d)\to L^2(\alpha;\RR^d)$,

To obtain a genuine tangent norm, this quadratic form must be nondegenerate after quotienting velocity fields that induce the same measure variation.

The associated quadratic path value is

If the dynamic problem is also sequentially closed and attains finite infima as in Proposition: Homogeneous Dynamic Actions Define Distances, $\mathsf D_Q$ is an extended distance. The usual $\Wass_2$ geometry corresponds to $Q_\alpha=\Id$ in this simplified notation. Thus $Q_\alpha$ records how the chosen geometry deforms the Euclidean $L^2(\alpha)$ tangent norm: no deformation for $\Wass_2$, and a nontrivial tensor for generalized Riemannian geometries. Generalized Dynamic Wasserstein Flows later reuses the same tensor as a preconditioner for metric descent.

Local Velocity Actions¶

Many dynamic distances are local with respect to a reference measure $\lambda$. Write $\alpha=a\lambda$. A velocity action is specified by a pointwise integrand

where $a\in\RR_+$ is a density value and $w\in\RR^d$ is a velocity value, and defines

For a fixed reference $\lambda$, this covers density-dependent mobilities and congestion constraints. If, in addition, $A$ is positively 1-homogeneous in its first variable, $A(\eta a,w)=\eta A(a,w)$ for $\eta\geq0$, then the same formula is intrinsic: replacing $\lambda$ by another dominating measure gives the same value. The usual Benamou--Brenier action is the model case

Homogeneous Momentum Actions¶

The same action can be written in momentum variables, and this is the form in which convexity and metric properties are easiest to read. Set $\omega=\alpha w$, so that $\omega$ is a vector-valued measure. When the local description is written with the same reference $\lambda$, so that $\alpha=a\lambda$ and $\omega=m\lambda$, the pointwise momentum perspective is

and the measure action relative to $\lambda$ is

with value $+\infty$ if $\alpha$ or the total variation $|\omega|$ is not absolutely continuous with respect to $\lambda$. This zero-density convention is the lower-semicontinuous one for the superlinear actions used below; other growths use the corresponding recession extension. If $A$ is positively 1-homogeneous in $a$, then $J_A$ is jointly 1-homogeneous: $J_A(\eta a,\eta m)=\eta J_A(a,m)$. In that intrinsic case the value of $\mathbb J_{A,\lambda}$ is independent of the dominating reference measure, and we write simply $\mathbb J_A$. For $A_2(a,w)=a\norm{w}^2$, one recovers the quadratic perspective

which is the integrand used in the convex Benamou--Brenier formulation.

The next proposition isolates the assumptions under which the momentum formulation generated by $A$ defines a path metric rather than only a variational principle.

Without homogeneity or nondegeneracy, the same momentum action remains useful as a variational principle, but its $r$-th root need not define a distance.

Concave-Mobility Actions¶

One can instead keep a quadratic momentum action and change the mobility. Dolbeault, Nazaret and Savaré introduced this construction as a class of generalized transport distances adapted to nonlinear diffusion Dolbeault et al., 2009.

The convexity of $J_\theta$ is the special case $L(u)=\norm u^2$ of Proposition Proposition: Concave Mobilities Give Convex Momentum Actions. This is why concavity of the mobility, rather than convexity, is the structural condition that makes the continuity-equation formulation convex.

Fix now a reference measure $\lambda$. If $\alpha=a\lambda$, the induced squared tangent action is

with the closed convention above. It is $+\infty$ when $\alpha\not\ll\lambda$ or $a$ leaves $I$ on a set of positive $\lambda$-measure. Equivalently, on the set where $\theta(a(x))>0$,

Hence this is a local Riemannian case whenever the multiplier $a(x)/\theta(a(x))$ is finite and positive, in the sense of (53). The associated tensor is the multiplication operator

defined $\alpha$-a.e. on the set where $\theta(a)>0$. Except for linear mobilities $\theta(a)=ca$, and in particular the normalized case $\theta(a)=a$ which recovers $\Wass_2$, the pointwise velocity action $A_\theta(a,w)$ is not positively 1-homogeneous in the density variable $a$. Consequently the construction is not intrinsic under a change of $\lambda$: the resulting distance depends on the chosen reference measure and is finite only between endpoints that can be joined by a finite-action curve with $\alpha_t\ll\lambda$.

The subscript $\lambda$ recalls that the action is measured through the density $a=\d\alpha/\d\lambda$. Equivalently, because this action is quadratic in the momentum, $\mathsf W_{\theta,\lambda}^2$ is the path value (51) with $\mathbb A=\mathbb A_{\theta,\lambda}$.

The choice $\theta(a)=a$ recovers $\Wass_2$. Other choices encode different geometry: $\theta(a)=a^\gamma$ with $0<\gamma\leq1$ changes the cost of moving dilute mass, while $\theta(a)=a(1-a/M)$ on $[0,M]$ models a volume-filling or exclusion effect. The distance is comparable with $\Wass_2$ on classes where $\theta(a)$ is bounded above and below by positive multiples of $a$; otherwise zero-mobility barriers can make some pairs infinitely far apart.

Dynamic Spectral Wasserstein Distances¶

The static spectral distances of Spectral and Robust Wasserstein Distances penalize a coupling through the covariance of its displacement. A dynamic version keeps the continuity equation but replaces the pointwise kinetic energy by a gauge of the whole velocity covariance. The resulting action is nonlocal in space: velocity directions are charged globally through their covariance, rather than independently at each point.

The trace gauge gives the usual Wasserstein tangent action, while the operator gauge $\gamma(M)=\lambda_{\max}(M)$ charges only the largest directional velocity variance.

In density--momentum variables, this corresponds to the measure action

or, when $\alpha=\rho\,\d x$ and $\omega=m\,\d x$,

This functional is convex in the density--momentum fields $(\rho,m)$ by the matrix perspective, together with the monotonicity and convexity of $\gamma$. It is nevertheless not, in general, obtained by integrating a pointwise action density, because the covariance is computed globally before applying $\gamma$. Among spectral gauges, the linear ones are exactly the scaled traces $\gamma(M)=c\operatorname{tr}(M)$ with $c>0$. Their velocity and momentum densities are

The case $c=1$ recovers Benamou--Brenier. A functional $M\mapsto\operatorname{tr}(GM)$ with nonscalar $G\succeq0$ is Loewner-monotone but not orthogonally invariant, hence anisotropic rather than spectral.

The following result, in the form used for normalized spectral flows in Peyré, 2026, shows that this dynamic construction is not merely infinitesimal: it exactly recovers the static displacement-covariance formulation.

The use of this geometry for normalized flows, including the operator-gauge connection with Muon-type normalization, is developed in Dynamic Spectral Wasserstein Flows.

Kernelized Benamou--Brenier Distances¶

A different way to deform the Benamou--Brenier geometry is to keep the local continuity equation but to measure velocities in a reproducing-kernel Hilbert space rather than in $L^2(\alpha)$. This construction is motivated by Stein variational gradient descent, studied later in Stein Variational Gradient Descent: the kernel makes the velocity field computable from particles, while its regularity is inherited from that of the kernel. The price is a much more restrictive transport geometry.

This is the vector-valued analogue of the scalar RKHS norm used for MMDs in Dual RKHS Norms and Maximum Mean Discrepancies. The action itself is independent of $\alpha$; the measure only enters through the continuity equation $\partial_t\alpha_t+\diverg(\alpha_t v_t)=0$, which says how the common velocity field moves all particles. This type of Stein geometry was introduced in the analysis of SVGD by Liu and Wang Liu & Wang, 2016Liu, 2017 and later developed geometrically in Duncan et al., 2023Nüsken & Renger, 2023. The important caveat is that the admissible tangent space is the restricted RKHS class, not the whole Wasserstein tangent space; its regularity depends on $k$.