Computational Optimal Transport

• Optimal transport generalizes least-effort paths for single items to moving multiple items or distributions between configurations.For an individual, the least-effort route is typically a straight line in the plane or a geodesic in other metric spaces; distributional transport requires more advanced formalism, first hinted at by Monge [1781].Source: Introduction · p. 7
• Practical concerns shaped early optimal-transport work, including Monge [1781], Tolstoı [1930], Hitchcock [1941], and Kantorovich [1942].The text connects distributional transport to everyday movement of people, commodities, and information, and to economic problems such as logistics, production planning, and network routing.Source: Introduction · p. 7
• After the 1980s, work by Brenier [1991] and others revealed connections between optimal transport, convexity, partial differential equations, and statistics.At the turn of the millennium, researchers in computer, imaging, and data sciences also recognized the theory's relevance.Source: Introduction · p. 7

4 Introduction

• The book emphasizes data-science applications and fast computational schemes, while presenting the theoretical insights that support optimal transport’s practical effectiveness.Its focus includes imaging sciences and machine learning.Source: §4 · p. 8
• Chapters 2–4 and 9–10 focus on probability-vector geometry, with gray boxes extending the exposition to discrete measures and dark-gray boxes to arbitrary measures.Chapters 5–8 address the interplay between continuous and discrete measures for a more mathematically inclined audience.Source: §4 · p. 8
• The survey leaves many active topics untreated, including distributionally robust optimization [Shafieezadeh Abadeh et al., 2015, Esfahani and Kuhn, 2018, Lee and Raginsky, 2018, GAO et al., 2018], Wasserstein Langevin convergence, and low-dimensional Monge–Ampère methods [Froese and Oberman, 2011, Benamou et al., 2014, Sulman et al., 2011].The Monge–Ampère methods are only briefly mentioned in Remark 2.25.Source: §4 · pp. 8–9

Theoretical Foundations

• The framework progresses from probability vectors (histograms) to discrete measures and then arbitrary measures; histograms suffice to implement the algorithms in Chapters 3 and 4.Arbitrary-measure formulations also apply to moving point clouds and settings with samples from continuous densities.Source: Theoretical Foundations · p. 11

2.1 Histograms and Measures

• A discrete measure assigns weights to locations, while arbitrary Radon measures can represent discrete, continuous, or mixed mass distributions.For discrete measures considered here, all weights are assumed positive to avoid accounting for locations with no mass.Source: §2.1 · p. 12

2.2 Assignment and Monge Problem

• At n = 70, the assignment search space exceeds 10^100 permutations, so efficient optimization algorithms are necessary.Assignment compares equal-size sets through a bijection; the broader Monge formulation transports each source point through a map that conserves mass.Source: §2.2 · pp. 13–14
• A Monge map may not exist for discrete measures with incompatible weights, while optimal transport push-forwards measures rather than warping image functions.Even for measures with densities, the push-forward density includes a Jacobian, so it is not simply the pull-back of the target density.Source: §2.2 · pp. 14–16, Figure 2.3 · p. 17

2.3 Kantorovich Relaxation

• Kantorovich transport allows source mass to split across targets, yielding a symmetric linear program over a convex set of couplings.For arbitrary measures, couplings are joint distributions whose marginals equal the source and target measures.Source: §2.3 · pp. 17, 18, 21
• For equal-size uniform histograms, Kantorovich relaxation is tight: an optimal coupling can be chosen as a permutation matrix for an optimal matching.This follows because the extreme points of the feasible coupling polytope are permutation matrices.Source: §2.3 · p. 19
• On compact spaces with continuous cost, the Kantorovich relaxation has a solution because its feasible couplings are weakly compact and the cost functional is continuous.Source: §2.3 · p. 22

2.4. Metric Properties of Optimal Transport

• For p≥1, Wp is a metric when the cost is the pth power of a ground metric, and Wasserstein convergence characterizes weak convergence under the stated moment conditions.For 0<p≤1, Wp^p rather than Wp defines a metric.Source: §2.4 · pp. 23, 24, 26

2.5 Dual Problem

• The Kantorovich dual maximizes marginal-potential values subject to f_i+g_j≤C_ij, equals the primal cost, and its optimality conditions locate transport support where potentials saturate cost.For arbitrary measures, the support lies where f(x)+g(y)=c(x,y).Source: §2.5 · p. 30
• For quadratic Euclidean cost in R^d, if one input measure has a density, the optimal Kantorovich plan is unique and supported on the graph of a Monge map.The map is the gradient of a convex function that pushes the input measure to the other measure.Source: §2.5 · p. 31
• When the source has a density, the Kantorovich relaxation is tight and the optimal coupling is the map π = (Id, ∇ϕ)♯α for a convex potential ϕ.This extends the notion of increasing functions to higher dimensions.Source: §2.5 · p. 32
• Brenier’s theorem makes OT a basis for quantile functions in arbitrary dimensions, supporting quantile regression [Carlier et al., 2016].Source: §2.5 · pp. 32–33
• The Monge–Ampère formulation is technically challenging: solutions may be nonsmooth, and anisotropy, discretization, and boundary conditions complicate numerical solvers.For arbitrary input measures such as Dirac masses, analysis may require Alexandrov solutions; standard finite differences can be inaccurate for anisotropic convex functions.Source: §2.5 · p. 33

2.6 Special Cases

• For a binary off-diagonal cost, the 1-Wasserstein distance equals the 1-norm of the mass difference; the analogous discrete-measure distance equals total variation.Source: §2.6 · p. 34
• In one dimension, empirical measures with ordered points admit a local ℓp-distance formula, and Wasserstein geometry is isometric to an Lp space; W1 is a norm.The ordered-vector formula is local because the order can change as point values change.Source: §2.6 · pp. 34, 36
• One-dimensional OT equalizes image histograms by reassigning sorted grayscale values, producing an image whose pixel distribution matches the target and a displacement interpolation between inputs.Source: §2.6 · p. 35
• For Gaussian measures, the optimal map is linear and the Wasserstein distance is expressed through the Bures metric, which reduces to Hellinger distance for diagonal covariances; in 1-D, W2 is Euclidean in mean–standard-deviation coordinates.Source: §2.6 · p. 38
• For elliptically contoured densities, the Gaussian linear transport map and Bures-based distance formula extend with a generator-dependent scaling.The scaling is 1 for Gaussians and 1/(d+2) for uniform distributions on ellipsoids.Source: §2.6 · pp. 38–39

Algorithmic Foundations

• OT is a linear program equivalent to a minimum cost network flow problem [Korte and Vygen, 2012, p. 213, Lem. 9.3], and the simplex method enabled rigorous solutions after earlier formulations lacked provably correct algorithms.The text traces early formulations to Tolstoı, Hitchcock, and Kantorovich, and notes the connection between Dantzig’s motivation and transportation problems [Dantzig, 1949, p. 210].Source: Algorithmic Foundations · pp. 41–42

3.1 The Kantorovich Linear Programs

• The primal OT problem becomes a standard-form linear program by vectorizing P and encoding its row- and column-sum constraints with a matrix; one constraint is redundant.The formulation retains the redundant constraint, with the text noting that degeneracy may arise in computations.Source: §3.1 · p. 42

3.2 C-Transforms

• Freezing one dual vector makes its C-transform the largest feasible choice for the other and reduces the dual to a piecewise-affine concave maximization, but alternating transforms quickly reaches a plateau.The transforms also satisfy identities such as f C C̄ C = f C.Source: §3.2 · p. 44, §3.3 · p. 45

3.3 Complementary Slackness

• Feasible primal and dual solutions are both optimal when complementary; at optimality, positive transport can occur only where Cij = fi + gj.Complementarity requires equality on every positive entry of P, while dual feasibility requires Cij ≥ fi + gj everywhere.Source: §3.3 · pp. 45–46

3.4.1 Tree Structure of the Support of All Vertices of U(a, b)

• An extremal transport plan has acyclic support on the bipartite graph and therefore at most n + m −1 nonzero entries.The graph encodes source and target nodes, with each positive entry corresponding to an edge in the plan’s support.Source: §3.4.1 · p. 46, §3.4 · p. 47

3.4.2 The North-West Corner Rule

• NW corner fills each entry with the maximum feasible mass in sequence, producing a vertex in at most n + m operations; row and column permutations yield many such plans, usually a strict subset of all vertices.The rule can initialize a primal algorithm such as the network simplex.Source: §3.4.2 · p. 48, §3.5 · p. 49

3.5 A Heuristic Description of the Network Simplex

• The network simplex associates a complementary dual pair with a feasible primal plan; if the pair is feasible the plan is optimal, otherwise the plan is modified toward dual feasibility.The complementary equalities are imposed on the positive-support edges.Source: §3.5 · p. 49

3.5.1 Obtaining a Dual Pair Complementary to P

• For a vertex, complementary dual variables satisfy one equality per support edge; this system is underdetermined, with additional free variables when the support graph has multiple components.The passage’s example has 11 dual variables and 8 support-edge equalities.Source: §3.5.1 · pp. 49–50
• For each tree in G(P), the network simplex fixes a root dual variable to 0 and traverses the tree to determine a complementary dual pair.The traversal can use breadth-first or depth-first search.Source: §3.5.1 · p. 50, Figure 3.4 · p. 51

3.5.2 Network Simplex Update

• When a violating edge is added, the network simplex either recomputes dual variables without changing P or augments flow around a cycle and removes a minimum-flow negative edge.The cycle update increases positive-edge flows and decreases the flows on negatively impacted edges until one reaches zero.Source: §3.5.2 · pp. 51–52

3.5.3 Improvement of the Primal Solution

• A cycle update can only improve the primal cost when θ > 0; when θ = 0, the graph changes but P does not.If adding the violating edge creates no cycle, the primal solution remains unchanged.Source: §3.5.3 · p. 52

3.6. Dual Ascent Methods

• Dual ascent maintains dual-feasible iterates and improves their objective using sparse updates guided by a max-flow labeling procedure.Unlike simplex-type methods, it need not keep iterates at extreme points; its constructed primal solution is nonnegative but need not satisfy marginal constraints.Source: §3.6 · pp. 53–55

Complementary slackness.

• The auction algorithm maintains ε-complementary slackness at each iteration while growing a partial assignment toward a permutation.Its updates adjust dual prices and assignment indices as the assigned set expands.Source: Complementary slackness · pp. 57–58

3.7 Auction Algorithm

• The auction algorithm takes at most N = n∥C∥∞/ε steps and returns an assignment whose cost is nε-suboptimal.A naive operation count is N^3∥C∥∞/ε; ε-scaling can reduce it to N^3 log ∥C∥∞.Source: §3.7 · p. 59

3.7. Auction Algorithm

• The auction algorithm is presented as an alternative use of C-transforms, with its inequality chain invoking ε-CS, cancellation, and permutation suboptimality.The derivation also begins from the suboptimality of g.Source: §3.7 · pp. 59–60
• On low-dimensional regular Euclidean grids, multiscale strategies can significantly speed up classical linear solvers [Schmitzer, 2016a, Oberman and Ruan, 2015].Source: §3.7 · p. 60

Entropic Regularization of Optimal Transport

• Entropic regularization enables GPU-suited alternating matrix-vector updates (or matrix-matrix updates for shared supports), can avoid storing the cost matrix, and yields a smooth, automatically differentiable approximate distance.For some applications, the matrix-vector products require only kernel evaluations rather than a stored cost matrix.Source: Entropic Regularization of Optimal… · p. 61

4.1 Entropic Regularization

• The ε-strongly convex objective has a unique optimum, while entropy moves a simplex LP solution inward toward its entropic center, motivated by diffuse traffic patterns.The cited passage traces the regularization idea to transportation models intended to predict more diffuse traffic than sparse OT solutions.Source: §4.1 · p. 62, Figure 4.1 · p. 62
• As ε→0, the optimizer converges to the maximum-entropy OT solution; as ε→∞, it converges to a⊗b, and couplings become less sparse as ε grows.The product coupling corresponds to independent random variables with the prescribed marginals.Source: §4.1 · pp. 62–63
• The regularized coupling is a projection of the Gibbs kernel and defines the static Schrödinger problem [Léonard, 2014, Rüschendorf and Thomsen, 1998], converging to maximum-entropy OT as ε→0 [Léonard, 2012, Carlier et al., 2017].The passage also notes an alternative dynamic formulation over paths connecting the measures.Source: §4.1 · pp. 63, 65
• For arbitrary measures, relative entropy against α⊗β extends regularization, and replacing the reference marginals by measures with the same null sets leaves the solution unchanged.The reference measure matters through its support, not through the particular marginals chosen.Source: §4.1 · p. 65
• The regularization penalty can be interpreted as mutual information, which is nonnegative and equals zero exactly when the coupled variables are independent.Here the mutual information is the relative entropy of the coupling with respect to the product of its marginals.Source: §4.1 · pp. 65–66
• For one-dimensional marginals, entropic regularization generates copulas interpolating between independence, ξ(s,t)=st, and full dependence, ξ(s,t)=min(s,t).The passage associates the independent case with the product coupling and the fully dependent case with a coupling supported on a Monge map.Source: §4.1 · p. 66

Sinkhorn’s algorithm

• Alternating updates of scaling vectors u and v enforce the two marginals, and Sinkhorn iterations converge to the same optimal coupling despite scaling-vector nonuniqueness.The coupling evolves from the Gibbs kernel toward the solution by progressively shifting mass away from the diagonal.Source: §4.2 · p. 68

Applications and extensions

• Sinkhorn offers an efficient, scalable approximation to OT through parallelization across problems, and entropy regularization defines a differentiable loss unlike the linear-programming formulation, as shown by Cuturi [2013].When a = b, averaged projection iterations can preserve symmetry.Source: §4.2 · p. 69

Overall complexity

• O(n^2 log(n)τ^-3) operations suffice for a τ-approximate solution of unregularized OT after a rounding step, according to Altschuler et al. [2017].The rounded coupling satisfies ⟨P̂, C⟩ ≤ L_C(a, b) + τ.Source: §4.2 · p. 69

Numerical stability

• When Gibbs-kernel entries underflow to zero, Sinkhorn updates can divide by zero; computing scaling multipliers in the log domain can partly address this failure.Convergence also deteriorates as ε → 0.Source: §4.2 · p. 70

Proximal point algorithm

• A KL proximal-point scheme can approximate unregularized OT by repeatedly applying Sinkhorn with a decaying regularization parameter ε/ℓ.The approach can be efficient in small spatial dimensions when combined with adaptive sparse-grid multigrid [Schmitzer, 2016b].Source: §4.2 · p. 71

Other regularizations

• Replacing entropy with a strictly convex penalty such as a quadratic can produce sparse couplings, but requires positivity constraints and Dykstra projections and is less parallelizable than Sinkhorn.For quadratic regularization, the required Euclidean projection with positivity constraints can be solved using simplex projection algorithms.Source: §4.2 · pp. 71–72

Barycentric projection

• Barycentric projection converts a coupling into a transport map; when the unregularized coupling is supported on a unique Monge map, positive ε gives a smooth approximation useful in imaging.The projection has also been used to approximate principal geodesics in Wasserstein space.Source: §4.2 · p. 73

Convergence

• Positive-kernel contraction in the Hilbert projective metric yields global linear convergence of Sinkhorn’s iterations, while the marginal-violation rate degrades as ε → 0.Marginal constraint violations can serve as stopping criteria; the global rate may be pessimistic when an unregularized Monge map exists, and can worsen exponentially with ε for near-random costs.Source: §4.2 · pp. 74–76

4.3. Speeding Up Sinkhorn’s Iterations

• The Sinkhorn-map Jacobian is positive and maps 1_n to 1_n, giving it a dominant eigenvector 1_m with eigenvalue 1.The optimal coupling is P = diag(e^f/ε)K diag(e^g/ε).Source: §4.3 · p. 77
• Because f is defined up to a constant, the second eigenvalue 1 − κ governs the local linear rate.Source: §4.3 · p. 77
• In simple cases, including a smooth Monge map when ε = 0, κ scales like ε; the bistochastic assignment case is discussed by.Source: §4.3 · p. 77

4.3 Speeding Up Sinkhorn’s Iterations

• The main computational bottleneck is multiplication by K and K⊤, with naive implementation costing O(nm).The section describes cases where this cost can be significantly reduced.Source: §4.3 · p. 77
• Matrix-matrix operations let Sinkhorn iterations for multiple histogram pairs run in parallel, improving on the parallelism of its intrinsically GPU-friendly matrix-vector operations.The basic iterations are already GPU-friendly; the matrix-matrix formulation offers even better parallelism for computing distances at large scale.Source: §4.3 · pp. 77–78
• When indices are d-tuples and the cost is additive across their sub-indices, the Sinkhorn kernel has a separable multiplicative structure.That structure is what enables the slice-wise kernel application.Source: §4.3 · p. 78
• For n=m, applying the separable kernel along each slice evaluates Ku exactly in O(n^(1+1/d)), instead of O(n^2).This avoids instantiating the full n × n kernel matrix, whose size can grow exponentially with dimension.Source: §4.3 · p. 78

Separable convolutions

• Reshaping planar histograms as matrices lets Sinkhorn apply separate 1-D convolutions instead of a naive matrix-vector product; related shape applications handle grids as large as 200^3 = 8 × 10^6 [Solomon et al., 2015, Bonneel et al., 2016].For higher dimensions, the same 1-D convolutions are applied to tensor slices, an operation described as especially efficient on GPUs.Source: §4.3 · p. 79

Approximated convolutions

• For translation-invariant kernels, Fourier transforms provide nearly linear-time convolution approximations, but can be unstable for small ε; autoregressive filters are another option.The passage cites Deriche filtering as an example and points to a comparison of fast filtering methods.Source: §4.3 · p. 79

Geodesic in heat approximation

• For geodesic costs on nonplanar domains, Varadhan’s formulas approximate the kernel using the Laplacian for p = 1 or the heat kernel for p = 2.Discretizing the Laplace–Beltrami operator and precomputing a reordered sparse Cholesky factorization enables linear complexity per Sinkhorn iteration on 2-D manifolds; performance degrades with manifold dimension.Source: §4.3 · pp. 79–80

Extrapolation acceleration

• Linear or nonlinear extrapolation can improve Sinkhorn’s local linear convergence rate from O((1 − κ)^ℓ) to O((1 − √κ)^ℓ) for some κ > 0.The method exploits Sinkhorn’s fixed-point form to improve conditioning near the solution.Source: §4.3 · p. 80

4.4 Stability and Log-Domain Computations

• Alternating block coordinate ascent on the entropic OT dual is mathematically equivalent to Sinkhorn iterations.The dual updates cancel the gradients in f and g in turn.Source: §4.4 · pp. 81–82
• The soft-min converges to the ordinary minimum as ε → 0, but Sinkhorn iterations do not converge at ε = 0 because alternate minimization fails for the constrained unregularized dual.The soft-min is a differentiable approximation of the minimum, and the corresponding operations are equivalent to the entropic c-transform.Source: §4.4 · pp. 82–83
• Log-domain Sinkhorn iterations remain stable for arbitrary ε > 0, avoiding the original updates’ overflow or underflow risk at small ε, but require nm exponential computations per step.The log-sum-exp trick stabilizes soft-min evaluation; subtracting previously computed scalings yields a stabilized iteration.Source: §4.4 · pp. 81, 83
• Line-by-line soft-min computations are typically slower than matrix multiplications and cannot be efficiently parallelized across several marginals; decaying-ε multiscale solvers can speed computation on low-dimensional Euclidean domains.The cited multiscale approach uses sparse grids [Schmitzer, 2016b].Source: §4.4 · p. 83
• For generic measures, entropic transport smooths the constraint R(c), recovered as ε → 0, but potential existence is harder to prove because potentials are not automatically Lipschitz.Existence can be shown using convergence of Sinkhorn iterations.Source: §4.4 · pp. 83–84
• The alternative unconstrained dual is better conditioned because its functional’s Hessian is uniformly bounded by ε, but its logarithm prevents stochastic optimization.The passage contrasts this with the other unconstrained formulation, which achieves the same optimal value.Source: §4.4 · p. 84

4.5 Regularized Approximations of the Optimal Transport Cost

• The entropic dual is a smooth, unconstrained concave maximization problem that approximates the original Kantorovich dual.Source: §4.5 · p. 84
• Every optimal entropic-dual solution pair belongs to R(C), the set of feasible Kantorovich potentials.Source: §4.5 · p. 84

4.5. Regularized Approximations of the Optimal Transport Cost

• The regularized cost C(a, b) is jointly convex in a and b for ε ≥ 0.For ε > 0, its gradient is expressed using optimal solutions of the entropic dual whose coordinates sum to 0.Source: §4.5 · p. 85
• proposed lower and upper bounds for approximating Wasserstein distance between histograms by evaluating primal and dual objectives at Sinkhorn solutions.The section defines the corresponding primal and dual Sinkhorn divergences.Source: §4.5 · p. 85
• The primal-dual Sinkhorn bounds are valid only upon algorithm convergence, but using L finite iterations instead yields a lower bound on the regularized cost.The finite-step algorithmic Sinkhorn functional lower bounds the regularized cost whenever L ≥ 1.Source: §4.5 · pp. 85–86
• Primal Sinkhorn iterates are not feasible, so their cost cannot be treated as an approximation of the unregularized OT cost; rounding can instead yield an upper bound on the regularized cost.The rounding scheme can be computed with parallel matrix operations for multiple histogram pairs.Source: §4.5 · p. 86
• Unlike the regularized expression, the finite dual Sinkhorn divergence is generally nonconvex, but it is differentiable with respect to its arguments, including C, a, and b.Automatic differentiation can be used to compute these derivatives.Source: §4.5 · p. 86

4.6 Generalized Sinkhorn

• Regularized OT is a special case of structured convex optimization; indicator functions for the marginals recover the hard constraints, and the unique solution retains the scaling form (4.12).The generalized formulation uses F = ι{a} and G = ι{b} to recover U(a, b).Source: §4.6 · pp. 86–87

4.6. Generalized Sinkhorn

• Sinkhorn iterations extend to this generalized problem using a proximal operator for the KL divergence [Peyré, 2015, Frogner et al., 2015, Chizat et al., 2018b, Karlsson and Ringh, 2016].The passages identify the proximal operator but do not provide its formula.Source: §4.6 · p. 87, G (KTu) KTu , (4.51) · p. 87

G (KTu) KTu , (4.51)

• For some F and G, generalized Sinkhorn iterations converge linearly and extend to arbitrary measures.The source points to Chizat et al. (2018b) for details.Source: G (KTu) KTu , (4.51) · p. 87
• When the relevant proximal operators are closed-form or efficient to compute, separable functions reduce the updates to scalar optimization problems, and log-domain computations can stabilize them.The passages describe separable functions as sums of component functions and note that each component problem is usually simple.Source: G (KTu) KTu , (4.51) · p. 87
• The algorithm approximates solutions to unbalanced OT problems and gradient-flow problems.The cited examples are equations (10.7) and (9.26), respectively.Source: G (KTu) KTu , (4.51) · p. 87
• Generalized Sinkhorn iterations are a special case of Dykstra’s algorithm and alternate maximization on the dual problem [Dykstra, 1983, 1985, Bauschke and Lewis, 2000, Censor and Reich, 1998].The source describes Dykstra’s algorithm as extended to Bregman divergence.Source: G (KTu) KTu , (4.51) · pp. 87–88
• The formulation extends to more than two functions and one coupling, including multimarginal Sinkhorn and regularized barycenter iterations.The passages identify algorithms (10.2) and (9.18) as particular cases.Source: G (KTu) KTu , (4.51) · p. 88

5.1 c-Transform and ¯c-Transform

• The c-transform canonically extends dual potentials to whole spaces and improves them to c-concave forms; unlike entropic Sinkhorn updates, unregularized alternate maximization immediately cycles.Entropic updates replace minima with soft-min operations.Source: §5.1 · p. 90

5.2 Semidiscrete Formulation

• The semidiscrete problem becomes finite-dimensional in dual weights, whose Laguerre cells define a piecewise-constant transport map; squared Euclidean cost yields efficiently computable power diagrams.When all weights are constant, the cells reduce to a Voronoi partition.Source: §5.2 · pp. 91–92
• Semidiscrete solvers have applications in graphics, incompressible-fluid simulation, and illumination design, with the fluid formulation enforcing incompressibility at each iteration.It does so by requiring the possibly weighted point cloud to approximate a uniform distribution in the domain.Source: §5.3 · p. 93

5.3 Entropic Semidiscrete Formulation

• Entropic smoothing turns the semidiscrete objective into a smooth finite-dimensional problem related to multiclass logistic regression, enabling Newton or quasi-Newton methods.Its Hessian is bounded by 1/ε; without regularization, second-order methods are justified when Laguerre-cell volumes stay bounded below and the continuous measure has a density.Source: §5.3 · pp. 94–95

5.4 Stochastic Optimization Methods

• Stochastic optimization uses independent samples from α without discretizing it; stochastic gradient methods are provably convergent, and averaging adapts to local curvature for faster convergence.For continuous-continuous OT, stochastic dual optimization requires suitable regularization such as ε > 0; neural-network parameterizations lack approximation guarantees.Source: §5.4 · pp. 96–98

W1 Optimal Transport

• Metric-cost transport can have nonunique optimal couplings, while its distance cost is more robust to outliers and noise than a quadratic cost and yields a norm that compares signed distributions.The lack of a uniquely defined Monge map can make this setting ill-suited to measure interpolation.Source: W1 Optimal Transport · p. 100

6.1 W1 on Metric Spaces

• On metric spaces, f = gc holds exactly when Lip(f) ≤ 1, in which case fc = −f; the resulting single-potential formulation identifies W1 as a norm also applicable to signed measures.The dual reformulation uses functions with Lipschitz constant at most 1.Source: §6.1 · pp. 101–102
• In one dimension, ordering the points reduces the W1 constraints to a linear program, and W1 also has a closed-form expression using cumulative functions.For X = R with d(x, y) = |x − y|, the points can be ordered before solving.Source: §6.1 · p. 102

6.2 W1 on Euclidean Spaces

• In Euclidean space, the global Lipschitz constraint becomes a uniform gradient bound, while the dual Beckmann formulation represents transport as a flow subject to divergence constraints.After finite-element discretization, first-order methods such as Douglas–Rachford or ADMM can solve the resulting nonsmooth convex problems; on uniform grids, FFT handles projection onto the divergence constraint.Source: §6.2 · pp. 102–103

6.3 W1 on a Graph

• On graphs, W1 can be computed with local discrete gradient and divergence operators rather than a full coupling; its Beckmann formulation is a min-cost flow linear program with dedicated simplex solvers.The graph setting uses the shortest-path metric, and edge flows satisfy the discrete divergence constraint.Source: §6.3 · pp. 103–104

7.1 Continuous Formulation

• Dynamic OT minimizes the length of a measure-valued curve subject to mass conservation; a momentum-variable reformulation makes the problem convex and requires momentum to vanish where mass is zero.When an optimal Monge map exists, the interpolant is McCann’s interpolation; more generally it can be obtained by pushing an optimal coupling along straight paths.Source: Dynamic Formulations · p. 106, §7.1 · pp. 106–108

7.2 Discretization on Uniform Staggered Grids

• Staggered uniform grids discretize time, density, and momentum for dynamic OT, while graph formulations use graph divergence and may use a logarithmic mean to preserve a heat-equation gradient-flow structure.The grid scheme follows Papadakis et al. [2014].Source: §7.2 · pp. 109–110

7.3 Proximal Solvers

• Because vanishing interpolated mass makes the discretized objective nonsmooth, interior-point solvers offer high precision but costly iterations, whereas proximal first-order methods suit large structured problems.Douglas–Rachford iterations converge for 0 < α < 2 and γ > 0, and are closely related to ADMM.Source: §7.3 · pp. 110–112

7.4 Dynamical Unbalanced OT

• Dynamical unbalanced OT adds a source term to the continuity equation and penalizes it like momentum, trading off mass transport against creation or destruction.As τ → 0 with equal endpoint masses, WFR tends to classical W; as τ → +∞, it approaches the Hellinger metric over densities.Source: §7.4 · pp. 112–114

7.5 More General Mobility Functionals

• General mobility functions θ(a, b) = a^(s−p)∥b∥^p allow p ≥ 1 and 1 ≤ s ≤ p; for 1 < s ≤ p, solutions must have Lebesgue densities.The case s = 1 gives classical Wp, while s = p corresponds to a dual Sobolev norm.Source: §7.5 · p. 114, §7.6 · p. 115

7.6 Dynamic Formulation over the Paths Space

• The path-space formulation represents transport by a distribution over curves matching endpoint measures; classical optimal paths yield displacement interpolation, while entropic paths converge to classical solutions as ε → 0.Entropic interpolation is a Gaussian blurring of traveling Diracs along geodesics.Source: §7.6 · pp. 115–116
• Figure 7.5 shows Brownian bridge paths associated with the Schrödinger entropic interpolation over path space.The samples are colored blue at t = 0 and red at t = 1.Source: Figure 7.5 · p. 116

7.6. Dynamic Formulation over the Paths Space

• The entropic interpolation can be described using a regularization of the Benamou–Brenier dynamic formulation (7.2).This formulation is associated with [Gentil et al., 2015, Chen et al., 2016a].Source: §7.6 · p. 117

Statistical Divergences

• The chapter focuses on sample-based approximation of Wasserstein distance for continuous-support distributions; finite-ground-space work includes Sommerfeld and Munk [2018] and Bigot et al. [2017a].It also shows that optimal transport cannot be efficiently approximated by a Hilbertian metric on a suitable feature representation of probability measures.Source: Statistical Divergences · p. 118

8.1 ϕ-Divergences

• ϕ-divergences compare measures pointwise without transporting mass, and are jointly 1-homogeneous, convex, and weakly* lower semicontinuous.They are built from entropy functions and can serve as loss functions for estimating probability distributions.Source: §8.1 · pp. 119–120

8.2 Integral Probability Metrics

• Integral probability metrics are dual norms defined by prescribed families of test functions, yielding weaker norms as the function family is restricted.Their role as alternatives to Wasserstein distances is to compare measures through test functions rather than pointwise mass ratios.Source: §8.2 · p. 124

8.2.1 W1 and Flat Norm

• The text refers to 1-Lipschitz functions; extending Formula (6.2) with constraint |fk| ≤ 1 computes the flat norm, which metrizes weak convergence on the whole space M(X).The flat norm is sometimes called the “Kantorovich–Rubinstein” norm and has been used as a fidelity term for inverse problems in imaging.Source: §8.2.1 · p. 125

8.2.2 Dual RKHS Norms and Maximum Mean Discrepancies

• Maximum mean discrepancy is an RKHS dual norm that metrizes weak convergence when its RKHS unit ball spans a dense subset of continuous functions.For translation-invariant kernels on Rd, this density condition is equivalent to a nonvanishing Fourier transform; Gaussian kernels require bandwidth selection and can be ill-suited to multiscale data.Source: §8.2.2 · p. 126, §8.2 · p. 127

8.2. Integral Probability Metrics

• The H−1 norm is undefined for singular measures when d > 1, and its weighted version is described as a linearization of Wasserstein W2.H−1 is the dual of H1 and is defined on zero-mean measures with densities.Source: §8.2 · p. 128
• W2(α, β) is equivalent to the weighted negative Sobolev functional ∥α − β∥H−1(α), but replacing it with the unweighted H−1(Rd) expression requires bounded densities on a shared support.Negative Sobolev spaces metrize weak convergence when r>d/2, requiring higher regularity as dimension increases.Source: §8.2 · p. 128

8.3. Wasserstein Spaces Are Not Hilbertian

• Energy distance is scale-free and avoids a bandwidth parameter, but requires 0<p<2; at p=2 it degenerates to distance between means and is no longer a norm.Its scaling under dilation therefore differs from Wasserstein's perfect linear scaling.Source: §8.3 · p. 129
• For d≥2 with squared Euclidean ground cost, p-Wasserstein distances for p=1 and p=2 are not Hilbertian, as positive centered-matrix eigenvalues demonstrate.The proposition's counterexample uses measures supported on four vectors in R2; the figure shows positive maximal eigenvalues across the displayed p values.Source: §8.3 · pp. 130–131, Figure 8.6 · p. 131

8.3.1 Embeddings and Distortion

• Quasi-isometric ℓ1 embeddings are available for p-Wasserstein spaces with 0<p≤1, while sliced OT and entropic soft-min constructions provide positive-definite-kernel variants.Embedding constants grow quickly with dimension; for p=1 the embedding applies only to discrete measures, while wavelet-based embeddings can be computed approximately in linear time on uniform grids.Source: §8.3.1 · pp. 131–132, §8.3.2 · p. 132

8.4.2 Empirical Estimators for ϕ-divergences

• Empirical-measure plug-in estimates cannot consistently recover total variation or other ϕ-divergences; density estimation instead smooths samples using a bandwidth or nearest neighbors.The nearest-neighbor approach estimates the radius enclosing k samples, adapting locally to sampling density.Source: §8.4 · p. 133, §8.4.2 · p. 134, §8.5 · p. 135

8.4. Empirical Estimators for OT, MMD and ϕ-divergences

• For bounded-support measures in Rd, OT estimation rates depend on dimension and can reflect intrinsic support dimension; concentration around the mean is also established.The cited rate is tight when one measure has a density; the supplied passage does not include the rate formula.Source: §8.4 · p. 133

8.4.1 Empirical Estimators for OT and MMD

• For dual RKHS norms (MMD), sample complexity does not depend on ambient dimension, unlike Wasserstein distances.A direct empirical-measure estimator is considered for both Wasserstein distances and MMD norms.Source: §8.4.1 · p. 132, §8.4 · p. 133

8.5 Entropic Regularization: Between OT and MMD

• Entropic regularization interpolates between Wasserstein and energy distances, and its sample-complexity rates interpolate between those of OT and MMD.The divergence tends to Wp as ε→0, while its pth power tends to the squared energy distance as ε→+∞.Source: §8.5 · p. 136

Variational Wasserstein Problems

• Wasserstein losses are used for texture synthesis, image segmentation, inverse problems, shape matching, and geometry-aware histogram factorization.Examples include image restoration and tomographic inversion; Wasserstein distances can replace Kullback–Leibler reconstruction loss when histogram-bin geometry is informative.Source: Variational Wasserstein Problems · pp. 137–138

9.1 Differentiating the Wasserstein Loss

• The Wasserstein loss is convex when parameters are histograms or simplex weights on fixed atoms, but general parameterizations are not necessarily convex.The section presents gradient approximation methods because variational Wasserstein problems require access to Wasserstein-distance gradients or approximations.Source: §9.1 · p. 138, Variational Wasserstein Problems · p. 138

9.1.1 Eulerian Discretization

• Eulerian discretization optimizes histogram weights on fixed locations, using entropic OT to obtain a differentiable objective and gradients for local gradient descent.The dual entropic OT potential supplies the histogram gradient; centering the potentials ensures mass conservation during gradient descent.Source: §9.1.1 · p. 140

9.1.2 Lagrangian Discretization

• Lagrangian discretization represents a measure by equally weighted moving particles, making the transport cost depend on particle locations and enabling gradients through the cost matrix.It is usually suitable for fitting continuous-valued observations in high-dimensional spaces, whereas Eulerian schemes suit low-dimensional supports or intrinsically discrete problems.Source: §9.1 · p. 139, §9.1.2 · p. 140

9.1.3 Automatic Differentiation

• When only a few Sinkhorn iterations are affordable, reverse-mode automatic differentiation of the algorithm’s output is preferable to using inexact dual-potential or coupling gradients.Differentiating Sinkhorn divergences costs about as much as evaluating the functional, but reverse-mode memory use grows with the number of iterations.Source: §9.1.3 · pp. 141–142

9.2 Wasserstein Barycenters, Clustering and Dictionary Learning

• Wasserstein barycenters extend means to probability measures; under squared Euclidean cost, one input with a density guarantees uniqueness, and the barycenter’s mean and support obey geometric constraints.For squared Euclidean cost, the barycenter’s mean is the barycenter of input means and its support lies in the convex hull of the input supports; empirical barycenters are consistent as sampled inputs increase.Source: §9.2 · pp. 143–145
• Barycenter fixed-point iterations converge under mild conditions, while a one-measure, squared-Euclidean barycenter constrained to at most k support points is equivalent to k-means.Entropic regularization smooths the resulting assignments, as in soft-clustering methods such as Gaussian mixtures.Source: §9.2 · p. 144

Special cases.

• For Gaussian inputs under squared Euclidean cost, the Wasserstein barycenter is itself Gaussian, and its mean equals the mean of the inputs.Source: Special cases · p. 146
• The Gaussian barycenter covariance is the unique positive-definite fixed point of a map, and fixed-point iterations converge in practice, with an alternative map also converging.The fixed-point characterization is studied in; the alternative map is proposed by Álvarez-Esteban et al. [2016].Source: Special cases · p. 146
• For 1-D distributions, monotone rearrangement makes the Wp barycenter nearly closed-form; discrete measures with fewer than n support points can be handled in O(n log(n)) operations.The computation uses a simple sorting procedure.Source: Special cases · p. 147
• Entropic barycenters can be computed by gradient descent, semi-dual descent, or generalized Sinkhorn iterations, and are used for 2-D and 3-D shape interpolation.The semi-dual approach accommodates additional regularization, such as smoothness; Sinkhorn iterations arise by rewriting the problem as a weighted KL projection.Source: Special cases · pp. 148, 150, Figure 9.3 · p. 150, Figure 9.4 · p. 151
• Wasserstein propagation generalizes barycenters from a star graph to distributions at unknown nodes, and entropy regularization extends Sinkhorn iterations to graph transport problems.Replacing the graph with a continuous domain yields measure-valued harmonic functions used in image and surface processing.Source: Special cases · pp. 151–152

9.3 Gradient Flows

• Explicit gradient steps require smoothness and can be unstable, whereas implicit steps handle nonsmooth functions and are stable; for convex F, they converge for any τ > 0.In the illustrated p-norm example, explicit steps are unstable for p = 1 and p = +∞, while implicit steps remain stable.Source: §9.3 · pp. 153, 155, Figure 9.7 · p. 155
• The Wasserstein gradient flow of negative entropy is linear heat diffusion, and JKO flows also model nonlinear evolution equations and constrained motion.Examples include porous-medium and total-variation flows, quantum drifts, heat evolution on manifolds, and congestion-constrained crowd motion.Source: §9.3 · p. 154
• A Lagrangian particle discretization evolves particle locations by coupled ODEs; with a quadratic potential and entropy, the limiting stationary distribution is Gaussian.For entropy, a nearest-neighbor distance estimator can provide a consistent density-based estimate from particles.Source: §9.3 · pp. 155–156
• Geodesic convexity ensures a unique continuous gradient flow that is the limit of discrete implicit steps and converges to a stationary distribution; entropy is a typical example.Linear functionals are geodesically convex when their potential is convex, and certain generalized entropy functionals are geodesically convex under stated conditions.Source: §9.3 · p. 157
• In gradient-flow discretizations, steps that are too small can cause locking, while convergence analysis is cited for a one-dimensional discretization method [Matthes and Osberger, 2014, 2017].The locking phenomenon arises from the polyhedral nature of the linear program defining the distance.Source: §9.3 · p. 158
• Unbalanced OT gradient flows allow mass creation or destruction, which is crucial for modeling physical, biological, or chemical phenomena, including Hele-Shaw cell growth.Such a flow can be approximated using the generalized Sinkhorn algorithm.Source: §9.3 · p. 158, §9.4 · p. 159

9.4 Minimum Kantorovich Estimators

• Maximum-likelihood estimation can fail for singular distributions or inaccessible densities, a typical issue for low-dimensional push-forward generative models that are easy to sample.Such models may have singular support and intractable densities, even though samples are generated by applying the parameterized map to reference samples.Source: §9.4 · pp. 159–160
• Weak metrics can replace the likelihood loss; evaluating an OT loss at fixed parameters requires semidiscrete transport, while parameter minimization is typically highly nonconvex.The continuous dual potential needed for a subgradient can be difficult to compute, especially in high dimensions; parametric approximations can incur unanalyzed errors and the formula is numerically unstable.Source: §9.4 · p. 160
• For OT loss, a primal optimal coupling gives a parameter-gradient formula that avoids computing gradients of the dual potentials.The alternative is based on computing an optimal coupling for an equivalent primal problem.Source: §9.4 · pp. 160–161
• Minimum Kantorovich estimators have been used for restricted Boltzmann machines, approximate Bayesian computation, GANs, and VAEs.Sinkhorn divergences are also used for shape matching and registration; OT learning pipelines may require adapting the ground cost, and OT can model and solve transfer or domain adaptation [Courty et al., 2017b,a].Source: §9.4 · pp. 161–162

10.1 Multimarginal Problems

• Multimarginal OT couples S histograms, extends naturally to measures, and supports entropic regularization with Sinkhorn scaling.Applications include quantum-chemistry approximations, particularly density functional theory [Cotar et al., 2013, Gori-Giorgi et al., 2009, Buttazzo et al., 2012]; for discrete inputs, the barycenter is also discrete.Source: §10.1 · pp. 163–165

10.2 Unbalanced Optimal Transport

• Unbalanced OT relaxes the initial Kantorovich formulation by penalizing marginal deviation using divergence Dϕ; τ1, τ2 control mass-variation penalties versus transport, and as equal τs tend to +∞, the balanced case recovers hard-marginal OT.The formalism recovers earlier formulations using an ℓ2 norm or an ℓ1 norm as in partial transport [Figalli, 2010, Caffarelli and McCann, 2010].Source: §10.2 · pp. 166–167
• Finite relaxation lets unbalanced OT vary mass during interpolation, avoiding bump-splitting and yielding smoother local matching than classical OT under hard mass conservation.The text notes that unbalanced OT can be useful for noisy or uncertain data, or when total mass carries meaningful signal; applications include image-label comparison [Frogner et al., 2015] and cell dynamics.Source: §10.2 · p. 168

10.3 Problems with Extra Constraints on the Couplings

• Capacity constraints upper-bound coupling density, while martingale transport constrains conditional means and requires convex order for feasibility.Entropic penalization and Dykstra’s algorithm provide approximate solutions; KL projection for the martingale constraint may require subiterations [Benamou et al., 2015; Guo and Obloj, 2017].Source: §10.3 · pp. 169–170

10.4 Sliced Wasserstein Distance and Barycenters

• A Lagrangian sliced-Wasserstein gradient flow converges numerically to the target with no local minimizer and yields a Monge-like matching, with applications to color transfer and texture synthesis.The method uses stochastic gradient descent over randomly sampled directions and is related to alternate minimization.Source: §10.4 · p. 171
• The Radon-transform approach computes sliced barycenters from separate one-dimensional barycenters and inverse reconstruction, but the result is not generally an exact Wasserstein barycenter.The text attributes the limitation to the Radon transform not being surjective, while noting that the approximation is numerically effective on fixed grids and often close in practice.Source: §10.4 · p. 172

10.5. Transporting Vectors and Matrices

• The sliced Wasserstein distance is Hilbertian, unlike W2 in dimensions d ≥2, enabling positive-semidefinite kernels and kernel-trick tasks including regression, classification, clustering, and dimensionality reduction.The passage gives parameter ranges of 1 ≤ p ≤ 2 for exponential kernels and 0 < p < 2 for energy-distance kernels.Source: §10.5 · p. 173
• For vector-valued measures, dual norms extend W1, the flat norm, and MMD; cone-valued OT covers positive measures and matrices, though static tensor formulations lack a sound theoretical framework.The cited tensor interpolation method generalizes Sinkhorn algorithms using quantum relative entropy.Source: §10.5 · pp. 174–175
• For positive trace-one matrices, the Bures metric generalizes Hellinger distance on diagonal matrices, matrix KL is convex in both arguments, and a Sinkhorn equivalent converges only in some cases.Convex dynamic formulations define OT distances between these matrices; the Sinkhorn equivalent seems empirically to always work.Source: §10.5 · p. 175

10.6.1 Hausdorff Distance

• Hausdorff distance defines a distance between compact sets and is bounded above by W∞ for measures with those sets as supports.The Hausdorff distance is connected to W∞ because the support of a measure coupling is a set coupling between the supports.Source: §10.6.1 · pp. 176–177

10.6.2 Gromov–Hausdorff distance

• Gromov–Hausdorff distance compares metric spaces by minimizing Hausdorff distance over isometric embeddings, but its optimization is highly nonconvex and globally intractable.Approximation schemes have supported shape-matching applications in vision and graphics [Mémoli and Sapiro, 2005, Bronstein et al., 2006]; smoothing the Hausdorff maximization motivates Gromov–Wasserstein distance.Source: §10.6.2 · pp. 177–178

10.6.3 Gromov–Wasserstein Distance

• Gromov–Wasserstein compares metric-measure spaces using within-space similarity matrices, but its nonconvex optimization is NP-hard in general.It defines a distance up to measure-preserving isometries and is tightly connected to Gromov–Hausdorff distance.Source: §10.6.3 · pp. 178–179

10.6.4 Entropic Regularization

• Entropic regularization enables successive linearizations solved by Sinkhorn iterations, interpreted as mirror descent, to reach a stationary point and compute soft domain maps.The approach was proposed initially by Gold and Rangarajan [1996, Rangarajan et al., 1999] and revisited for graphics by Solomon et al. [2016a].Source: §10.6.4 · p. 180, Figure 10.10 · p. 181