Optimal Transport: From Exact Linear Programs to Sinkhorn and Neural Potentials

The Feasibility-vs-Fidelity Tension

One tension drives the entire computational history of optimal transport: exact solutions are provably correct but computationally prohibitive, while tractable approximations each pay a distinct fidelity price. Every algorithmic generation from 1955 to 2024 is best understood as a specific answer to the question "how much fidelity can we afford to trade for the compute budget we actually have?" The five major lineages (classical LP, entropic/Sinkhorn, sliced projection, unbalanced/Gromov-Wasserstein extensions, and neural/flow-matching) converge, in recent work, on one unifying principle: recover a coupling or transport map that minimizes a cost, subject to the marginal and feasibility constraints your budget can enforce.

The table below places every major method on this trade-off surface. Methods appear as responses to the limitations of their predecessors; each section below follows that logic rather than a strict chronology. Textbook coverage is in the Mathematical Foundations section; algorithmic details are in Computational Methods and Extensions.

Villani's two volumes (Villani 2003, 2009) remain the definitive references; Santambrogio (Santambrogio 2014, 2015) is the most compact entry point with a calculus-of-variations emphasis; Ambrosio-Gigli-Savaré (Ambrosio et al. 2008) develops the gradient-flow theory; and Peyré-Cuturi (Peyré and Cuturi 2020) is the standard computational reference. Rachev's 1985 monograph (Rachev 1985) introduced the unified "Monge-Kantorovich" nomenclature to the Western probability community. Each reference carries characteristic strengths and well-known limitations^123456789, walked through method-by-method after the table.

The table groups the literature into five generations, each a response to the feasibility-vs-fidelity tension. The sub-sections below follow the same logic. The table below summarizes the structural trade-offs of each lineage, with the same caveat about its error column.

Generation 1 (1955-2000): Exact LP and Analytic Foundations

Classical exact-LP methods (Hungarian, network simplex) solve the discrete Kantorovich problem to numerical zero error, but the $O(n^3 \log n)$ cost makes them prohibitive beyond $n \approx 10^4$ source/target points¹⁰; this cubic wall is the principal motivation for every subsequent computational development. The same era produced the analytic theory: Brenier's (Brenier 1991) polar factorization established that, for absolutely continuous source measures and squared-Euclidean cost, the optimal map exists, is unique, and equals the gradient of a convex potential; Gangbo-McCann (Gangbo and McCann 1996) generalized this to any strictly convex superlinear cost and identified $c$-cyclical monotonicity as the universal characterization of optimal couplings. These foundational results give exact answers in closed form for tractable cases (Gaussians via Bures-Wasserstein, 1D via quantile composition) but provide no algorithm for the high-dimensional empirical regime that ML faces.

The fidelity of exact-LP is maximal (zero relative error) but the feasibility cost is a cubic barrier. Generation 2 resolves this by accepting a tunable bias in exchange for near-linear parallel complexity.

Generation 2 (2013-2018): Entropic Regularization and the Sinkhorn Revolution

Entropic-OT via the Sinkhorn iteration (Cuturi 2013) is the algorithmic pivot of modern computational OT: the addition of a small KL-regularizer $\varepsilon \mathrm{KL}(\gamma \| \mu \otimes \nu)$ makes the problem strictly convex, the optimal coupling factors as $P^* = \mathrm{diag}(\mathbf{u}) K \mathrm{diag}(\mathbf{v})$, and matrix-scaling iterations converge linearly. The complexity drops from $O(n^3 \log n)$ to $\tilde{O}(n^2/\varepsilon^2)$ (Dvurechensky et al. 2018), with massive parallelism across rows and columns, and the iterations are differentiable through the converged dual potentials; these three properties (speed, parallelism, differentiability) are jointly responsible for OT's penetration into deep learning. The price is the entropic bias: $W_\varepsilon \neq W$ and the entropic plan is "blurred" by an amount controlled by $\varepsilon$¹¹. Sinkhorn divergences $S_\varepsilon = W_\varepsilon(\mu,\nu) - \tfrac12 W_\varepsilon(\mu,\mu) - \tfrac12 W_\varepsilon(\nu,\nu)$ (Genevay et al. 2018) correct this debiasing at the cost of three solver runs; epsilon-scaling (geometric annealing of $\varepsilon$) recovers near-exact OT in roughly an order of magnitude fewer total iterations.

Entropic OT still scales as $O(n^2)$ in the cost matrix, which dominates memory and compute when $n$ is large and the ground cost is expensive to evaluate. Generation 3 sidesteps the cost matrix entirely by projecting to 1D.

Generation 3 (2018-2019): Sliced OT and Dimension-Free Projections

Sliced and projection-based methods (Kolouri et al. 2018; Le et al. 2019) sidestep the $O(n^2)$ cost matrix by averaging univariate Wasserstein distances over directions of the unit sphere $S^{d-1}$: $SW_p^p(\mu,\nu) = \int_{S^{d-1}} W_p^p(\theta_\#\mu, \theta_\#\nu)\,d\theta$. Each 1D projection sorts in $O(n \log n)$, so $L$ projections cost $O(Ln \log n)$; for typical $L = 50\text{-}1000$ this is two to three orders of magnitude faster than full OT. The dimension-free $O(n^{-1/2})$ sample complexity of sliced-Wasserstein (vs the $O(n^{-1/d})$ curse of full OT¹²) is its main statistical advantage and the reason it scales gracefully from 1D to 100D. The trade-off is projection bias¹³: SW is a true metric but only equivalent to $W_2$ for Gaussian distributions; on data concentrated on low-dimensional manifolds, projection averages over informative-and-uninformative directions equally and underestimates the true distance.

Both sliced and entropic methods still assume the source and target live in the same metric space and carry equal total mass. Many practical problems violate both assumptions. Generation 4 relaxes them.

Generation 4 (2017-2022): Unbalanced OT, Gromov-Wasserstein, and Neural Potentials

Unbalanced OT (Chizat et al. 2018) relaxes the equal-mass marginal constraint of classical OT; required whenever source and target have different totals (single-cell snapshots, partial-shape correspondences, class-imbalanced datasets)¹⁴. The relaxation replaces hard marginals with soft KL or TV penalties:

$$\begin{equation}\label{eq:uot-rw} \min_{\gamma \geq 0} \int c\,d\gamma + \tau_1 \mathrm{KL}(\gamma_1 \mid \mu) + \tau_2 \mathrm{KL}(\gamma_2 \mid \nu) \end{equation}$$

which admits a scaling iteration analogous to Sinkhorn. Gromov-Wasserstein (Mémoli 2011) handles the dual difficulty: source and target live in different metric spaces (graphs of different sizes, point clouds in different ambient dimensions), so there is no shared cost $c(x,y)$; GW instead matches pair-distances via a quadratic objective $\min_\gamma \sum_{ijkl} |d^X_{ik} - d^Y_{jl}|^2 \gamma_{ij}\gamma_{kl}$. The objective is a quadratic assignment problem, NP-hard in general and non-convex in the coupling; entropic GW solvers depend strongly on initialization, with random / spectral / identity warm-starts producing different local optima¹⁵.

The unbalanced and GW extensions are not Pareto-dominated by classical OT; they enable problems (unequal-mass biology, graph matching) that classical OT cannot even formulate. But all of the above methods require either a discrete cost matrix or a fixed-support approximation, blocking application to the continuous high-dimensional distributions that arise in generative modeling. Generation 5 abandons the matrix entirely.

Neural-OT methods (Korotin et al. 2023; Seguy et al. 2018; Makkuva et al. 2020) parameterize Kantorovich potentials $f, g$ (or the Brenier potential $\phi$) as neural networks and solve the dual saddle-point min-max objective by stochastic gradient ascent-descent. The ICNN saddle-point solver is due to Makkuva et al. (Makkuva et al. 2020); the architecture it builds on, the input convex neural network (ICNN) of Amos, Xu and Kolter (Amos et al. 2017), was introduced for structured prediction and control and contains no transport objective of its own. ICNNs enforce the convexity required by Brenier's theorem via non-negative weights and convex non-decreasing activations; this guarantees the gradient $T(x) = \nabla \phi_\theta(x)$ is a valid Brenier map. The price is the adversarial training instability of all GAN-class objectives¹⁶: hyperparameter sensitivity, mode collapse, and oscillation. These are mitigated by careful learning-rate scheduling, gradient penalties, and warm-starting from Sinkhorn-derived initial maps. JKOnet (Bunne et al. 2022) uses ICNNs to discretize Wasserstein gradient flows via the JKO scheme, learning energy landscapes from population snapshots.

Generation 5 (2022-2024): Flow-Matching, Schrödinger Bridges, and Unification

OT flow-matching and Schrödinger-bridge methods are the 2022-2024 generative-modelling synthesis. Flow-matching trains a velocity field whose ODE integration transports noise to data; the OT-conditional variant uses mini-batch OT couplings (rather than the diffusion's stochastic noise schedule) to define straighter probability paths, which simpler velocity fields can learn and which integrate in fewer ODE solver steps¹⁷. Schrödinger bridges solve the entropy-regularized OT problem in path-measure space via iterative proportional fitting (the dynamic analog of Sinkhorn), interpolating between deterministic OT ($\varepsilon \to 0$) and pure diffusion ($\varepsilon \to \infty$).

The five lineages (exact-LP, entropic, sliced, unbalanced/GW, and neural/flow) converge on one principle: recover a coupling or transport map that minimizes a transport cost under the marginal and feasibility constraints your compute budget can enforce. Brenier potentials, Schrödinger bridges, and continuous-time normalizing flows are different parameterizations of the same underlying displacement-interpolation theory introduced by Brenier (Brenier 1991) and Benamou-Brenier (Benamou and Brenier 2000) three decades earlier. The figure below maps all five lineages onto a feasibility-vs-fidelity trade-off plane.

References