Beyond Comparing Measures

This chapter leaves the setting of scalar measures on a common ambient space. Vector- and matrix-valued OT transports mass with internal degrees of freedom, Gromov--Wasserstein compares metric-measure spaces without a prescribed correspondence, and quantum OT replaces scalar couplings by positive operators. In each case, the transport plan must also encode structure carried by the support, the fibers, or the non-commutative state space.

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

Vector and Matrix-Valued Measures¶

Scalar OT transports a nonnegative density. In imaging, color processing, spectral analysis, diffusion tensor imaging and quantum-inspired models, the object attached to a point can instead have several nonnegative components or a positive semidefinite matrix. The first step beyond scalar OT is the positive vector-valued case: the fiber remains linear and commutative, but the transport cost may couple its channels.

Positive Vector-Valued Measures¶

The simplest way to keep internal structure in transport is to attach several nonnegative masses to each spatial point and to decide whether these channels move independently or interact through the cost.

This models multi-channel densities such as spectral bins or several species transported on the same domain. In a conservative model the mass of each channel is preserved, so one assumes $\al_0^k(\X)=\al_1^k(\X)$ for every $k$. The natural vector-valued extension therefore starts from the positive cone $\RR_+^m$.

For the dynamic formulas below, assume that $\X$ is either $\RR^d$, the flat torus, or a bounded convex domain in $\RR^d$ equipped with no-flux boundary conditions. First suppose that the endpoints and the curve have densities. The direct analogue of Benamou--Brenier fixes a vector density $u_t(x)\in\RR_+^m$ and a spatial flux $V_t(x)=(V_{t,1},\ldots,V_{t,d})\in(\RR^m)^d$, where $V_{t,\ell}^k$ is the velocity of channel $1\leq k\leq m$ in spatial direction $1\leq\ell\leq d$. The conservative vector transport cost associated with an action density $\Phi$ is

subject to the endpoint constraints $u_0\d x=\al_0$, $u_1\d x=\al_1$ and the componentwise continuity equation

Thus each component satisfies its own continuity equation, but the cost may still couple the components. Singular curves are handled as in scalar dynamic OT by replacing densities and fluxes by measures and using the lower semicontinuous perspective recession convention.

The notation $\mathcal W_\Phi$ does not by itself imply a distance for an arbitrary action $\Phi$. Sufficient structural conditions are reversibility $\Phi(u,-V)=\Phi(u,V)$ and quadratic time scaling $\Phi(u,sV)=s^2\Phi(u,V)$ for $s\geq0$, together with lower semicontinuity and coercive nondegeneracy. The latter means that curves with vanishing action cannot join two distinct endpoints. Time reversal then gives symmetry, while concatenation with an optimal allocation of time gives the triangle inequality for the square root of (2). Thus $\mathcal W_\Phi$ is an extended distance on each finite-action component. Equality of all channel masses is necessary for finiteness; without these assumptions, the formula is only a dynamic transport cost and may be asymmetric or degenerate.

A simple quadratic family is obtained from a mobility matrix $\mathsf M(u)\in\mathbb S_+^m$:

with the usual convention that the value is finite only when each $V_\ell$ belongs to the range of $\mathsf M(u)$. If $u\mapsto\mathsf M(u)$ is linear and takes values in $\mathbb S_+^m$, this action is jointly convex and positively one-homogeneous in $(u,V)$. Indeed, the matrix fractional map $(M,z)\mapsto z^\top M^\dagger z$ is jointly convex on its effective domain, and $\mathsf M(su)=s\mathsf M(u)$ for $s>0$. For $m=1$ and $\mathsf M(u)=u$, one recovers exactly the scalar Benamou--Brenier action. For

the channels move independently. Non-diagonal mobilities are the simplest way to couple the coordinates while keeping the same componentwise conservation law. For instance, with $q=m^{-1/2}(1,\ldots,1)$ and $\kappa\geq0$,

increases the mobility in the common channel direction $q$ while leaving transverse directions controlled by the diagonal part. The local cost of moving one component can therefore depend on the densities and velocities of the other components, even though each component mass remains conserved.

For the linear mobility actions above, these requirements hold: the action is even and quadratic in $V$, while the range convention and finite-dimensional coercivity make it separating on its effective domain. Their lower-semicontinuous measure extensions therefore define extended distances on the finite-action components with prescribed channel masses. The diagonal choice, writing $m_k=\al_0^k(\X)=\al_1^k(\X)$, gives the weighted product metric whose squared value is $\sum_{k:m_k>0}m_k\Wass_2^2(\al_0^k/m_k,\al_1^k/m_k)$.

The conservative positive-cone model above is the basic extension of Benamou--Brenier. Adding a source term $\partial_t u+\nabla\cdot V=S$ and a convex perspective penalty in $S$ gives unbalanced or reaction--transport variants Maas et al., 2015Maas et al., 2016Dolbeault et al., 2009Mielke, 2013. These generalized transport models include dissipation and density modulation. The figure below contrasts the exact diagonal case $\kappa=0$, where each positive channel is transported by its quantile map, with a large-$\kappa$ illustrative common-mode interpolation in which the channels move more coherently. The endpoints are two-mode mixtures: at each spatial mode the two channels have Gaussian profiles with the same center but different amplitudes.

Figure Div contrasts the exact diagonal case $\kappa=0$, where each positive channel is transported by its quantile map, with a large-$\kappa$ illustrative common-mode interpolation in which the channels move more coherently.

show_book_figure("vector-valued-measure-geodesics")

One-dimensional positive $\RR_+^2$-valued transport displayed by arrow glyphs at eight time levels. Each endpoint is a mixture of two localized Gaussian modes, and, inside each mode, both channel profiles have the same center. Each arrow is proportional to the local fiber value $(u_t^1(x),u_t^2(x))$, and time runs vertically from the red source to the blue target. Left: for $\kappa=0$, the diagonal mobility gives two independent scalar quantile geodesics. Right: a large-$\kappa$ common-mode interpolation bends the display toward $q=2^{-1/2}(1,1)$, illustrating the effect of a mobility that favors coherent channel motion while keeping the same componentwise continuity equation.

The interactive demo keeps the same glyph idea and lets the coupling strength bend the fibers toward a common channel direction.

Positive Matrix-Valued Measures¶

The next simplest fiber is the positive matrix cone. This is the simplest tensor-valued model beyond vectors: the diagonal entries behave like positive channels, while the eigenvectors encode local orientations.

Equivalently, $\mathcal A$ is a symmetric matrix of finite signed measures such that $\mathcal A(E)\succeq0$ for every Borel set $E$. If $\mathcal A$ has density $A(x)\succeq0$, then $\operatorname{tr}A(x)$ is the scalar amount of mass at $x$, while, wherever $\operatorname{tr}A(x)>0$, the normalized matrix $A(x)/\operatorname{tr}A(x)$ records an internal covariance or orientation. This is the matrix analogue of the positive vector case: diagonal matrices encode nonnegative vector components, and non-diagonal matrices add a local eigenbasis.

The conservative Benamou--Brenier model fixes a matrix density $A_t(x)\in\mathbb S_+^m$ and symmetric matrix fluxes $P_t(x)=(P_{t,1},\ldots,P_{t,d})\in(\mathbb S^m)^d$. With no flux through the boundary of $\X$, the full matrix mass $\int_\X A_t(x)\d x$ is conserved, so the endpoints must have the same total matrix. The model minimizes the matrix-perspective action

subject to $A_0\d x=\mathcal A_0$, $A_1\d x=\mathcal A_1$ and to the matrix-valued continuity equation

Here $A^\dagger$ denotes the Moore--Penrose inverse, with the usual lower-semicontinuous perspective convention: the action is finite only when the columns of each $P_{t,\ell}$ belong to the range of $A_t$. The matrix fractional map $(A,P)\mapsto\operatorname{tr}(P^\top A^\dagger P)$ is jointly convex on this effective domain and positively one-homogeneous in $(A,P)$. This gives the simplest non-trivial matrix-valued transport model: spatial motion is conservative, but the fiber carries orientation through the eigenvectors of $A_t(x)$.

Here $\mathcal W_{\mathrm{mat}}$ is an extended distance. The action is invariant under $P\mapsto-P$, quadratic under time rescaling, and zero only for $P=0$ on its effective domain. These properties give symmetry, the triangle inequality, and definiteness. The value is finite only between endpoints joined by a finite-action curve; equality of their total matrices is necessary but need not by itself guarantee finiteness. Consequently, $\mathcal W_{\mathrm{mat}}$ is a genuine distance on each finite-action component and is $+\infty$ between distinct components.

The restriction to a fixed diagonal basis gives eigenvalue transport; it should be read as a commuting submodel, not as a claim that non-diagonal excursions can never change the unrestricted value. The genuinely matrix-valued case starts when the eigenspaces vary with $x$ or along the interpolation, so that the transported object carries both mass and orientation. Static matrix-valued Monge--Kantorovich problems and dual test-function metrics were developed in Ning & Georgiou, 2014Jiang et al., 2012Ning et al., 2015; dynamic versions and related non-commutative geometries appear in Chen et al., 2016Chen et al., 2020Carlen & Maas, 2014Peyré et al., 2019. The figure below shows the analogous independent/coupled contrast for positive $2\times2$ matrix fibers, using two localized matrix modes whose eigenvalue profiles share a common center at each mode.

Figure Div shows the analogous independent/coupled contrast for positive $2\times2$ matrix fibers, using two localized matrix modes whose eigenvalue profiles share a common center at each mode.

show_book_figure("matrix-valued-measure-geodesic")

Positive $2\times2$ matrix-valued transport on a one-dimensional base. Each endpoint is a mixture of two localized matrix modes; within one mode, both eigenvalue profiles are Gaussian bumps with the same center. Each ellipse is the glyph of a positive semidefinite matrix $A_t(x)$, with axes given by eigenvectors and eigenvalues. Left: the matrices are diagonal in a fixed basis, giving the commuting tensor analogue of independent vector channels. Right: a coupled illustrative interpolation bends packet motion toward the trace-density transport and uses non-commuting eigendirections; the superposition remains positive semidefinite and produces spatially varying orientations.

Wasserstein Over Wasserstein¶

The Wasserstein construction can itself be iterated: once $(\X,d)$ is a metric space, $(\Pp_p(\X),\Wass_p)$ is again a metric space and can serve as the ground space for another Wasserstein distance. This is useful when the objects to compare are random probability measures, or mixtures whose components are meaningful objects rather than only a collapsed density. We first record the structure of this new ground space, then distinguish a law of measures from its collapsed mixture, and finally transport such laws.

Wasserstein Spaces as Ground Spaces¶

The standard setting is that of Polish spaces, introduced in Definition: Polish Metric Space. These assumptions provide separability, completeness, tightness criteria and regular conditional probabilities. The next proposition shows that Wasserstein spaces preserve this well-behaved structure.

Random Measures and Collapsed Mixtures¶

Fix $1\leq p<\infty$. Elements of $\Pp_p(\Pp_p(\X))$ are probability laws over probability measures, or random probability measures. A measurable parametric family gives

This law belongs to $\Pp_p(\Pp_p(\X))$ precisely when, for one and hence every $x_0\in\X$,

The term “barycentric mixture” should not be confused with the Wasserstein barycenter of Definition Definition: Optimal-Transport Barycenter in Section OT Barycenters. The collapse (19) is linear in $\mathfrak A$, whereas a Wasserstein barycenter gives a nonlinear way to reduce a law over measures to one measure by minimizing its average transport cost.

If $\mathfrak A\in\Pp_p(\Pp_p(\X))$, then $\bar\alpha_{\mathfrak A}\in\Pp_p(\X)$ because

Transport Between Random Measures¶

Having equipped $\Pp_p(\X)$ with the ground distance $\Wass_p$, one applies the ordinary Kantorovich construction to probability laws on this space:

For $p=2$, Gaussian mixtures provide an explicit example with two geometries. A mixture can be viewed as a collapsed density on $\X$, or as a component law over Gaussian atoms in the Bures--Wasserstein space. For two component laws

the component-level problem uses the cost

If $\P^\star$ is an optimal coupling between the weights $a$ and $b$, and if $A_{ij}$ is the Brenier linear part from $\Sigma_i$ to $\Lambda_j$, each active pair follows the Gaussian geodesic

Collapsing these component geodesics gives

This component-level interpolation generally differs from the true $\Wass_2$ interpolation between the collapsed mixture densities.

Figure Div contrasts this component-level geodesic with the true one-dimensional Wasserstein interpolation of the collapsed mixtures, making the internal mass splitting absent from the former visible.

show_book_figure("kantorovich-wow-mixtures")

Two interpolations between the same asymmetric three-component one-dimensional Gaussian mixtures. The red endpoint has a broad central component carrying most of the mass, while the blue endpoint has two dominant sharp side modes. Left: Gaussian components are transported as atoms using their Bures--Wasserstein distance. Right: the collapsed densities are interpolated by the true one-dimensional quantile formula for $\Wass_2$. The central mass is split and recombined in the collapsed geometry, making the two paths visibly distinct.

The interactive comparison keeps both geometries side by side: component-level transport moves Gaussian atoms, while collapsed transport rearranges the full density.

This viewpoint also clarifies lower bounds for Gromov--Wasserstein distances: a metric-measure space can be mapped to a law of local distance profiles, and these laws can be compared by Wasserstein-over-Wasserstein.

Gromov--Wasserstein¶

Gromov--Wasserstein compares spaces through their internal distance structures rather than through a fixed ambient ground cost. This is the right extension for graphs, shapes and point clouds whose points are not pre-aligned.

Discrete Formulation¶

Optimal transport needs a ground cost $\C$ to compare histograms $(a,b)$, and thus cannot be used directly if the histograms are not defined on the same underlying space, or if one cannot pre-register these spaces to define a ground cost. Instead, assume that two matrices $D\in\RR^{n\times n}$ and $D'\in\RR^{m\times m}$ represent relationships between points. A typical scenario is when these matrices are powers of distance matrices. Define the quadratic distortion and its minimum by

where $p\geq1$ and $\Delta$ is usually $\Delta(u,v)=|u-v|$. This is a non-convex quadratic problem over the transport polytope. In the uniform case with $m=n$ and $\P$ constrained to be a permutation matrix, it becomes a Quadratic Assignment Problem, already NP-hard in full generality Loiola et al., 2007. The relaxed coupling formulation can therefore be read as a soft graph-matching model Lyzinski et al., 2016.

Figure Div shows this intrinsic matching principle under progressively stronger deformations: correspondences are selected from within-space distance patterns rather than an ambient cross-space cost.

show_book_figure("gromov-isometry-matching")

Gromov--Wasserstein correspondences under increasing deformation. The red and blue point clouds are not compared through an ambient Euclidean cross-cost; instead, the GW coupling compares their internal pairwise distances. A perfectly isometric copy admits a clean structural match, while mild and deliberately stronger deformations progressively bend the correspondence.

The interactive demo uses a fixed structural correspondence and lets the deformation change the pairwise-distance residual. This isolates the quantity minimized by the GW objective.

When $D,D'$ are genuine distance matrices, the construction below defines a distance between metric spaces equipped with a probability distribution, up to measure-preserving isometries Mémoli, 2011Sturm, 2012Schmitzer & Schnörr, 2013. The same construction also explains why GW satisfies the triangle inequality after quotienting by isometries, and its relation to Hausdorff and Gromov--Hausdorff distances is discussed at the end of the section.

General Setting¶

The continuous formulation abstracts the discrete distance matrices into metric-measure spaces, so that GW compares intrinsic geometries independently of labels, parametrizations or ambient coordinates.

Compactness is often assumed below to avoid additional tightness and integrability arguments.

For metric-measure spaces $\mathbb X=(\X,d_\X,\alpha)$ and $\mathbb Y=(\Y,d_\Y,\beta)$, define

The metric structure also gives geodesics. Sturm’s construction allows one to speak about interpolation, barycenters and gradient flows directly on the space of metric-measure spaces, even though the intermediate space lives on a product support and is therefore expensive numerically Sturm, 2012.

Figure Div complements the global GW objective with local diagnostics, displaying where a mildly non-isometric correspondence creates the largest pairwise-distance residuals.

show_book_figure("gromov-nonisometric-distortion")

Local distortion in a mildly non-isometric GW match. The left panel colors transport segments by the average residual induced by the displayed hard correspondence. The right panel shows the pairwise-distance residual matrix $|d_\X(x_i,x_{i'})-d_\Y(y_{\sigma(i)},y_{\sigma(i')})|$, with darker entries marking larger local distortion. This matrix is the local contribution minimized by the discrete GW objective for the displayed correspondence.

Mapping each point to its distance-profile distribution and pushing forward the ambient measure produces a distribution over distributions. This is a direct instance of the Wasserstein-over-Wasserstein geometry defined in (21) and developed in Section Wasserstein Over Wasserstein.

The next figure exposes the two nested transport problems for the planar shapes $\mathbb X$ and $\mathbb Y$: sorting each distance profile computes the one-dimensional costs, then an outer assignment couples the resulting profile laws.

Figure Div makes the two transport levels explicit for the planar shapes $\XX$ and $\YY$: sorting each distance profile computes the one-dimensional costs, then an outer assignment couples the resulting profile laws.

# Exact sorted profiles; the histograms are display summaries only.
show_book_figure("gromov-memoli-distance-profiles")

Mémoli distance profiles expose a computable lower bound for intrinsic GW comparison. The planar shapes $\mathbb X$ and $\mathbb Y$ are represented by cat and bunny silhouettes, centered and normalized to unit diameter. Matching colors identify representative anchor pairs selected by the optimal profile assignment with $C_{ij}=\Wass_2^2(\alpha_{x_i},\beta_{y_j})$, their connecting segments and their two side histograms. The histograms are display summaries only: every profile cost is computed by sorting the complete distance profiles. The resulting outer assignment realizes the Mémoli profile lower bound.

This lower bound is useful computationally because the profile cost matrix $\C_{ij}=\Wass_p(\alpha_{x_i},\beta_{y_j})^p$ is an ordinary OT cost between points. Solving this easier OT problem gives a geometry-aware initialization for the non-convex GW iterations.

Relation With Wasserstein-Procrustes¶

The profile lower bound is intrinsic. In Euclidean applications, it is naturally paired with an extrinsic upper certificate obtained by registering the two measures before applying the ordinary Wasserstein distance. Proposition Proposition: Wasserstein-Procrustes Upper Certificate, proved in Quotient Wasserstein and Wasserstein-Procrustes, supplies exactly this certificate. The converse need not hold, because a small GW value may be achieved by an intrinsic correspondence that is not induced by any ambient rigid motion. Combining the profile lower bound with the Procrustes upper certificate gives the sandwich

where $\mathbb X=(\RR^d,\norm{\cdot},\alpha)$ and $\mathbb Y=(\RR^d,\norm{\cdot},\beta)$. The left term is intrinsic and inexpensive; the right term is an ambient rigid-registration certificate.

Entropic Regularization and Fused GW¶

For the common squared distortion $\Delta(u,v)^2=(u-v)^2$, one often seeks a stationary point of the entropic relaxation

For symmetric distance matrices, define the half-gradient

A standard fixed-point linearization Peyré et al., 2016 computes

The factor $\epsilon/2$ is essential because $\C(P)$ is one half of the quadratic gradient. Each update is an ordinary entropic OT problem and can therefore be solved with Sinkhorn iterations. If the iterates converge to a positive fixed point, it satisfies the stationarity conditions of the regularized GW objective. The basic fixed-point iteration is not a descent method in general and has no global guarantee for this non-convex problem; line searches or proximal variants are needed when monotone decrease is required.

Fused Gromov--Wasserstein augments the structural term with a feature transport cost Vayer et al., 2019. In the discrete case, given a cross-feature cost $M\in\RR^{n\times m}$ and a parameter $\lambda\in[0,1]$, one minimizes

The endpoints $\lambda=0$ and $\lambda=1$ recover feature-only OT and pure GW respectively; intermediate values trade attribute matching against structural matching. The first term compares node attributes in the usual OT sense, and the second compares intrinsic geometry; this is useful when two spaces have both distances and features, and the two sources of information may disagree.

Figure Div isolates this tradeoff on a small graph pair by comparing feature-only, structure-only and fused correspondences.

show_book_figure("fused-gromov-feature-geometry")

Feature information and intrinsic geometry in fused Gromov--Wasserstein. Small inner disks encode binary node features. Feature-only OT follows the attributes even when this crosses the shape structure, pure GW follows the intrinsic ordering, and fused GW balances the feature term with the pairwise-distance distortion.

Hausdorff and Gromov--Hausdorff Viewpoints¶

If $A,B$ are compact subsets of a common metric space $(\mathcal Z,d_\mathcal Z)$, their Hausdorff distance is

The Gromov--Hausdorff distance removes the common ambient space by minimizing this quantity over all isometric embeddings into a third space:

Equivalently, it is half the minimal distortion of a correspondence between $\X$ and $\Y$ Gromov, 2001Mémoli, 2007. This is a worst-case set distance: every point must be matched with small distortion. Gromov--Wasserstein replaces correspondences by probability couplings and worst-case distortion by averaged distortion. It is therefore better adapted to noisy sampled shapes and weighted graphs, but it can ignore small sets of mass that would dominate the Hausdorff distance.

Quantum Optimal Transport¶

Quantum optimal transport replaces probability vectors by density matrices and scalar couplings by positive operators on a tensor product space. This is the right language when the transported objects are matrix-valued signals, covariance-like descriptors or quantum states, and it exposes a precise bridge between OT, non-commutative entropy and operator scaling Ning & Georgiou, 2014Chen et al., 2016Chen et al., 2020Peyré et al., 2019Caglioti et al., 2020Chakrabarti et al., 2019.

Finite-Dimensional States and Couplings¶

For $A\in\mathbb H_n$ and $B\in\mathbb H_m$, the tensor product $A\otimes B$ is the linear operator determined by

In product bases $(e_i)_{i=1}^n$ and $(f_j)_{j=1}^m$, the operator $T$ can equivalently be viewed as the four-index tensor

The full and partial traces then read

In particular, $\operatorname{tr}(T)=\operatorname{tr}(\operatorname{Tr}_B T) =\operatorname{tr}(\operatorname{Tr}_A T)$. The feasible set below is never empty, since $A\otimes B$ has partial traces $A$ and $B$ whenever $A$ and $B$ are density matrices.

The feasible set is never empty, since $A\otimes B$ has marginals $A$ and $B$.

The dual potentials have the usual scalar gauge freedom: replacing $(F,G)$ by $(F+tI_n,G-tI_m)$ leaves both the constraint and the value unchanged because $\operatorname{tr}(A)=\operatorname{tr}(B)=1$.

Entropic Regularization and Bregman Iterations¶

For $\epsilon>0$ define

This is the non-commutative analogue of entropic OT: the Shannon entropy of a coupling is replaced by the trace entropy of a density matrix Peyré et al., 2019Chakrabarti et al., 2019.

Writing $K=\exp(-C/\epsilon)$, the objective differs by a constant from $\epsilon$ times the quantum KL divergence

The exact quantum analogue of Sinkhorn is an implicit alternating Bregman projection scheme onto the affine marginal sets

In the diagonal case this proposition gives the usual multiplicative Sinkhorn updates. In the non-commutative case, however, the exact block equations

do not admit scalar division formulas, because the exponential of $F\otimes I_m+I_n\otimes G-C$ cannot be separated unless the local potential commutes with the cost.

Gurvits Scaling and Quantum Sinkhorn¶

The exact Bregman scheme above has a major computational obstruction: unlike classical marginal rescaling, projection onto either non-commutative marginal constraint requires solving a nonlinear matrix equation and has no closed form in general. Gurvits scaling instead provides an explicit fixed-point construction of a feasible quantum coupling, used as a surrogate for entropic QOT. Its iterates do not minimize the entropic QOT objective and, outside special commuting settings, do not generally arise as block minimizations of a global energy. In the diagonal case, however, the construction coincides with classical Sinkhorn and hence recovers its variational interpretation.

The algorithm often called quantum Sinkhorn comes from the operator-scaling literature of Gurvits and subsequent developments Gurvits, 2003Gurvits, 2004Georgiou & Pavon, 2015Garg & Oliveira, 2018. It replaces the true Gibbs coupling (64) by the symmetric factorization

where $U=\exp(F/(2\epsilon))$, $V=\exp(G/(2\epsilon))$ and $K=\exp(-C/\epsilon)$. For two matrices $M,N$, their commutator is $[M,N]\eqdef MN-NM$. Thus $[Z,C]=0$ means that $Z$ and $C$ commute; in this case $T_s(F,G)=T_e(F,G)$, whereas otherwise $T_s$ is a Strang-type symmetric surrogate.

To make the Choi identification unambiguous, use the product bases fixed above and associate with $K$ the positive linear map $\mathcal K:\mathbb H_m\to\mathbb H_n$ and its Hilbert--Schmidt adjoint $\mathcal K^\star:\mathbb H_n\to\mathbb H_m$ through

Equivalently,

The adjoint relation is $\operatorname{tr}(Y\mathcal K(X)) =\operatorname{tr}(\mathcal K^\star(Y)X)$. In the standard Choi convention, the completely positive map with Choi matrix $K$ is $X\mapsto\mathcal K(X^\top)$; convention (72) merely places this basis-dependent transpose in the Choi identification rather than in the marginal equations. With this convention, the marginal equations for the symmetric coupling are exactly

and can be enforced by the congruence normalizations

These inverse square roots are well-defined when $K\succ0$ and $U,V,A,B\succ0$. Under the standard strict-positivity and scalability hypotheses, the alternating normalizations converge to the prescribed marginals Georgiou & Pavon, 2015Garg & Oliveira, 2018. At finite tolerance they return an approximate coupling. When all matrices are diagonal, the updates reduce to classical Sinkhorn scaling; when the targets are proportional to identities, they match the usual bistochastic operator-scaling normalization up to trace convention.

Dynamic Time Warping¶

Dynamic time warping (DTW) compares ordered feature sequences when the same phenomenon may be observed under different clocks. It is historically rooted in speech recognition Vintsyuk, 1968Sakoe & Chiba, 1978, and is now a standard tool for time-series alignment, retrieval and classification Berndt & Clifford, 1994Müller, 2007. Like OT, it minimizes an aggregate feature mismatch over correspondences; unlike OT, those correspondences must respect chronology.

Ordered Alignments Versus Transport Couplings¶

For two empirical measures, Kantorovich OT minimizes a linear cost over the convex polytope of nonnegative matrices with prescribed row and column sums. DTW instead minimizes over the finite, non-convex set of connected monotone paths through the pairwise cost matrix. Every time index must be visited, but it may be visited repeatedly; the row and column sums therefore record endogenous visit counts rather than prescribed masses. Normalizing a path matrix produces a coupling only for these path-dependent marginals, not for fixed input histograms. Conversely, ordinary OT between the unordered empirical feature measures forgets chronology and may match indices in a crossing order. Temporal penalties, causal constraints, and joint OT--DTW models interpolate between the two viewpoints; spatio-temporal alignment, for example, combines regularized OT for spatial comparison with soft-DTW for chronological alignment Janati et al., 2020. Here “dynamic” refers to Bellman’s dynamic programming on the index grid, not to the transport PDEs of Paragraph.

Discrete Variational Problem¶

Let $x=(x_i)_{i=1}^n$ and $y=(y_j)_{j=1}^m$ be two sequences in a feature space $\mathcal Z$, and set $\C_{ij}=c(x_i,y_j)$ for a nonnegative cost $c:\mathcal Z\times\mathcal Z\to\RR_+$. A warping path is a sequence

that starts at $(1,1)$, ends at $(n,m)$, and has increments

Denote the set of such paths by $\Omega_{n,m}$ and the corresponding incidence matrix by $(A_\omega)_{ij}=\mathbf1_{\{(i,j)\in\omega\}}$. Its length satisfies $\max\{n,m\}\leq L\leq n+m-1$, its total mass is $\sum_{ij}(A_\omega)_{ij}=L$, and its two marginals are precisely the row and column visit counts.

The definition is symmetric when $c$ is symmetric, but it is generally not a metric: repetitions can give zero cost to distinct sequences, and the triangle inequality can fail. Step weights, slope constraints and a Sakoe-Chiba band are common variants that penalize excessive repetition or restrict the admissible temporal distortion Sakoe & Chiba, 1978Müller, 2007.

Dynamic Programming¶

The monotone path structure converts the exponentially large variational problem (79) into a shortest-path computation on an acyclic grid.

Continuous Time Warping¶

Discrete DTW depends on the sampling density because every visited cell contributes once. A direct continuous registration minimizes $\int_0^1 c(x(t),y(\gamma(t)))\d t$ over nondecreasing endpoint-fixing maps $\gamma$, but this one-clock formulation is asymmetric and privileges the parameterization of $x$. Continuous DTW instead traverses both clocks and measures mismatch per unit length in the parameter square Buchin et al., 2022.

For simplicity, let $x,y:[0,1]\to\mathcal Z$ use normalized clocks, and let $\Gamma_\uparrow$ contain pairs $(\phi,\psi)$ of absolutely continuous, nondecreasing surjections of $[0,1]$ onto itself. Equivalently, the endpoint conditions are $\phi(0)=\psi(0)=0$ and $\phi(1)=\psi(1)=1$. The continuous DTW functional is

Because $\dot\phi,\dot\psi\geq0$ almost everywhere, the last factor is the $\ell^1$ line element of the monotone path $s\mapsto(\phi(s),\psi(s))$. Formula (81) is invariant under increasing reparameterizations of the auxiliary variable $s$ and does not privilege either clock; when $c$ is symmetric, the resulting functional is also symmetric in $x$ and $y$. After parameterization by $\ell^1$ arc length, it is simply the line integral of the feature mismatch along a monotone path. For physical clock intervals $[0,p]$ and $[0,q]$, the same formula uses $\phi:[0,1]\to[0,p]$ and $\psi:[0,1]\to[0,q]$. Exact computation is substantially harder than the discrete recurrence: for arc-length-parametrized one-dimensional polygonal curves and the standard cost $c(u,v)=|u-v|$, Buchin, Nusser and Wong propagate piecewise-quadratic boundary costs in $O((n+m)^5)$ time Buchin et al., 2022. This complexity statement does not apply to an arbitrary feature cost $c$.

Soft-DTW and the Sinkhorn Analogy¶