Source-linked AI summary
Sinkhorn Distances: Lightspeed Computation of Optimal Transportation Distances
Marco Cuturi
TL;DR
Optimal transportation distances are costly to compute for histograms with dimensions above a few hundred. The paper introduces entropic regularization solved by Sinkhorn-Knopp scaling, achieving over 100,000× speedups over EMD solvers and better MNIST performance than competing distances.
Problem
Computing optimal transportation distances becomes costly for histograms with dimensions above a few hundred, limiting their use in large-scale data analysis.
Method
The paper adds an entropic regularization term to the optimal transportation problem and solves the resulting strictly convex problem with Sinkhorn-Knopp matrix scaling.
Results
More than 100,000× faster than EMD solvers on large dimensions, Sinkhorn distances also beat all evaluated alternatives, including EMD, on MNIST classification.
Takeaways & Limitations
Entropic regularization yields Sinkhorn distances that combine substantially faster computation with competitive or improved performance in the reported applications.
Takeaways & Limitations
For small α, Sinkhorn distances can violate the coincidence axiom because dM,α(r, r) > 0 for some r.
Abstract
from arXiv · showhide
Optimal transportation distances are a fundamental family of parameterized distances for histograms. Despite their appealing theoretical properties, excellent performance in retrieval tasks and intuitive formulation, their computation involves the resolution of a linear program whose cost is prohibitive whenever the histograms' dimension exceeds a few hundreds. We propose in this work a new family of optimal transportation distances that look at transportation problems from a maximum-entropy perspective. We smooth the classical optimal transportation problem with an entropic regularization term, and show that the resulting optimum is also a distance which can be computed through Sinkhorn-Knopp's matrix scaling algorithm at a speed that is several orders of magnitude faster than that of transportation solvers. We also report improved performance over classical optimal transportation distances on the MNIST benchmark problem.
1. Introduction
Optimal transportation distances are parameterized and effective for histogram analysis, but their computational cost limits large-scale use. The paper addresses this bottleneck by introducing entropic regularization, yielding strictly convex optimization and faster computation.
- Motivation: Optimal transportation distances are uniquely parameterized by a ground metric, which helps handle high-dimensional histograms.They are contrasted with Hellinger, χ2, Kullback-Leibler, and Total Variation distances, which are not parameterized.
- Motivation: Histograms and bags-of-features are central to data analysis, while optimal transportation distances perform well in practice, especially for comparing images.This has motivated substantial theoretical and practical interest.
- Computational bottleneck: More than a few seconds may be required to compute one distance between histograms with dimensions in the few hundreds, hindering large-scale data analysis and machine learning.The computational burden can undermine the relevance of optimal transportation distances in machine learning.
- Proposed approach: The paper regularizes optimal transportation with an entropic term motivated by the maximum-entropy principle.The authors argue that this regularization is intuitive given the geometry of the optimal transportation problem and established in transportation theory.
- Proposed approach: The entropic regularization turns the linear program into a strictly convex problem.The paper presents this transformation as one of the regularization’s optimization virtues.
2. Reminders on Optimal Transportation
Optimal transportation maps histogram r to c by minimizing the cost over matrices with prescribed row and column sums. When M is a metric matrix, the optimum defines a distance, but computing it for general M has worst-case complexity O(d^3 log d).
- Transportation polytope: U(r, c) is the transportation polytope containing nonnegative d × d matrices whose row and column sums are r and c.For multinomial variables with distributions r and c, each P ∈ U(r, c) represents a possible joint probability table.
- Optimal transportation: Given cost matrix M, a transportation matrix P incurs cost ⟨P, M⟩, and optimal transportation minimizes this quantity over U(r, c).The optimum is denoted d_M(r, c), and an optimal table P⋆ can be obtained with network simplex or other approaches.
- Metric conditions: d_M(r, c) is a distance when M is a metric matrix satisfying zero diagonal and the triangle inequality.The metric-matrix conditions are m_ii = 0 and m_ij ≤ m_ik + m_kj for all applicable indices.
- Computational limitation: O(d^3 log d) is the worst-case complexity for computing the optimum with known algorithms for general M, and practical complexity is super-cubic.Restrictions on M can enable much faster approximate solutions, but may reduce performance and applicability.
3. Sinkhorn Distances
Sinkhorn distances restrict transportation plans by an entropy or KL-divergence constraint around the independence table, yielding a max-entropy alternative to classical optimal transport. They remain symmetric and satisfy triangle inequalities for every α, with classical transport and an independence-kernel form as limiting cases.
- Definition: The feasible set U_α(r, c) comprises transportation tables whose KL divergence from the independence table rc^T is at most α, equivalently whose entropy is sufficiently high.These tables also represent joint probabilities with sufficiently small mutual information.
- Motivation: Entropic regularization favors plausible, smooth low-cost transportation plans over the quasi-deterministic extreme plans selected by classical linear optimization.Classical optima occur at vertices with at most 2d − 1 non-zero elements.
- Limiting cases: For α large enough, the Sinkhorn distance d_M,α equals the classical transportation distance d_M.The entropy-constrained feasible set then becomes the full transportation polytope.
- Limiting cases: When α = 0 and M is a Euclidean distance matrix, d_M,0 = r^T M c and defines the negative definite independence kernel.Specifically, e^(−t r^T M c) is positive definite for every t > 0.
- Metric properties: For every α ≥ 0 and M ∈ M, d_M,α is symmetric and satisfies all triangle inequalities, while coincidence requires multiplying by 1_{r≠c}.For small enough α, d_M,α(r, r) > 0 when h(r) > 0, so the unmodified function may fail the coincidence axiom.
4. Computing Sinkhorn Distances with the Sinkhorn-Knopp Algorithm
The entropy-regularized transportation problem is linked by duality to a multiplier-based divergence whose optimizer has a matrix-scaling form. Sinkhorn-Knopp iterations compute this optimizer efficiently, while bisection selects the multiplier for the desired entropy level; the paper thereafter calls the cheaper dual-Sinkhorn divergence the Sinkhorn distance despite lacking a proof that it is a distance.
- Dual formulation: For every α and histogram pair (r, c), duality associates a λ ∈ [0, ∞] such that d_M,α(r,c) = d_M,λ.The dual-Sinkhorn divergence can be computed at much cheaper cost than classical optimal transportation for reasonable λ.
- Optimizer structure: For λ > 0, strict convexity makes P^λ unique, with entries of the form u_i e^(-λm_ij) v_j for uniquely determined scaling vectors up to a multiplicative factor.This form is the matrix-scaling characterization used to solve the regularized problem.
- Sinkhorn-Knopp algorithm: Given e^(-λM) and marginals r and c, Sinkhorn-Knopp fixed-point iterations converge to P^λ while enforcing the prescribed row and column sums.The vectorized implementation sets K = exp(-λM), iteratively rescales using r and c, and handles null coordinates by retaining strictly positive entries of r.
- Multiplier selection: Bisection finds λ by increasing it until the solution entropy reaches h(r) + h(c) − α, using the fact that entropy decreases monotonically with λ.The procedure iteratively computes P^λ and adjusts λ to meet the entropy target.
- Scope and limitation: The paper uses the cheaper dual-Sinkhorn divergence thereafter and calls it the Sinkhorn distance, although it is not provably a distance.The authors leave exact computation of d_M,α for future work.
5. Experimental Results
Experiments on MNIST evaluate Sinkhorn distances under cross-validated settings and compare their accuracy, convergence toward EMD, and computational speed. Sinkhorn distances achieve comparable classification evaluation while offering dramatically faster computation, especially on GPUs, although practical regularization leaves them about 10% above EMD.
- MNIST setup: MNIST digits3 are represented as normalized intensity histograms on a 20 × 20 grid, with training subsets ranging from 3 × 10^3 to 25 × 10^3 points.The ground metric is the Euclidean distance between grid pixels, producing a 400 × 400 distance matrix.
- Evaluation protocol: Classification errors are reported as means and standard deviations over 24 experiments using 4-fold cross-validation repeated 6 times.Kernel parameters are selected separately for each distance and training fold from {1, q10(d), q20(d), q50(d)}.
- Parameter selection: λ = 9 is selected for Sinkhorn distances by internal cross-validation from {5, 7, 9, 11} × 1/q50(M), balancing diagonal dominance against closeness to classical optimal transportation.Smaller values of the Independence-kernel parameter a are preferred.
- Convergence to EMD: 10%: even at large λ, Sinkhorn distances remain about 10% above EMD on 402 distinct MNIST pairs.The relative gap decreases as λ increases, but numerical closeness to EMD is not expected or considered desirable for practical λ.
- Computational speed: 100.000×: on large dimensions and the same CPU, Sinkhorn distances are more than 100.000 faster than EMD solvers at threshold 0.01; GPUs add another order-of-magnitude speed-up.Rubner et al.’s implementation cannot run for histograms larger than d = 512.
- Computational speed: Iterations increase as e−λM becomes diagonally dominant, so a fixed iteration count depending only on λ is recommended for parallel computation and can yield still larger GPGPU speedups.The convergence tolerance is 0.01 in the norm of successive iterates of x ∈ R^d.
6. Conclusion
The paper shows that entropic regularization of optimal transport enables effective speed-ups independent of the ground metric’s structure. Preliminary evidence suggests Sinkhorn distances perform at least as well as EMD and may perform better in applications.
- 6. Conclusion: Entropic regularization guarantees effective speed-ups regardless of the structure of the ground metric M and suggests Sinkhorn distances may outperform EMD in applications.The authors identify new research directions and potential applications at the intersection of optimal transportation theory and machine learning.
7. Appendix: Proofs
The appendix proves that the independence-related expression defines a negative-definite kernel and establishes the entropy inequality underlying Lemma 1. It also gives a preprocessing strategy and specifies the domain where the independence kernel is positive definite.
- Proof of Property 1: r^T Mc is negative definite because it decomposes into a sum of negative-definite kernels.The decomposition uses r^T u + c^T u and -2r^T Kc, with K positive definite.
- Implementation remark: Precomputing the norm vector and a Cholesky factor enables faster computation of the independence kernel.The suggested preprocessing stores Lr and precomputes the diagonal term r^T u for each histogram.
- Implementation remark: The independence kernel is positive definite on histograms with the same 1-norm but not on arbitrary vectors.This condition limits the positive-definiteness guarantee to equal-mass histograms.
- Proof of Lemma 1: The constructed matrix S is a transportation matrix between x and z, and its entropy satisfies h(S) ≥ h(x)+h(z)−α.The proof represents the construction through a joint distribution T over three variables.
- Proof of Lemma 1: The construction makes X → Y → Z a Markov chain, so data processing yields I(X;Y) ≥ I(X;Z) and the stated entropy inequality.The factorization p(X,Y,Z) = p(X,Y)p(Y,Z)/p(Y) establishes the Markov-chain structure.