Breaking the Quadratic Space Barrier for Optimal Transport via Cost Matrix Compression

Optimal transport (OT) is a fundamental tool for comparing probability distributions, but the quadratic memory requirements of dense solvers limit scalability, especially on GPUs. We compress the ground-cost matrix using a sampling scheme inspired by Thorup–Zwick. We retain exact distances to nearby points and approximate the remaining distances through $\sqrt{n}$ sampled landmarks. For arbitrary ground metrics, this gives at most factor-$3$ distortion using $O(n^{3/2})$ expected space. More generally, a $k$-level scheme uses $O(k{n}^{1+1/k})$ expected space with distortion $2k-1$. We obtain two results using this compression scheme. First, we adapt the push–relabel framework of Lahn, Raghvendra, and Zhang to compute an $\epsilon$-additive approximation under the compressed costs using $O(n^{3/2})$ expected space, $O(\log n/{\epsilon }^2)$ parallel time, and $O(n^{3/2}/{\epsilon }^2)$ expected solver work. Under the original metric, the guarantee combines factor-$3$ distortion with $\epsilon$ additive error. In typical instances, especially when the input distributions are similar, much of the transport mass flows along local pairs whose costs are preserved exactly. Consequently, the compression yields high accuracy while substantially reducing storage. Our GPU implementation scales to approximately $600$k points on RTX 2060, compared to approximately $15$k for dense Sinkhorn implementations, with solution costs between $1.01$ and $1.35$ times the optimum on the evaluated high-dimensional datasets. Second, the compressed representation yields a $(15+\epsilon )$-approximation for incremental metric bipartite matching with $\tilde{O}(n)$ expected amortized insertion time and $O(n^{3/2})$ expected space after $O(n^2)$ preprocessing, for fixed $\epsilon >0$. This improves the superlinear update-time bound of Seth et al. (ICLR, 2026) for constant-factor approximation.