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. 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.

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 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 (15). The theorem of Benamou and Brenier states that this geodesic energy is exactly the squared Wasserstein distance Benamou & Brenier, 2000.

Convex Moment-Based Reformulation¶

Although (20) 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 (25). 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 (28), 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.

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 distances on spaces of probability measures by modifying the action minimized in the Benamou--Brenier formula. The objects introduced here are metric: they specify admissible curves, tangent variables and path energies. 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. When this action is normalized as a squared infinitesimal speed, it generates the length-space value

Equivalently, one may quotient by velocity fields that induce the same first-order variation of the measure. The formula above should be read as a dynamic definition of the distance, not as a property automatically satisfied by an arbitrary discrepancy. 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 least-action distance generated by this tensor is

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. Let $I\subset[0,+\infty)$ be a convex interval and let $\theta:I\to[0,+\infty)$ be concave. Define

The corresponding velocity action is

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

and it is set to $+\infty$ when $\alpha\not\ll\lambda$. 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 (40). 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 associated value is therefore written

where 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 (38) 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.

Let $\gamma$ be a monotone spectral gauge on $\mathbb S_+^d$. For a probability measure $\alpha$ and a velocity field $v\in L^2(\alpha;\RR^d)$, define the spectral tangent action

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. With the length-distance notation introduced in (38), the associated dynamic action distance is

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$. It becomes local only for linear spectral gauges. For instance, if $\gamma(M)=\operatorname{tr}(GM)$ with $G\succeq0$, then the velocity and momentum densities are

and the trace gauge, $G=\Id$ in $\gamma(M)=\operatorname{tr}(GM)$, recovers the usual Benamou--Brenier action.

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 smooth and computable from particles, at the price of defining a much more restrictive transport geometry.

Let $k$ be a positive definite kernel on $\RR^d$ with scalar RKHS $\RKHS_k$. The vector-valued RKHS is

This is the vector-valued analogue of the scalar RKHS norm used for MMDs in Dual RKHS Norms and Maximum Mean Discrepancies. The specific kernelized tangent action is

where the general distance formula (38) is understood with the restricted admissible tangent class $v_t\in\RKHS_k^d$. 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 smooth 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 smooth RKHS class, not the whole Wasserstein tangent space.

One should read $\mathcal W_k$ as an extended distance on finite-action components, not as a replacement for $\Wass_2$ on all of $\Pp_2(\RR^d)$. A useful sufficient condition for finiteness is that the endpoints lie on the same RKHS-flow orbit: if there exists $v\in L^2([0,1];\RKHS_k^d)$ whose flow map $\Phi_t$ solves $\dot\Phi_t=v_t\circ\Phi_t$ and satisfies $\alpha_1=(\Phi_1)_\sharp\alpha_0$, then $\mathcal W_k^2(\alpha_0,\alpha_1)\leq\int_0^1\norm{v_t}_{\RKHS_k^d}^2\d t$. In particular, for strictly positive definite kernels, two discrete measures with the same weights and distinct moving support points are at finite distance whenever their atoms can be connected by noncolliding smooth paths, because RKHS interpolation constructs vector fields realizing the prescribed atom velocities along the paths.

The same condition also explains the limitation. If $k$ is smooth enough that $\RKHS_k^d$ embeds into Lipschitz vector fields, finite-action curves are induced by regular flows. Atomic measures remain atomic with the same number of atoms, so a Dirac mass cannot be transported at finite kernelized action to a measure with a density. This lack of splitting is precisely what makes the geometry useful for deterministic particle methods, and also what makes it a nontrivial extended object rather than a full probability metric.

Nonlocal Wasserstein Distances¶

The local dynamic distances above transport mass through a vector field on the base space and a classical continuity equation. Nonlocal geometries use a different tangent model: the elementary motion is an exchange across an edge or a jump from $x$ to $y$. The common data are a reversible kernel $K$, a symmetric edge or jump measure $\mathsf J$, a pairwise increment $\bar\nabla$, and a pair-space action $\mathbb A_K$. The tangent variable is therefore attached to pairs of points, not to a single point $x$, so these constructions are not simply obtained by choosing another pointwise local action-density $A(\rho(x),m(x))$.

There are two complementary versions. On a finite state space, the goal is to put a Wasserstein-like geometry on the probability simplex so that the entropy gradient flow is exactly a prescribed reversible Markov chain Maas, 2011Mielke, 2013Chow et al., 2012. On a continuum space, the same edge calculus becomes a jump calculus over $\mathcal X\times\mathcal X$, which models nonlocal motion, heavy-tailed jumps, and fractional-type diffusion; this is the construction of Erbar Erbar, 2014, building on nonlinear mobilities Dolbeault et al., 2009, with subsequent metric and asymptotic refinements Slepčev & Warren, 2022.