A Scalable Constant-Factor Approximation Algorithm for $W_p$ Optimal Transport
Authors: Pankaj K Agarwal
Oliver Chubet
Sharath Raghvendra
Keegan Yao
Venue: ICLR 2026 Poster (2026)
Keywords: optimal transport
p-Wasserstein
min-cost matching
Primary area: optimization
TL;DR: We give a randomized $(4+\varepsilon)$-approximate $W_p$-optimal transport plan in $O(n^2+\varepsilon^{-1}n^{3/2}\log\Delta)$, where $\Delta$ is the ratio of the diameter to smallest distance in the input.

Abstract: Let $(X,d)$ be a metric space and let $\mu,\nu$ be discrete probability distributions supported on finite point sets $A,B \subseteq X$.
For any $p \in [1,\infty]$, the {\it $W_p$-distance} between $\mu$ and $\nu$, $W_p(\mu, \nu)$, is defined as the $p$-th root of the minimum cost of transporting all the probability mass from $\mu$ to $\nu$, where moving a probability mass of $\delta$ from $a \in A$ to $b \in B$ incurs a cost of $\delta d(a,b)^p$.
We give a (Las Vegas) randomized algorithm that computes a $(4+\varepsilon)$-approximate $W_p$ optimal-transport (OT) plan in $O(n^2 + (n^{3/2}\varepsilon^{-1}\log n\log\Delta)^{1+o(1)}\log U)$ time with probability at least $1-1/n$, for all $p \in [1,\infty]$, where $\varepsilon > 0$ is an arbitrarily small constant and $\Delta$ is the ratio between the largest and smallest interpoint distances in $A\cup B$.
The previous best result achieved an $O(\log n)$-approximation in $O(pn^2)$ time, for constant values of $p$.
Our algorithm significantly improves the approximation factor and, importantly, is the first quadratic-time method that extends to the $W_\infty$-distance.
In contrast, additive approximation methods such as Sinkhorn are efficient only for constant $p$ and fail to handle $p=\infty$. \changed{Our algorithm also extends to a query model where, for any integer $k > 1$, we give an algorithm that preprocesses $X$ into clusters in $O(n^2+kn^{1+1/k}\log n\log\Delta)$ time, after which a $O(k)$-approximate $W_p$ distance between any two distributions $\mu$ and $\nu$ with $X$ as support can be computed in $(n^{1+1/k}\log n\log\Delta)^{1+o(1)}$ time with probability at most $1-1/n$.}
Finally, for $p=\infty$, we show that obtaining a relative approximation factor better than $2$ in $O(n^2)$ time would resolve the long-standing open problem of computing a perfect matching in an arbitrary bipartite graph in quadratic time.

Published as a conference paper at ICLR 2026
A SCALABLE CONSTANT-FACTOR APPROXIMATION ALGORITHM FOR W, OPTIMAL TRANSPORT
Pankaj K. Agarwal",, Oliver Chubet², Sharath Raghvendra', and Keegan Yao
Duke University, North Carolina State University
ABSTRACT
Let (X, d) be a metric space and μ., v be discrete probability distributions sup- ported on finite point sets A, B C X. For any E [1, 20). the Wp-distance between μ and v, W,(,, v). defined as the p-th root of the minimum cost transporting all the probability mass from μ to D, where moving probability mass of from a E to E incurs a cost of ōd(a, b)p. We give a (Las Vegas) randomized algorithm that computes (4 + aaappooxime W, optimal- (n3 g- log log A) A)++0(1) log U) time with probability at least 1/n, for all p E [1, co). where e > 0 is an arbitrarily small con- stant and A is the ratio between the largest and smallest interpoint distances in A AUB. B. The previous best result achieved an O(log r1) -ppaappooinntiion in (pp=²) time, for constant values of p. Our algorithm significantly improves the approxi- mation factor and, importantly, is the first quadratic-time method that extends to the -istaaneece. In contrast, additive approximation methods such as Sinkhorn are efficient only for constant p and fail to handle p cc. Our algorithm also extends to a query model where, for any integer > 1, we give an algorithm that preprocesses X into clusters in (n kn 1/1 log n log A) time, after which a (ap))-appoxximat W, distance between any two distributions μ and with X support can be computed in (n +1// log n log A) time with probability most 1 1/nnn. Finally, for = oc, we show that obtaining a relative approximation factor better than 2 in O(n²) time would resolve the long-standing open problem of computing a perfect matching in an arbitrary bipartite graph in quadratic time.
(OT) plan in On
INTRODUCTION
Let μ and v be discrete probability distributions supported on sets A and B, respectively, with [A] + |B] = For each pair ba E A x B, let da, b) denote the ground distance between a and transport plan a function
B > Rzo that assigns a mass to each pair (a, b) such that EsEB (a, b) < p(a) and E.eA oa, b) I v(b). Given parameter p 1, suppose the cost of moving probability mass of 5 from point E A to point b E B is given by d(a, b)". The W, cost of any transport plan a between μ and v is defined as
1/1
uy(o)= =
Σ BEALEB
aaa, b) d(a.b))
The W, cost above for finite p extends naturally to the Wa cost of a transport plan o, defined
Woo (a) lim wp(o)= a. AxB (a.
d(a, b).
In the W, optimal transport (OT) problem, we wish to compute the transport plan o² that transports the entire mass and has the smallest W, cost. We refer to the cost u,'(a*) as the Wy-distance and denote it by W,(,,v) If μ and are discrete distributions, each supported n points with every
FFollowing convention from Theoretical Computer Science, all authors ordered alphahetically.

---

Published as a conference paper at ICLR 2026
point assigned a mass of 1/n, then the corresponding Wp optimal transport problem is called the Wp matching problem.
The W distances, for varying values of p. possess several appealing properties that make them favorable in many applications. OT plans under the W1 distance measure total displacement and are therefore useful to capture structural properties such semantic relationships from word embed- dings (Kusner et al. (2013). OT plans arising from the W2 distance have nice structural qualities such monotonicity (Brenier (1991): Aurenhammer et al. (1998)) and translation invariance (Co- hen & Guibas (999)) , and tend to preserve the geometry of the distributions. Furthermore, recent work in machine learning and topological data analysis uses Wo distance to establish consistency and convergence properties of topological summaries (Vishwanath et al. (2020): Damrich et al. (2024), ) and to design topological layers in neural networks (Kim et al. (2020p) Due to these fa- vorable properties, Wr distances have been used in applications across machine learning (Chang et al. (2033): CCuang et al. (2022)), computer vision (Backurs al. (2020): ?), and natural language
processing (Alvarez-Melis & Jaakkola (2018) Yurochkin et al. (2019)).
From an algorithmic standpoint, the exact computation of the W, distance between discrete distributions can be formulated as a minimum-cost flow (MCF) problem, which can be solved in (+o((1) time using recent advances in MCF algorithms (Chen et al. (2022)). While these results mark important theoretical progress, they are highly complicated, making them unsuitable for prac- tical implementations. Indeed, even the simpler task of designing truly quadratic-time exact algo- rithm for deciding whether a dense graph admits a perfect matching remains a longstanding open problem in graph theory (Behnezhad et al. 2024)).
Because practical exact algorithms remain expensive, research has shifted toward scalable ap- proximation methods. In this setting, an algorithm is said to compute relative a-approximation for Wp(n, 2) it outputs a transport plan a whose cost satisfies (p(o) S R. (1,, v). A seminal result by] Charikar (2002) introduced an Olog apprrooximaaionn for Wi by embedding the ground metric into a hierarchically well-separated tree; a greedy transport procedure on the tree yielded an exact solution in (n²) time, producing an overall Olog n) appaappoxximation. This work opened the door to more refined methods, and subsequent efforts developed near-linear-time (1 >-aapproximation algorithms under additional assumptions, such as when the ground distance is Euclidean in fixed dimensions (Agarwal et al. (2022) 2202); & (2023))) and more recently, sub-quadratic algo- rithms for higher-dimensional Fuclidean settings (Andoni & Zhang (2023); |Beretta et al.] (2025). )- However, none of these techniques extend naturally to the case p ' 2. Building on this line of work, Lahn et al. (2025) recently presented a relative Ollgg n) apprrooxxiiiaooon algorithm for any finite p 2 and for any metric, with runtime O (n² log log A log n), where log U is the bit-length of the input probabilities and A is the spread of A U B (the ratio of its largest to smallest nonzero pairwise distance).
One influential direction of work was introduced by Cuturi regularization of OT. guarantees solutions within an additive error of & diam(A U B), where diam(A (A U B) denotes the maximum ground distance between points in A UJ B. Although weaker than a relative (1 + -)-aprroximation it applies across all metrics and inspired a series of additive approximation algorithms, including parallelizable variants (Altschuler et al. (2017); Dvurechensky et al. y); Jambulapati et al.] (2019): > Lahn et al. (2019, 2023)). Nonetheless, their runtimes re- main on the order of n² n01) for W1 and worsen to n² // Otp) for larger p. with no extension to the case p 00.
2013) who proposed entropic
Our results. In this paper, we present two constant-factor approximation algorithms for the W problem. These are the first truly quadratic-time (assuming A + U 20(n)(*) approximation algorithms for the WP problem over any ground metric, applicable to all p E [1, oc]: that is, including p co. This improves the Oloo n) approximationn of Lahn et al. (2025) constant factor, while extending the guarantee to every p.
Let I and v be two discrete distributions supported on a set of n. points in an arbitrary metric space. Let p E [1) be a parameter, and let E >> be an arbitrar- ily small constant. A (4 + -ppproxmmttt OT plan under the W, metric can be computed in Onn² + -n²/>=-1 log A logn) +o(1) log U) time with probability at least 1 - 1/n, where A is the spread of the support set and U is the ratio of the maximun to minimum (non-zero) probability in μ or
Theorem 1.1.
2

---

Published as a conference paper at ICLR 2026
For any integer k >> 1, we construct a directed spanner of size O log A) (for k = 2), and more generally OKk log A) for arbitrary k. We prove that the shortest-path distance from any b E B to any a E A in the spanner approximates d(a, b). within a multiplicative factor of (4 + =)P for k = 2, and (2 + =))² for general k. Running a minimum-cost flow algorithm this directed spanner then yields an approximate optimal transport plan; in partieular, instantiating it with the algorithm of ohen et al. (2022) and setting A = 2 implies Theorem [.]
Our main technical contribution is a new technique for approximating .(.,.),, drawing on ideas from Bourgain's 's multi-level sampling [ourgain (1985). Multi-level sampling has inspired numer- ous clustering- based constructions, including metric spanners and distance oracles. Classical span- ner constructions Har-Peled et al. (2023): Bassana & Sen (2007); Cohen (1998) require super- quadratic preprocessing time and approximate only the underlying metric. However, for p > 1, the transport cost d(. which is not a metric (it fails the triangle inequality), and therefore cannot be handled by these undirected metric spanners. Our spanner differs from prior construc- tions in two key ways. First, we introduce carefully chosen Steiner points at multiple scales of the clustering, which substantially simplify the structure of the spanner and reduce its construction time to O(n² + kn +/1 log A) which is a significant improvement over classical spanners. Second, to approximate d(, we assign a direction to each edge, producing a directed spanner. Although our directed spanner need not be strongly connected, it guarantees that there exists a path from any b E B to any a E A. We show that it approximates the cost d(a, b). from any b E B to any a E A within a factor of (4 +++ =)P (when k = 2) and (2k(1 + =))² (for arbitrary k).
Query-model interpretation. Theorem 1.1 naturally extends to a query model, where the support set X is fixed in advance and many distributions supported on X must be compared. This set- ting arises frequently in natural language processing. For example, a document can be represented as a probability distribution over fixed word embeddings (Kusner et al., (2015)); the similarity be- tween two documents is measured by the W, distance between these distributions, referred to as the Word Mover's Distance (WMD). Kusner et al. (2015) nevertheless note that the quadratic (or worse) running time of optimal transport is a major barrier to the scalability of WMD while process- ing large set of documents. In contrast, for any integer k > 1, our algorithm performs a one-time O(n2 1 log A) preprocessing of X, after which each W, (or >o query between any dis- tributions supported on X can be answered by using any MCF solver on a sparse directed graph of size Okk 1+1/ log A) for an (k)) apprrxpppxiiittion. Thus, in applications like WMD where many OT computations share the same support, our method yields sub-quadratic query time.
Since the interior-poin point method of Chen et al. (2022) does not admit a practical implementa- tion. also develop simple combinatorial algorithm for optimal transport. A crucial ingredi- ent in implementing combinatorial algorithms efficiently is the ability to answer dynamic weighted bichromatic closest-pair (BCP) queries quickly. Existing work has used Bourgain's 's sampling-basee clusterings to construct distance oracles, which return approximate point-to-point distances (Tho-] rup & Zwick (205)); Mendel & Naor (2007): Awerbuch al. (1998). However, while such oracles are well suited for distance estimation, it appears extremely challenging to adapt them to support dynamic BCP queries, which require repeatedly identifying the minimum over all cross pairs un- der changing weights. To address this gap. we tailor our clustering scheme so that it naturally supports dynamic weighted BCP queries. Importantly, although pairwise-distance computations are expensive in our framework, our algorithms avoid them entirely and rely solely on the BCP primitive. Using these data structures with standard implementation of combinatorial OT algo- rithms (Vaidya 1989)), we obtain a simple combinatorial algorithm that computes, for any integer k >>> 1 an apprrxxaaooiiioo to the Wp optimal-transpon cost in O(n log A log U log" n) time (Theorem 1.2]. as well as a quadratic time (4 + appppooxiee matching algorithm (Theo- rem I.3 complementing the quadratic OT algorithm obtained through the interior-pois method of Chen et al. (2022) (Theorem I.I).
Theorem 1.2. Let (X, d) be an arbitrary metric space, let and v be two probability measures supported on finite point sets A, C X with A U B| n, and let p E [1, oo] be a parameter For any integer k > 1, an apprrxxititiino to the Wp optimal-transport cost W v) can be computed in ( log" log A logg U) time with probability at least 1- 1/n, where A is the spread of B and log U is the maximum to minimum (non zero) probubility at μ and v.
Theorem 1.3. Let A and B be two point sets of size n each in an arbitrary metric space, and let p (1, oo) be parameter. A (4 + -ppproximmaee W, matching of A and B can be computed
3

---

Published as a conference paper at ICLR 2026
in O 3/2 log" log A) time with probability at least 1 - 1/n, and (k)) ((K))apppoximatee Wp matching of A and B can be computed in On +n3/241/kg------ 3/2 log" n log A) time with probability at least 1 1/n.
To our knowledge, Theorems 1.2 and 1.3 provide the first practical and implementable approxi- mation algorithm for computing the W oc: OT in sub-n? time and -matching in quadratic time. We also establish conditional lower bounds that suggest our results cannot be significantly improved without major breakthrough in the graph matching problem, namely, computing perfect matching in any bipartite graph in (n²) time.
Theorem 1.4. there exists a quadratic-tine algorithm that computes a rreative approxima- tion or diamm/2- additive approximation of the W Wo--matt -mathhing problem, where diam is the diameter of the point set and e >> 0 is any arbitrarily small constant, then a perfect matching in an arbitrary hipartite graph can be computed in On?)) time if one exists.
We conclude with a primitive implementation of our simple combinatorial algorithm alongside some experimental results suggesting that the algorithm computes good quality matthhnhin for p E [1, oo] in Section 4 While we prove that the approximation factor of our algorithm is (4 + =) in the worst case, our experimental results indicate that our algorithm computes even better approx- imations of wmtmhhhhnngs in practice.
2 DISTANCE APPROXIMATION AND PROXIMITY QUERIES
Let P be set of points equipped with metric d: P x P + Rz0 and let k > 1 be integer. We introduce a k level clustering-base basedd distance function that approximates d(:,.) within a multi- plicative factor of 21(11 E). Our construction is similar in spirit to the construction of k-spanners and distance oracles discussed in Section | Our distance can be represented using roughly kn k//+1/K space (as opposed to ((n²) space to store all pairwise distances), and is used both to build a span- and to maintain weighted bichromatic closest pairs. For simplicity in presentation, we describe our clustering construction and associated algorithms for k = 2 in the paper. The extension to any arbitrary integer & > 2 deferred to the Appendix
We begin with a few notations. Given a point a E P and a subset Q C P, the distance from x to Q is defined as d(x, Q) mingeq dzz, q). For point q E P and subset Q C P. define the Voronoi set of q to be
Vga. Q) = ty E P d(y. g) < d(y, Q)}.
That is, V(g, Q) consists of the points in P for which g is closer than any point in Q.
Two-level clustering. We construct a two-level clustering of points of P. Set P. = P. Next, we choose a subset P1 S Po by sampling each point in Po independently with probability n - The expected size of P1 is E[] = 1/2
Let A be the spread of the point set P. Without loss of generality, assume minp. d d(p. 9) implying that the diameter ma d(p.g) is A. We choose & 0 be a sufficiently small
constant. Set t [log }) A1, ro 0, and r. = (1+ + 4)' for S i i t. We generate two types of clusters: (i) For each E Po \ Pi and for every < f, define C = fo E V(q, (q, P1) d(z, q) S r.). (ii) For each 9 E P1 and for every < f, define Cg = {o E Po d(z,g) I re). We refer to i as the index of the cluster [.) Let C = [ g Po. S 1) be the collection of all clusters. Note that point p P may belong to many clusters. The number of clusters that contain p is called the degree of p and is denoted as (p). While the degree of any particular point may be (n) in the worst case, we prove that, with high probability, the degree of all points in P is much smaller.
Lemma 2.1. For all p E P. E (p)) O(ynE- log A). Furthermore, there is constant > 0 such that Pr[Vp E P, degc(p) S cynE= log n log A) 2 1-1/n.
Proof. We first bound the number of clusters of type (ii). Each point of Po is selected independently
with probability 1/v so E[[P|]] Vn. By Chernoff bound, [P]] I 2/77 with probability at least +. Each point E P1 generates O( log A) clusters of type (ii). and any point may participate in each such cluster. Therefore, the expected number of type (ii) clusters that p
4

---

Published as a conference paper at ICLR 2026
participates in is O(ynE- loge" log A). and with probability at least 1/n2, the number of type (ii) clusters that participates in is bounded by O(yn=" loo A).
We now bound the contribution from clusters of type (i). For each q E Po \, PI. we compute its Voronoi set V (g, P1) and then construct O(= log A) clusters from that region. therefore suffices to bound the number of Voronoi sets that contain fixed point p E P. Fix p E P. Let wi...... W be the points of Po ordered in non-decreasing distance from p. Let be the smallest index such that w, E PI.
For any index >> 8, we have d(p, wx) > d(p w.). Since w, E Pt, the point p is closer to w, than to we, and hence p & V(we, Pi). Thus, p can belong V(uz, Pi) only if f 8, and therefore the number of Voronoi sets containing p is at most s. Since each point is selected into P1 independently with probability 1/v the random variable is geometric with success probability 1/ Hence, the expected value of is E(s] = vr. Moreover,
Vnl J/)
Pr(s > t = 1
S Se-/aa
Choosing 2/n log n, get Pr(s 2//log n| I c c2 = Therefore, the number
of type (i) clusters containing p is O(se- log A), which is O(vne- log A) in expectation and O(yn log e log A) with probability least 1/n2. Combining the bounds for type (i) and type (ii) clusters and applying union bound completes the proof.
Our algorithms and their running-time bounds assume that every vertex E P satisfies (9) O(vnE log n log A). Since this bound holds with probability least - 1/n, the stated exccution-time guarantees also hold with the same probability. Also, because each point par- ticipates in O log n log A) clusters, the total space required to explicitly store all clusters O -- log A) in expectation and (n e- log log A) with probability at least 1/n.
Constructing clusters. Constructing clusters of type (i) proceeds as follows. For each point q Po Pt, we first compute closest point in P1 by scanning P1. The Voronoi set V(4, Pi) is then determined by checking, for every q' E Po. whether o' is closer to g than to its nearest neighbor in Pt. This requires O(n) time per q. Once V(q. Pi) is obtained, the associated clusters are formed by partitioning its points according to their distance from q. which can be done with single additional scan. Thus, clusters of type (i) can be constructed in (n) time per E Po \ Pt. and O(n²) time overall. For clusters of type (ii), each point q E P1 induces partition of Po based distance from q. which takes (n) per g and (n²2) time in total.
Cluster-induced distance approximation. Using this two-level clustering, we define dis- tance function de: P x P > Rzo that approximates the ground metric within a factor of (4 + e). For any pair of points x. y E P, let i be the smallest index of cluster that contains both and y- Then we set de (2. y) = 2r,.
Lemma 2.2. d(x, y) S de v) s (4 (4+ c)d(z,))
Proof. We say a point x is separated from y by P1 if there is a point a E P1 such that d(ær, a) <
d(x y)- From this and triangle inequality, d(y, a) < d(xr, a) + d(x, y) I 2d(x, y). Let be the
smallest index such that E [.]. Then, d(y, a) S < (1 + {d a). So we have de(x, v) S 2r; < 2 (1 + &) d(y,a) < 4 (1 + 5) d(x, If is not separated from y. then E V(,, Pi). Let
be the smallest index such that Cy contains x. Then, dz v) 2r; S 2 (1 d(r
LI
2.1 PROXIMITY QUERIES
Next, show how to the above clustering to construct a directed spanner and dynamic data structures that support proximity queries with respect to dc(.,)) +). Using these structures, in Scction obtain efficient approximation algorithms.
Directed spanner. Let A, B C P be two disjoint subsets of P, let P P + Reo be a metric, and let p [1, xc). We construct graph G = (V, E) and set of edge weights up: E + R20 such that the shortest path from any a E A to any b E B in the graph G with respect to weights wp is approximately dP da b).
5

---

Published as a conference paper at ICLR 2026
For each cluster C' C, create two vertices ac, by. Set V AAUBU fac,be |C E C). For each cluster C' E C, we add the following three sets of edges to E:
(i) Add the edge ac + by and set wp(ac + bc) (2r;) ) the index of is
(ii) For every a E C, we add the edge a + ac and set wppa + ac) =
(iii) For every E BnC, we add the edge by + b and set wbc - b) = O.
Clearly |V] = O(ne -1 log since |C| = O(ne log A). Since the degree of each point is
O(ynE- log n log A), the number of edges is O(n 0n log n log A). Define dG,p AxB + Ry0 as the shortest path distance in G with respect edge weights wpp.
Lemma 2.3. The weighted graph G with weights wy satisfies da I da.pla,)) b) - (4 + z)p (a, b) for all a, b E A x B and for any p E [1, o0).
Dynamic bichromatic closest pair. Let be point set, p 2 an integer, and let A, B C P be two disjoint subsets. Given weight function w: A U B - R, we define the weighted distance d : A x B + R as
(a, b) da(a, w(a) w(b).
Our goal is to maintain BCP (A, B) = arg (a B d (a, b) under insertions and deletions of points in A and B, where all inserted points come from P. We now describe simple data structure that supports maintaining BCP,,... (A, B) dynamically.
We construct the clustering C over the entire set P. For each cluster C E C, we maintain the points in n C' in a max-heap keyed by their weights, and the points in A n C' in a max-heap keyed by their weights. Let a (resp., b") denote the point the root of the heap for A n C' (resp., B E C). If C has index i, we define dc = (2r;)²- - w(a²) w(bC). We then maintain the set X ((a BO) C C} in global min-heap H. using Φc the key. The pair at the root of H is precisely BCP (A, B). The following observation is critical to the design of the BCP data structure.
Lemma 2.4. Let (a" be the pair stored the root of H. Then d d b²) men Ax n (a,
similar claim also appears in Lahn et al. (2025) as Lemma 2.3. We include the proof in Appendix C for completeness. Inserting or deleting point q A (resp., 9 B) quires visiting all clusters containing q and updating the corresponding max-heap. If the root of any such heap changes, we update the global heap accordingly. Since each point participates in dg = OlVrE= log n log A) clusters, the total update time is O(Vri= log" log A).
Lemma 2.5. Let P be set of n points in a metric space and let C be its two-level clustering. Let A, C P be two weighted point sets. There is a data structure that maintains the BCP under the distance function dw in O(1) time while allowing for insertion deletion of points, where inserting or deleting a point 9 E P takes O'dege (q) log n) time.
3 ALGRRIHTS FOR W,
In this section, we use the collection of clusters and data structures constructed in Section 2 to design two efficient algorithms for the optimal transport problem.
3.1 MINIMUM-CSS' FLOW BASED ALGORITHM
Let l and v discrete distributions with support sets and B; let |A| [B] n. We compute an approximate OT as follows. First assume 2 1 is a finite value. Let G = (V. E) be the directed graph constructed on A B described in Section 2.1] and let up be the corresponding weight function on E. We add source vertex and sink vertex f to the graph G. We also add an edge s > a for every a E with weight (p + a) = and an edge b + t for every E B with weight (b + t) = 0. This addition gives the graph single source and single sink to run minimum cost flow. Next, we assign capacities to each edge as follows: For each cluster C e C and corresponding edge ac by in F,, assign capacity of u(ac - bc) = 1. Additionally, for each

---

Published as a conference paper at ICLR 2026
a. b E C' we assign the capacity u(a ac) = u(be b) = 1. Finally, for each a E A and b E B, we assign the source and sink edge capacities as u(s + a) = pla) and u(b - t) = v(b).
We compute the capacitated min-cost max- flow f in this directed graph using the algorithm by Chen et al. (2022) in (IE))++()) log U time. Using the minimum cost flow f", we compute a transport plan a where (p(r) I (EeEr I (()))((e)) !ep as follows. Initially, o(a,b) b) = 0 for all a A, E B. While the total foww from s is positive, find any path π s + a + ac + by + b + t where (e) > 0 for every edge e on the path π and increment ( b) by = ff) I E x). Additionally decrement (e) by À for every edge until f" is zero everywhere. This concludes the construction of the transport plan Ø. It follows naturally from Lemma 2.3 that a is an approximate transport plan with respect to the W, distance. Since |ج|| = O 2 log n log A) edges with probability at least 1 1/n, the overall runtime of constructing a from l is O(n g log ri log A). This proves Theorem 1.1 for p E [1, oc).
n. We repeat
For 0o, we proceed as follows. We maintain the same source and sink vertices 8, I. as well as the same edge capacities as above. We then compute sequence of maximum flows instead of single minimum cost flow in G. and perform binary search on the radii of the clusters. By con- struction, there are at most O -15511m log A) different values of T;. For a fixed 1 S i S O( log A), define the graph G2 to be the graph G with all edges of cost wr (ac - be) 2r, removed. Compute maximum flow fr in G2 from to t in (n23/3g- log log A) A))+a(1) log U time using the algorithm of Chen et al. (2022). If E.EA Je(8 > a) = 1, then conclude that (,,(,,)) v) I 2r; and decrease i. Otherwise, conclude that W (1,, > 2r, and increase i.
Let it be the smallest value of i such that E.ea fi(s + x) = 1. Then we compute transport plan a from fi+ as above in the case when p < oc. Initially, a(a, b) = 0 for all a E A, E B. While the total flow from 8 to is positive, find any path क = -> a + ac - by + -> I where
fi- (e) >> O for every edge on the path π and increment aa, b) by À min [fi((( (e)
E x).
Additionally decrement (e) by À for every edge e E æ. We repeat until f:- is zero everywhere. This concludes the construction of the transport plan It follows naturally from Lemma 2.3 that is an approximate transport plan with respect to the Woo distance. Similar to the algorithm for finite p. we observe that the overall runtime of constructing a from si- is O(n 8 log n log A). This
proves Theorem
1.1 for p
3.2 COMBINATORIAL ALGORITHMS
Classical capacity-scaling algorithms for the OTT perform O(n U) Hungarian searches (i.e., dual-adjustment steps) to compute a solution. Each Hungarian search can be implemented by in- serting/deleting each point of A U B (resp., A) into BCP data structure at most once Lahn et al. (202);; [ttiinson & Vaidya 1995); Agarwal & Sharathkumar (0122). Since inserting or deleting a point q requires time proportional to O O'dege (4) log n),, a single Hungarian search can be executed in O(log n Eue (9)) time, which is O(k log" n) for any integer i >> 1. Theorem 1.2 follows from the fact that the algorithm performs On log U) such Hungarian searches. For the Wo-matching problem. the number of Hungarian searches can be reduced to ((vn) Agarwal & Sharathkumar (2012 2014), yielding the overall bound stated in Theorem 1.3 In this section, we describe the resulting combinatorial algorithm for k = 2, which runs in Ö(n²) time and computes a (4 + apprroximatee W. pyp-mhhhhing between two point sets A and B of size n.
We begin by constructing the two-level clustering C on the input point set AUB. We describe the algorithm and analyze its running time under the assumption that the degree of every points of BUB O(vn, -1 log n log A). Since this degree bound holds with probability at least 1//n, the same gigh-probability guarantee applics to the overall algorithm. The algorithm selects an appropriate scaling parameter > and simulates a single scale of the -abow-Tarjan cost-scaling algorithm for bipartite matching. We begin by defining scaled costs as a scaled version of the p-th power of the proxy distance o(a, b) = i (a b)]. The algorithm proceeds with these integer costs č(a, b).
Matchings and augmenting paths. A matching M is collection of vertex-disjoins edges. vertex not incident to any edge of M is said to be free. A matching is perfect if no vertex is free. Given matching M, alternating path is a path whose edges alternate between those in M and those
7

---

Published as a conference paper at ICLR 2026
outside M. An augmenting path is alternating path whose two endpoints are free. Augmenting along such a path flips the membership of its edges in M, thereby increasing the size of M by one.
1-feasible matching. Each vertex A. U B is assigned an integer dual variable y(u). A matching M and dual weights y(-) are 1-feasible
v(a) + u(b) < C(a, b) + v(a) + y(b) = ê(a, b)
for all (a, b) E A x B, for all (a,b) b) E M.
(2)
We define the slack of an edge (a, b) with respect to a matching M and dual weights y/-) as
0, s(a, b) ((aab) - v(a)) u(b) +1, if b) & M.
if (f, b) M,
An edge is admissible if s(a, b) = 0, and the set of admissible edges forms the admissible graph
We initialize the matching M = 0 and set all dual weights to zero, i.e., y(v) = 0 for every E AUB. Note that (M, y) is 1-feasible. Let BF B be the free vertices of B with respect to M. The algorithm maintains 1-feasible pair ( M, 1) consisting of a matching M and dual weights y(.), and executes iterations. Each iteration has the dual adjustment step, which builds an augmenting path of admissible edges, and augmentation step. which computes a maximal set of vertex- -iitx-- augmenting paths and augments the matching along these paths to increase the size of the matching. Next, we describe the dual adjustment and the augmentation steps.
Dual adjustment via BCP-based Hungarian Search. The Hungarian search procedure runs a Dijkstra-style shortest path search using slacks as edge lengths. This search is implemented using bichromatic closest pair (BCP) queries, described in Section 2.1] with implicit dual updates. The search maintains tree. Let U' C B: the set of vertices of B already added to this search tree and let V C A be the set of vertices of A not yet added to the search tree. Initially, U contains all free vertices in B, each with distance label la 0, and V A. Define effective weights
(bb = 5(y(b) (b) for b E U. w(a) vyaa) for a V.
At each iteration, select the edge (a,b) = arg mi (sl', (q(a",,') b') + fv} which show (in Appendix DD is exactly the weighted BCP maintained by our data structure, i.c. arg min,'ev, {d ,6) w(a') w(b')} Remove from V. set la = 4, + s(a b). and add it to the search tree. If is free, an augmenting path has been found and the search terminates. Otherwise, let b' be its matched partner; set by = la update w(b') y(b') ler, and insert D' into U. Once the search terminates, the dual weights are updated as
y(a) < y(a)= - la- + to for all a E S,
y(b) - v(b) + la- - lo for all b E T,
where _ A and T S B are the sets of vertices reached.
Efficiency of Hungarian Search. The BCP data structure maintains the weighted closest pair be- tween the sets U and V. Initially, we set U Br and V A. In each iteration, one point of A is removed from V, and at most one point of B \ U is inserted into U. Inserting a point into the data structure requires O(degc(u) log time. Since each vertex in A U B is inserted or deleted at most once, the total running time of Hungarian search is O (EUEAUB dege(u) logn). 0 A log" n).
Once the dual- adjustmettt phase reaches a free vertex in A, the search guar- antees the existence of at least one augmenting path in the admissible graph. The algorithm then computes a maximal set of vertex-disjoin augmenting paths via a sequence of partial depth-first searches (DFS). We initiate DFS from each free vertex of B, one at time. Let X C A denote the set of vertices in A that have not yet been visited by any DFS: initially, X = A. For every cluster C' E C, we maintain max-heap over the points in X n C. keyed by y(a) for E X. In addition, each point E AU maintains the list of clusters in which it participates.
Augmentation step.
The DFS alternates between unmatched admissible edges from B to A and matched edges from A back to B. When the search is at a vertex u E B, it scans the clusters containing u. For each such cluster C (say of index i), let be the root of the max-heap at C. It checks if (2r:))/A4 /8/ v(u)
8

---

Published as a conference paper at ICLR 2026
((0)+1 +10+11= 0, i... (u, o) is admissible. If not, we remove C. from the cluster list of u and continue to the next cluster. If an admissible edge to some vertex v is found, we remove v from X (i.e., delete v from all max-heape corresponding to clusters that contain v) and extend the alternating path. If is matched to some u', we append the edge (tz, ) to the path and continue the DFS from u'. If instead is free, an augmenting path has been discovered. We prove the correctness of the DFS procedure in Lemma D.2 After completing DFS, the algorithm proceeds with the next free vertex of B. The process terminates once DFS has been initiated from every free vertex in B. Finally, all discovered augmenting paths are flipped simultaneously to update the matching. For each augmenting path and for every vertex b E B lying on it, reduce the dual weight by one: y(b) + y(b) 1. This correction guarantees that all newly matched edges remain tight under the 1-feasible condition.
The algorithm executes dual adjustment and augmentation until all vertices are matched.
Efficiency of partial DFS. Initializing the max-heaps across all clusters requires 0(13/2=-1oog log A log" n) time. When exploring a vertex u, we sean the clusters in its cluster list. If cluster does not yield an admissible edge, it is removed from the list. Each cluster can be removed from u's list at most once, so the total number of such events for a vertex u is bounded by (e) for u and cumulative across all A, the time taken by these events can be bounded by O log A log" n). Finally, when an admissible edge is found a vertex u, u is deleted from X requiring max-heap updates which take O(dege(u) log time. Since each vertex is processed at most once in this manner, the overall cost remains at Onn3/2g-1 log A log" n). Thus, the total running time of the partial DFS phase is 0nn3/2g-1 log A log" n).
Overall Efficiency. We select the parameter d (using O(log guesses) so that the edge costs be- come integers and the sum of the costs of the edges in the optimal matching (referred to as total cost) is scaled to ((n/)). Scaling by å preserves the true optimum total cost, while rounding intro- duces at most an additive error of n in the total cost. Moreover, the 1-fcasible matching produced is itself within +n of the rounded optimum (Gabow & Tarjan (1989). Hence, the total deviation is at most 2n, and whenever the rounded optimum is at least 2n/e, the resulting solution is guar- anteed to be within a (1 + E) factor of the true optimum. Gabow and Tarjan showed that if the costs are integers and the value of the optimal solution is O(n/s), then single scale of their al- gorithm converges in O( O( /n/EE) phases. In particular, when the optimal total cost is S 2n/s, the algorithm terminates in O(vn/s) phases, and combined with the error bounds above, produces a (1 + é -approxximation. Each of the two steps-dual adjustment and augmentation--ff a phase can be implemented in On/3/2g log A log" n). Combined across ((/n/)) phases, the total execu- tion time is bound in Theorem 1.3
EMPIRICAL EVALUATION
This section contains an empirical evaluation of the clustering method from Section 2 as well as the approximation factor obtained from the primal-dus algorithm of Section(.22 We also compare the performance of our algorithm to that Lahn et al. (002).. Both methods compute an optimal matching with respect to proxy distance using the same primal-dual framework. Our approach employs a clustering- worst-case factor of (4 + =). whereas Lahn et al. (0225) uses an HST-based proxy distance (HST- dist) with an Oflog n. distortion guarantee. Computations were performed on a computer with an 8-core Apple MI CPU with 16GB RAM. Samples are drawn from uniform and truncated normal distributions on the unit cube in up to 10 dimensions, and we additionally evaluate the algorithm on samples of the MNIST Letas [1988,, treated as 784-dimensional data.
proxy distance (cluster-dis)) that approximates true distances within a
Proxy distance accuracy. We first evaluate the quality of the clustering by comparing the induced cluster distances to the ground metric. For each value of n, we measure both the maximum and the average distortion across all pairs. Figure confirms that the worst-case distortion never exceeds the theoretical (atex-approximmatio guarantee of Lemma and that the average distortion is often substantially smaller-typically close to factor of 2 indiiaatinn that the effective approximation in practice is far tighter than the worst-case bound. For comparison, we perform the same evaluation for the HST-dist of Lahn et al.] (2025) for p 2 and on the uniform distribution. As shown in Figure [el both the average and the maximum distortion incurred by HST-dist are considerably larger than those of cluster-dist, highlighting the improved accuracy of cluster-dist.
9

---

Published as a conference paper at ICLR 2026
Bumoer Doerebbms
227
Conpane
(a) Approximation Ratio (Normal) (b) Approximation Ratio (MNIST,
(c) Algorithm Efficiency
IVB
TES 210 PADNE 2Oo7
5900 3000 1500 19300 25300 30.00 3000 40306
(d) Average Degree
Max and Average Distortion
(f) Ratio of p-Wasserstein Cost,
p=
Figure 1: Empirical evaluation of the 2-layer clustering and mmtthhiin algorithm.
Algorithm accuracy. To evaluate the accuracy of the primal-dual matching algorithm, compare the computed matching cost to that obtained using the exact distance matrix. Figures (2) and (Ta) report the approximation ratio across values of 2,, 3, 4, 5, oo}. The ratios consistently remain well within the theoretical (4+)) factor, with typical values close to 1.5-2, again suggesting that the empirical performance is considerably better than the worst-case analysis. This trend is stable across both uniform and normal distributions, as well as MNIST data. We further provide the ratio of the optimal matching costs under the HST-dist and Cluster-dist for p 2 and the uniform distribution in Figure I The results demonstrate that the optimal matching cost computed using cluster-dist is consistently more than 4.5 times smaller than the cost obtained using HST-dist.
Clustering efficiency. Next, we examine the degree, ie. the number of clusters each point partici- pates Figure (Id) shows that the observed average degrees closely track the theoretical bound of Lemma 2.1 for dimension d < 10, and both distributions. This indicates that the two-layer clustering both space-efficient and stable across different settings.
Algorithm efficiency. We measure efficiency by the number of bichromatic closest pair (BCP) queries, which dominate the running time. As shown in Figure (Ic). the query counts scale as predicted and remain nearly identical across all choices of p. Combined with the Ö(n 1/) per-query complexity, this provides strong empirical evidence that the algorithm runs in quadratic time and scales smoothly with problem size.
Summary. Overall, the experiments demonstrate that the proposed method is both theoretically grounded and empirically robust. While the theoretical analysis guarantees only a (4 + 8) approx- imation, the observed approximation ratios are consistently much smaller, indicating that the algo- rithm performs substantially better practice. The clustering step is efficient both time and space, and its distortions far below the worst-case bound. These results suggest that approach is a practical alternative to additive methods such as Sinkhorn, particularly in regimes where existing techniques either fail apply (e.g., g.. p oo) require tigherrthan-quadratee time.
10

---

Published as a conference paper at ICLR 2026
5 REPRODUCIBILITY
All experimental results were produced using a proof-of-concept python implementation. The source code for the implementation can be found httpts::/aannnnnnoe/se 4open. science/ E/anon_man manzon_matchinn-28333
ACKNOWLEDGEMENT
Work by Pankaj Agarwal and Keegan Yao has been partially supported by NSF grants CCF-22- 23870 and IIS-2402823- and by US-Israel Binational Science Foundation Grant 2022/131. Work by Oliver Chubet and Sharath Raghvendra is supported by an NSF Grant CCF-2514753. We would like to thank the anonymous reviewers for their useful feedback.
REFERENCES
Pankaj K Agarwal and R Sharathkumar. Algorithms for the transportation problem in geometric settings. In Proc. 23rd Annual ACM-SIAM Sympos. Discrete Algorithms, pp. 306-317, 20,6317,,011 2012.
Pankaj K. Agarwal and R. Sharathkumar. Approximation algorithms for bipartite matching with metric and geometric costs. In Proc. 46th Anmual ACM Sympos. Theory of Comput. pp. 555-564, 2014.
Pankaj K Agarwal, Hsien-Chin Chang, Sharath Raghvendra, and Allen Xiao. Deterministic, near- linear e-approximation algorithm for geometric bipartite matching. In Proc. 54th Annual ACM Sympos. Theory of Comput., pp. 1052-1065, 2022.
Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Keegan Yao. Fast and accurate approximations of the optimal transport in semi-discrete and discrete settings. In Proc. 35th
Anmual ACM-SIAM Sympos. Discrete Algorithms, pp. 4514-4529. SIAM, 2024.
Jason Altschuler, Jonathan Weed, and Philippe Rigollet. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. In Proc. Advances in Neural Information Processing Systems 30, pp. 1964-1974, 2017.
David Mlvarez-Melis and Tommi Jaakkola. Gromov-wasserstiin alignment of word embedding spaces. arXiv preprint (ttv://800000//1 2018.
Alexandr Andoni and Hengjie Zhang. Sub-quadratic (1+ approxxmmmte euclidean spanners, with applications. In Proc. 64th Annual IEEE Sympos. Foundations of Computer Science, pp. 98-112, 2023.
David S. Atkinson and Pravin M. Vaidya. Using geometry to solve the transportation problem in the plane. Algorithmica, 13:442-461, 1995.
Franz Aurenhammer, Friedrich Hoffmann, and Boris Aronov. Minkowski-type theorems and least- squares clustering. Algorithmica, 20(1):61-76, 1998.
Baruch Awerbuch, Bonnie Berger, Lenore Cowen, and David Peleg. Near-linear time construction of sparse neighborhood covers. SIAM Journal Computing, 28(1):263-277, 1998.
Arturs Backurs, Yihe Dong, Piotr Indyk, Ilya Razenshteyn, and Tal Wagner. Scalable nearest neigh- bor search for optimal transport. In Proc. 37th Annual International Conference Machine Learning, pp. 497-506, 2020.
Surender Baswana and Sandeep Sen. A simple and linear time randomized algorithm for computing sparse spanners in weighted graphs. Random Structures & Allgriihhms, 30(4):532-563, 2007.
Soheil Behnezhad, Mohammad Roghani, and Aviad Rubinstein. Approximating maximum match- ing requires almost quadratic time. Proc. 56th Anmual ACM Sympos. Theory of Computing, pp- 444 454, 2024.
11

---

Published as a conference paper at ICLR 2026
Lorenzo Beretta, Vincent Cohen-Addad, Rajesh Jayaram, and Erik Waingarten. Approximating high-dimensional earth mover's 's distance as fast as closest pair. arXiv preprint arXiv:2508.06774, 2025.
Jean Bourgain. On lipschitz embedding of finite metric spaces in hilbert space. Israel Journal of Mathematics, 52(1):46-5,, 1985.
Yann Brenier. Polar factorization and monotone rearrangement of vector valued functions. Commu- nications on Pure and Applied Mathematics, 44(4)3375-477, 1991.
Wanxing Chang, Shi, and Jingya Wang. Csot: Curriculum and structure- optimal transport for learning with noisy labels. Proc. Advances in Neural Information Processing Systems, 36: 1528-8541, 2023.
Moses S Charikar. Similarity estimation techniques from rounding algorithms. In Proc. 34th Annual ACM Sympos. Theory of Comput., pp. 380-388, 2002.
Li Chen, Rasmus Kyng, Yang P Liu, Richard Peng, Maximilian Probst Gutenberg, and Sushant Sachdeva. Maximum flow and minimum-cost flow in almost-linear time. In Proc. 63rd Anmual IEEE Sympos. Foundation of Comput. Sci., pp. 612-623, 2022.
Ching Yao Chuang, R Devon Hjelm, Xin Wang, Vibhav Vineet, Neel Joshi, Antonio Torralba, Ste- famie Jegelka, and Yale Song. Robust contrastive learning against noisy views. In Proc. IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 16670-16681, 2022.
Edith Cohen. Fast algorithms for constructing t-spanners and paths with stretch t. SIAM Journal on Computing, 28(1):210- 236, 1998.
Scott Cohen and L Guibas. The earth mover's distance under transformation sets. In Proc. 7th IEEE International Conference Computer Vision, volume 2, pp. 1076-1083, 1999.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Proc. Advances in Neural Information Processing Systems, 26, 2013.
Sebastian Damrich, Philipp Berens, and Dmitry Kobak. Persistent homology for high-dimensional data based on spectral methods. Proc. Advances in Neural Information Processing Systems, 37: 41954-42014, 2024.
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's algorithm. In Proc. 35th Annual International Conference Machine Learning, pp. 1367-1376, 2018.
Emily Fox and Jiashuai Lu. A deterministic near-linear time approximation scheme for geometric transportation. In Proc. 64th Annual IEEE Sympos. Foundations of Computer Science (FOCS), pp. 1301-1315. IEEE, 2023.
Harold N Gabow and Robert E Tarjan. Faster scaling algorithms for network problems. SIAM J. Comput., 18: 18110111033, October 1989. ISSN 0097-5397.
Sariel Har-Peled, Manor Mendel, and Dániel Oláh. Reliable spanners for metric spaces. ACM Transactions on Algorithms, 19(1):1-27, 2023.
Arun Jambulapati, Aaron Sidford, and Kevin Tian. direct Ō(1/e) iteration parallel algorithm for optimal transport. Proc. Advances in Neural Information Processing Systems, 32, 2019.
Kwangho Kim, Jisu Kim, Manzil Zaheer, Joon Kim, Frédéric Chazal, and Larry Wasserman. Pllay: Efficient topological layer based persistent landscapes. Proc. Advances in Neural Information Processing Systems, 33:15965-15977, 2020.
Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Weinberger. From word embeddings to document distances. Proc. 32nd Anmual International Conference on Machine Learning, pp. 957-966, 2015.
12

---

Published as a conference paper at ICLR 2026
Nathaniel Lahn, Deepika Mulchandani, and Sharath Raghvendra. A graph theoretic additive ap- proximation of optimal transport. Proc. Advances in Neural Information Processing Systems, pp. 13813-13823, 2019.
Nathaniel Lahn, Sharath Raghvendra, and Kaiyi Zhang. A combinatorial algorithm for approximat- ing the optimal transport in the parallel and mpe settings. Proc. Advances in Neural Information Processing Systems, 36:21675-21681 2023.
Nathaniel Lahn, Sharath Raghvendra, Emma Saarinen, and Pouyan Shirzadian. Scalable approxi- mation algorithms for p-wasserstein distance and its variants. In Proc. 42nd Annual International Conference on Machine Learning, 2025.
Yann LeCun. The mnist database of handwritten digits. http://yanm. lecun. com/exdb/mmistl/ 1998.
Manor Mendel and Assaf Naor. Ramsey partitions and proximity data structures. Journal of the European Mathematical Society, 9(2):253-275, 2007.
Mikkel Thorup and Uri Zwick. Approximate distance oracles. Journal of the ACM (JACM), 52(1): 1-24, 2005.
Pravin M Vaidya. Geometry helps in matching. SIAM J. Comput., 18(6):1201-1222 1989.
Siddharth Vishwanath, Kenji Fukumizu, Satoshi Kuriki, and Bharath K Sriperumbudur. Robust per- sistence diagrams using reproducing kernels. Proc. Advances in Neural Information Processing Systems, 33:21900-21911, 2020.
Mikhail Yurochkin, Sebastian Claici, Edward Chien, Farzaneh Mirzazadeh, and Justin Solomon. Hierarchical optimal transport for document representation. Proc. Advances in Neural Informa- tion Processing Systems, 32, 2019.
13

---

Published as a conference paper at ICLR 2026
A CONDITIONAL HARDNESS OF THE Wo -MATCHING PROBLEM
In this section, we prove a conditional hardness result for approximating the Wo distance.
Suppose we given the bipartite graph G E), where V Vi UV2 with [il, [V2l n, Vi n V/2 Ø, and E C Vi x 1/2 We will reduce this problem of determining if G contains a perfect matching computing an approximate W -isdtittncee
Construct the following metric on V. For each distinct v1, 122 E V. (v1, v2) E E or (#2, 01) E E then define ((1 U2) L. Otherwise, define p(v1, 02) 2. For completeness, one can define p(v, v) = 0 for all E V. The finite metric space pV) can be constructed in O(n²) time, given G. is easy to see that (V, p) is a metric space: (i) by definition, p(e: v)
and p(u, v') > 0 for all r' * u; (ii) if p(11. 12) = 1 for some # U2, then either (v1. U2) or (U2, v1) is in E, in which case p(u2, 01) = 1; (iii) p(v1, v2) S and p(u), v2) + p(03, 02) 2 1+ 1 = 2 for vi, U2 E and v3 o v1, 02, implying triangle inequality in combination with observation (i) to prove the degenerate case v/3 = 01 or 02- = 022
Additionally define the distributions M41 V1 + [0, 1] and μ2: V/2 - [0, by Mr (1(11) = (2(22) = 1. for all 11 E Vi and all 12 E V/2 Assume we are given an approximation algorithm A for the Wo distance, and let a denote the transport plan from μ1 to #l computed by algorithm A.
Since p(11, v2) E [1,2) for all 1/, U2 V1 x V22 must also be the case that
Woo (a) ==
max 2.y1a (2,8))
p(z, y) E {1,2]
for any transport plan a and therefore Wxx (11, E (1, 2). We use this observation to prove the following crucial relationship between W (11, P2) and the maximum cardinality matching in G.
Lemma A.1. A perfect matching exists in G if and only if W Wo 12) 1.
Proof. If Wo (11. = min, wp((r) = 2. then it is impossible to construct a transport plan o where every edge has distance 1 and therefore at least one edge (u, u) where ou, v) > 0 satisfies p(u, v) = 2. By definition, the metric p is equal to 2 if and only if the pair is not an edge in G. Therefore, any matching must have size at most 1.
If Wo (141, (22) = 1, then it is possible to construct a transport plan a where every edge has distance and by standard network flow theory this transport plan a convex combination of matchings in (ViVV, Vi V2). By definition of p, we note p(U1, 02) = 1 if and only if (v), U2) E E. We conclude that any matching M where (11, 0 for every (U1, U2) E M is also a perfect matching in G.
Therefore, it suffices to determine if either W (1g (2) = 1 or W (M1,P2) = 2, and extract any matching from the resulting transport plan in O(n)) time in the former case.
We briefly describe this (standard) procedure of how to construct a perfect matching M in V11 x V22 using a in O(n²) time such that for every (vi, U2) M, we have (1 U2) 0. To do this we iteratively choose an arbitrary edge (u, v) E V1 x Vz such that auu, v) > 0, add (u, v) to M, remove u. from V1 and v from 1/2. Repeat until M is perfect matching. Since is known to be a convex combination of perfect matchings, we conclude that the algorithm always constructs a perfect matching.
What follows is an observation that even when given relatively weak approximation algorithms for the oo- Wasserstein distance, it is possible to distinguish between these two cases. We first prove this for the case that A is a relative approximation algorithm.
Lemma A.2. Suppose for some $ >> O there exists an algorithm running in O(n²) time which com- putes a (2- -appproximaaee transport plan under the W so metric between two discrete probability distributions with supports of size at most n. Then, for any bipartite graph G, one can compute perfect matching in G or conclude that none exists in Onn²) time.
Proof. We assume the existence of algorithm A which, when given input finite metric space (V. p) and distributions μ. v, computes a transport plan a such that w. (a) S (2 - =) W (H,2) in Onn²) time.
14

---

Published as a conference paper at ICLR 2026
First suppose that ( 2. By the approximation guarantee of the algorithm A,
conclude that W sc (12) 2 2-5 > 1 and therefore (x1 Wxo(H11.H2) P2) = 2. Now suppose that
2
Woo (a = 1. Then we have immediately found a minimizing transport plan and conclude that Wo (MI, u2)
We conclude that if A is a (2 -)approoxmmation algorithm then necessarily Woo (a.A) = Wo (2) for the constructed input metric space (V, p) and distributions μl, 42. The result follows after simple application of Lemma A.1
(22
We emphasize that the metric space used to prove Lemma A.2 has spread A 2. Therefore, Lemma A.2 is useful even the algorithm A has ((n2((A) runtime for any function f IR20 + Rz0 independent of n. Next, we prove the analogous claim for the case where A is additive approximation algorithm.
Lemma A.3. Suppose for some E > there exists an algorithm running in O( time which computes a ( adddiiivvee approximate transport plan under the We metric between two discrete probability distributions with supports of size at most n in a metric space with diameter A. Then, for any hipartite graph G, one can compute a perfect matching G conclude that none exists in (n²)) time.
Proof. We assume the existence of an algorithm A which, when given an input finite metric space (V,p) p) and distributions J4, P, computes transport plan such that wos (a) S W (14, v) ( in (n²) time.
First suppose that (oo(o = 2. By the approximation guarantee of the algorithm A. we conclude that Wo (41, 2 (x
=) = 2- 2-(1-2)) > and therefore Wo (41, = 2. Now suppose that W (o = 1. Then we have immediately found a minimizing transport plan and We conclude that if A is a -e) apprroximaaio algorithm then necessarily Weo (oA) Woo (#1, p2)) for the constructed input metric space (V) and distributions M2. The result follows after a simple application of Lemma A.1
conclude that W (F1, (11, M2)
Then a simple combination of Lemmas A.2 and
implies Theorem 1.4
B LEVEL CLUSTERING
In Section 2 we used a two-layered clustering to approximate d up to multiplicative factor 4 + using O(n 3/2) space. This approach can be generalized to k-level clustering. Extending k-level clustering has the benefit of a reduced expected degree of each point in Po. at the expense of an increase in the guaranteed stretch factor of the data structure.
In this section, we describe clustering based distance function that can be constructed in O(n² + 1+1k log A) time, approximates d up to a factor of 2k(1 + E) and that can be repre- sented using O(k O[kn++1/kg-11 log A) space. Similar to Section 2.1 we can use this clustering to construct a spanner and maintain (weighted) bichromatic closest pairs. Then once we have those data structures, the algorithms in Section 3 work in the same manner. Applying the combinatorial algorithm in Section 3.2 to the k-level clustering for i gives the (6 + appppoxxiiatoo result
in Theorem 1.3
kelleeel clustering. We now construct a ke-leeel clustering of points of P. Let Po P. For i i= I,...,, h 1, we next choose subset P. C P. by sampling each point in P: independently with probability n -1/k The expected size of P: is E E[[ = n 1-ih for each i < k 1.
Set /. = [logrt+) A1. Let ro = 0 and r; = (1 + e)² for S i S t. For each S 1 < k 1,
E Pt P41 and S i < t, define the cluster C, as
c {x E V(g, P1+1) d(x; 9) < r.).
Finally for each E Pa-1 and 1 S i < t, define the cluster C₁[1] as
c (x Pu d(x, q) $ r:).
Lemma B.I. The k-level clustering can be constructed in O(n² kn 1/kk - log A) time.
15

---

Published as a conference paper at ICLR 2026
Proof.
The samples 0 = Ps C For each from A: - 1 to 0, we do the following:
can be computed in O(k)) time.
Compute the distance of each P to the set P.+1, where the distance from any E P to Ps is defined to be co. This takes n/ Pi41l time, which is O(n 0 Osis S 2, and constant for k 1.
O(n +1/k) in expectation for
For each E P. P and for each P compute the Voronoi region (ax, P.) by adding x to (ar, P.+1)
) if (a p) < d(x, P.+1). This takes O(n)P:) time; that is, (nn++//k) for i> 0 and (n for i i= 0.
For each p E P. P.+1. each E Vp,, P.+1), and for each index 1 J 6, add C, if d(r,p) $ T. Throughout the entire construction, this takes time proportional to the sum of the degrees over p E P. P.O 1+1/k log A).
Thus, the clustering is computed O 4 log A) time.
q P,i t be the collection of all clusters. Again define the degree of point E P as dege (p) == HC E C p E CH. In an analogous manner to Lemma 2.1 we prove that the expected degree of each point in P small.
Let C [Cg
Lemma B.2. E(dege(p) O(kn g log A) for all p P.
Proof. We note that for any 0 S S t and for any E Pi P.+1 for i < k -- 1, we have C1li] S (gg, P.+1). There are at most O( log A) different values of Therefore suffices to prove for any p E P, the number of points q E P, P.+1 where E V(g, P.+1) O(kn'A) in expectation.
Let x E Po and fix i < k - 1. Let Wi....y Wny be the elements of P. in order of non-decreasing distances to Note that P.+1 C P.. Let be the smallest index such that W. P.+1-
For any index r > we have d(x, wr) > d(x, w.). Since w, P.+I. the point x is closer to W, than to Wrr, and hence
Vww,, P.+1). Thus can only belong to V(w,, Pi+1) if < 8 and therefore the number of Voronoi sets containing x at level is at most s. Since each point in P.+1 selected from P. independently with probability n 1--1 1/k the random variable s is geometric with success probability n n-1/k 1/k, Hence, the expected value of s; is E[s] l/k
Cluster-induced distance approximation. As in Section for any cluster = C, and any pair of points x, y E Cyli]. define the cluster-induced distance between æ and y as
da 3) 2r:-
Then define C(x, y) {C e 2, C} and de (x. 3) mingec(#..)) dz y). We prove that this minimum cluster-induced distance approximates d(x, y) within a factor or 21: e))
Lemma B.3. d(x,g) y) < de(x, y) < 2k(1 + +) d(x, y)
Proof. First we show that if æ and are separated by P., i.c. if there exist a, 3 E Pt such that d(x, a) < d(z, y) and d(y, (8 d(x, y), then d(y, Pt) 2 d(y, Pt+1) d(z, y). Let and 8 be the closest points to x and y in P141. respectively. Let and be the closest points to x and y in P. respectively. Then
d(z,a) a) 2 d(y, a) d(z.y)
(triangle inequality] [definition of [
2 d(y,3) d(x, y)
Let be inserted in round I, Then y & [2 if and only if d(y, Pi-1) < d(x,v) Let C2[4] be the cluster containing both x and y that minimizes d(y, z). Assume that C2 [1] was inserted in round Note j < and - < k - 2. Then by the claim above,
d(y. 2) < d(y, P(t)) (1 (2)
16

---

Published as a conference paper at ICLR 2026
s ((k - 2) + 1)d(x.))) v).
Then by the triangle inequality it follows that (xx,)) S (a, y). Note that d(x, 2) $ T3 < (1 + e) z), where C₂ has minimum index i containing both x and y. Without loss of generality assume ( 2) 2 d(y, 2). Therefore de(x,)) y) = 2r; S 2(1 + E)d(2, 2) S 2k(1 + 2)d(r,y).
C OMITTED PROOFS FROM SECTION 2
We omitted many of the proofs of technical lemmas from the main text for the sake of space. In this section, we conclude with any missing details.
We first prove that the shortest path distance dG.p in the directed graph G with weights up approximates d(
Let p oo) be an arbitrarily chosen value, and let a E A,b B be arbitrary points. We note by construction of G, for any path π from a to b in the graph G, there exists unique edge e (ac + bc) E E for some cluster C' C such that w(e) > 0.
Proof of Lemma
2.3
To prove the upper bound, we show that there exists a path where the edge = (ac + bc) satisfies up(e) S (4 + =)P da, b). By Lemma 2.2 observe that there exists cluster C E with index i where a, b E C and 2r; S (4 14++2 =)- d b). Since a E C, we observe that the edge a - ac exists
Similarly, since E C, the edge be - exists in G. Therefore, the path a - ac + be > b is path from to b in G. The weight of the edge ac - be is wp wp(ac + bc) = (2r,)²= S ((4+ + B) d(a, 6))" = (4 + E)P. (a,
To prove the lower bound, we show that for any such path from to in G, up(e) 2 (a, b). Suppose क = a - ac + by + b is a path from a to b G. By construction, observe that the path ac bc - is path in G for some cluster C E C with index i, then both a and b are contained in C. Since a and b are contained in C, there exists a point E C (in particular, choose the center of the cluster (i(ā)) such that d(a, g) & r; and d db, q) I ri. Conclude that d(a, b) S d(a, 7) + d (q. b) & 2r; by triangle inequality and therefore d' (a, b) S (2r;)² = wp(ac - bc). L.
We next prove that the global heap H constructed for the BCP data structure will contain the weighted bichromatic closest pair at its root.
x B, let C(a, b) denote its witnessing cluster, i.e., the cluster of smallest index that contains both a and b. Let i(C(a, b)) denote this smallest index. By definition of de.
Proof of Lemma 2.4 For any pair (a, b) E Ax
de(a, b) 2r,((((,))))
Hence,
- w(a) d da, b) = wwaa) - w(b).
Now fix a cluster E with index i. Let
C E arg max w(a); aEAnC
bc E arg max w(b). BEBNC
For any (a, E (AnC) x (BN C), we have
(2r,)* w((a) w(b) (2r.)² -w(a -w((()))) oc.
Thus, among all pairs whose witnessing cluster is C. the smallest weighted distance is exactly éc, achieved by (aC)))). Since every pair (a, b) has a witnessing cluster C(a,, b), the inequality above implies
(a b) 2 oc(a.))
Taking the minimum over all pairs (a, b) E A x B yields
min (a,))
(a b) 2 min ΦC-
CEC
17

---

Published as a conference paper at ICLR 2026
On the other hand, for every cluster C' C. the pair (8) belongs to A x B and satisfies dw (a b%) = Φc. Hence,
min d S min dc. (a,6)
CEC
Combining the two inequalities gives
min (a, min Φc- (a.
CEC
The global heap stores one candidate per cluster with key de and returns the smallest such value. Therefore, the pair at the root of the heap achieves the smallest weighted distance, as claimed. M
Finally, we prove that updates to the dynamic BCP data structure can done in O(y log n log A) time.
Proof of Lemma 2.5 Without loss of generality assume E A. There are two steps to insert point p into the BCPs, data structure. Similar analysis also extends deletion as well any point q E B. For each cluster C' E C that contains 9, add to the max-heap of AnC. If the root of this heap is updated, from the global min-heap, we delete the previous pair and instead insert the new (a (C) with key &c (2r.)² w(aa) w(b
Thus, the total work done (d (q) log n). By Lemma 2.1 the degree can be bounded by On 2/2 f/n//llog log - log A), leading to insertion/delltion time O( log" ne- log A)
We briefly note that the details of correctness and time complexity for the weighted nearest neighbor data structure in Section were also omitted. The analysis of the described data structure nearly identical that of the dynamic BCP data structure: nearest neighbors are stored the root of the heap, where the weight is most extreme, and insertion or deletion will modify every cluster containing the point.
D OMITTED PROOFS FROM SECTION 3
Lemma D.1. If
(a. b) arg min NEU [d w(a') w(b')).
then (a, b) also minimizes
arg V'EV.W'U (s(a', (") + }
min
Proof.
Substituting w(a') 8 y(a') and wnbb) 5 (y(b) (v), we obtain
da ()) - w(aa) w(b)== d" (e'aa',) b) - óy(a') 5y(b') + 8ly.
Adding 8, dividing by 6, and taking the ceiling does not change the minimizer. Hence (a, b) also minimizes
da 8) 5y(a') 5y(b') + 56s + +81
8
Since y(a'), y(b), and by integers, and by the definition of e(a', b'), this simplifies to o(a',0) (6) - y(a') - -(d) + fo + 1. Using b') = é(a', b') - yaa') (b') + 1, we conclude that (a,b) also minimizes s(dd,8')- - lv.
Lemma D.2. Let X C A be the set of vertices not yet visited by the DFS. Let E B be current leaf in the DFS tree. Let C be a cluster of index i. If E X is the max of the heap of C keyed by dual weights and [(2r;) / y(u) u(t) + 1 = 0. then (u, v) is admissible edge.
Furthermore, suppose there any admissible edge (u, v²) incident to u in cluster C. Let be the vertex of maximum dual weight. Then the edge (u, v) must also be admissible.
18

---

Published as a conference paper at ICLR 2026
Proof. Let X S A be the set of vertices not yet visited by the DFS. The first edge found in the DFS step must be a non- matching edge because the search begins at a free point of B. Let (u, v) be an edge found in the DFS step, where E B and because whenever an edge (b, a) E B x is added to the the partial DFS tree, a is matched to a point V' also add (a,b)) and remove from X.
X. The edge (u, v) must be non mon-matching edge
So the slack of (u, 2)) is
s(u, 1) = a(u, v) - v(u) - ((v) +1= (u, u)/67 v(2) 3(0) +
Let C' be the cluster of index i from which (u, v) was obtained. The algorithm additionally guaran- tees that
[(2r:)"/0] s(u) =(v) + 1
To check admissibility of (u, u) we need to show that [d] (u, v)/57 y(u) y(v) +1 0. So sufficient to show that [de (u, w)/67 [(2r)))/6]. By the definition of de we have that
du v) = min (2r
So it must be true that (u, u) [min eeCyl (2r,) /07 s ((2r))P/0). Then by I-feasibility we have
((u) + s(u) - S a(u, - (u, v)/07 S [(2r))//)] y(u) v(v)
So we have [d u) = [(2r:)P/0T, implying that (u, v) is admissible.
Now let (u, 0²) be any admissible edge incident to u in cluster C' of index i. Let be the vertex of maximum dual weight. We show that the edge (u, v) must also be admissible.
By the admissibility of (u, o') we have [(2r:)P/0 we have (2r,))/8' y(u) + 12 y(v). Since y(v) is the maximum dual weight, this implies that y(0') = v(e). So (u, v) is also admissible.
y(u) + = y(e'). Then by I-feasibility
E PROOF OF THEOREM
We omitted formal cost analysis for the transport plan computed by the minimum-coss flow based algorithm from the main text since intuitively the fact that shortest paths in G approximate d" should imply minimum cost flow will also approximate W,, In this section, we give formal cost analysis of the algorithm in Section 3.1 for completeness. We first prove the case when p is finite.
Proof of Theorem I.7 for < 00. By standard analysis of network flow algorithms, since the capac- ity of each edge s - x is p(z) and the capacity of each edge + t is v(V), we conclude that the resulting a is transport plan. Moreover, given /", the transport plan a can be constructed in O(1 El) time since every - / path in G has constant length and at every iteration the flow along at least one edge is decremented to zero. We conclude with a cost analysis of a.
First, we prove the lower bound. By construction of the graph G, we note that any path from to must be of the form s + a ~ b + t where a ~i 6 is a path in from a E A to B. By Lemma 2.3 conclude that Eeca-ab (p(e) > (a, b) 2 dp (a, 6) for any E A, E B and any path a No in G. Moreover, every edge of the form s + a and - t has cost zero. This implies
1/p
1/P
W,(n, S uplo)=
Σ
ala, b) d' (a.
Σ (re()))p(e)) KEE
dB
We prove the upper bound. Suppose we are given the optimal transport plan a² We construct a flow f with cost at most (4 +++ E) wp (o*). For each v, y E A x B, by Lemma 2.3 there exists path x y in G such that Eeermy (e((e) s (4 5) d(z, v). Let ((z y) denote this specific path for each pair y.
Initially, set if to be zero everywhere. Then for each pair x, where o a' y) > 0, we increment flow f on every edge along s x. y + and every edge of 77 ( y) by a² (x, y). The resulting flow
19

---

Published as a conference paper at ICLR 2026
is max flow since the capacity every edge x leaving is p(z) and the capacity every edge y + entering is v(y). Moreover by construction f and the fact that every edge s - and y + t has cost zero, the cost of the flow f is
1/v
Σ f(e) wp(e) E
Σ ".yRO" (m.y)>0
(z, do.,(x,)))
1/8
Σ
0²(zz v) ((4 +)p dz(2,
u.y:o (1
I (4+ up(a") (4 +) W)(,vv).
We note that the minimum cost flow must have
than f.
Then the cost analysis for when oc follows in similar manner.
Proof of Theorem I.7 for p = co. By standard analysis of network flow algorithms, since the ca- pacity of each edge + u(²) and the capacity of each edge y t is v(y), we conclude that a transport plan. Moreover, given fi+, the transport plan a can be constructed in E(I time since every least one edge is decremented to zero. We conclude with a cost analysis of a.
the resulting
path in G has constant length and every iteration the flow along
First, we prove the lower bound. By construction of the graph G, we note that any path from to must be of the form + a ~a b + t where a -. 6 is path in G from E A to 6 E B. By Lemma 2.3, conclude that Σ Eesambb (e) > dG.1 (a, 2 d(a, for any a A, B and any path ~ D in G. Moreover, every edge of the form a + ac, bc + b, + and b - has cost zero. This implies
max WookL (oo (wx())= a.bEAxB (aab))>
d(a, S max wi (e).
of fo (c)>0
We now prove the upper bound. Suppose construct a flow f with cost at most (4 (4+32 E) woo (o*). For each x, y E A B, by Lemma 2.3 there exists path
given the optimal transport plan σ², We
in G such that Eeerny wi (e) (4 + e) d(x,y). Let (<x, v) denote this specific path for each pair x, y.
Initially, set f to zero everywhere. Then for each pair x, y where o" y) > 0, we increment flow f on every edge along s + y - t and every edge of (I. y) by o d'(x,y). The resulting flow a max flow since the capacity on every edge s + x leaving s is p(z) and the capacity on every edge y + t entering is vly). Moreover by construction of and the fact that every edge - x, y + f, x - ac and be + y has cost zero, the cost of the flow f is
max wi(e)= eE E f()
max =yEAA Mo (x.))>>
da (2r
max a.yEAxB o (a.))>0 E) d(x.
(4+ (o) (4 + Wor (1,
We note that the cost of is at most (1 +) times that of fi- since i" is the smallest index where
fa- is a max flow, i.e. has total flow 1, by definition. Rescaling & by some constant gives the desired result.
F OMITTED FIGURE FROM SECTION 4
20

---

Published as a conference paper at ICLR 2026
Figure 2: Approximation Ratio (Uniform)
21


Reviews:
[{"review_id": "vp2w31hB9t", "reviewer": "Reviewer_coG2", "summary": "The paper presents the first scalable, constant-factor approximation algorithm for the $W_{p}$ optimal transport (OT) problem that runs in nearly quadratic time. This method significantly improves upon previous $O(\\log n)$-approximation algorithms and is the first efficient, quadratic-time method to extend to the $W_{\\infty}$ distance.\n\n** Problem formulation.** \n\nThe paper addresses the $W_{p}$ optimal transport problem for discrete probability distributions $\\mu$ and $\\nu$, supported on finite point sets $A$ and $B$ within a metric space $(X,d)$. The objective is to find a transport plan $\\sigma: A \\times B \\rightarrow R_{\\geq 0}$ that minimizes the $W_{p}$ cost, defined as $w_{p}(\\sigma):=(\\sum_{a\\in A,b\\in B}\\sigma(a,b)\\times d(a,b)^{p})^{1/p}$. This formulation also covers the $W_{\\infty}$ cost, which is the limit as $p \\rightarrow \\infty$ and represents the maximum distance $d(a,b)$ for which $\\sigma(a,b) > 0$.\n\n** Main results ** \n\nThe main result (Theorem 1.1) is a randomized algorithm that computes a $(4+\\epsilon)$-approximate $W_{p}$ optimal transport plan in $O(n^{2}+(n^{3/2}\\epsilon^{-1}\\log^{2}\\Delta~\\log U)^{1+o(1)})$ expected time for any $p \\in [1, \\infty]$. Additionally, a simpler $\\tilde{O}(n^{2})$ combinatorial algorithm is provided for the $W_p$ matching problem (when distributions are uniform).\n\n** Technique/algorithm **\n\nThe core technique is a two-layered clustering scheme inspired by Bourgain's multi-level sampling, which approximates the cost function $d(\\cdot,\\cdot)^p$. The primary algorithm constructs a directed spanner graph based on these clusters and computes an approximate solution by running a minimum-cost max-flow algorithm on this graph. A second, simpler algorithm for the matching problem uses the same clustering to build efficient dynamic data structures (for bichromatic closest pair and weighted nearest neighbor queries) which are then used within a Gabow-Tarjan cost-scaling matching framework.\n\n** Experiment sumamry **\n\nThe paper provides an empirical evaluation of the simpler combinatorial matching algorithm (from Section 3.2) on synthetic data. These experiments demonstrate that the algorithm's practical approximation ratio is consistently much better than the theoretical $(4+\\epsilon)$ worst-case bound (typically around 1.5-2.0) and that the runtime scales quadratically, as predicted by the analysis.", "questions": ".", "limitations": "", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The theoretical results are solid and also discusses its potential practical application", "weaknesses": ".", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gs8SO0UJBx", "reviewer": "Reviewer_6YEK", "summary": "The paper presents two constant factor approximation algorithms for optimal transport. The first one, relies on computing a min-cost flow, but does so on a smaller graph that approximates distances of the full graph. The second algorithm is for the restricted setting where the distributions are uniform.", "questions": "In line with the above, could please give a more detailed justification on the claimed running time?\n\nDo you have experiments comparing your matching algorithms to others both in terms of efficiency and accuracy?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "Optimal transport is a heavily studied and new algorithmic insights should always be welcome.", "weaknesses": "The paper emphasizes the \"scalable\" nature of the algorithm, but the general algorithm still relies on Chen et al. for computing a min-cost flow. Moreover, I believe some more detail is needed to justify the expected running time. A graph is constructed, which has m edges in *expectation* and then the algorithm by Chen et al. is applied to that graph. However, as far as I can see, this does not immediately result in an algorithm that runs in time $m^{1+o(1)}$. To conclude this, some concentration bound on the number of edges or some suitable strict upper bound on the number of edges would have to be used. In any case, this would require a little more justification.\n\nFor dense graphs, there are also other nearly linear time algorithms available (van den Brand et al. STOC 2021). Maybe the current submission can claim to save polylogarithmic factors for maximally dense graphs? But the paper presents itself as aiming to be more practical and a result like that would not necessarily qualify for a practical improvement.\n\nNot surprisingly, no experiments are included for the general optimal transport algorithm.\n\nI am not claiming that the theoretical contribution is without merit (after adding an appropriate justification for the running time bound), but that there is a gap between the claimed practical focus and this result.\n\nThe algorithm for uniform distributions does not rely on Chen et al. and so the critique above does not apply.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "KQgqPRdQmv", "reviewer": "Reviewer_okv2", "summary": "This paper studies the problem of computing the $p$-Wasserstein distance in general metric spaces for any $p\\in[1,+\\infty]$. The problem can be formulated as a minimum-cost flow problem and therefore can be solved exactly in $n^{2+o(1)}$ time using the sophisticated and impractical theoretical algorithm of Chen et al. (FOCS 2022). This paper considers approximation algorithms for the Wasserstein distance and obtains near-quadratic $\\tilde{O}(n^2)$-time algorithms with $4+\\epsilon$ approximation ratio, where $n$ is the support size of the two input distributions. This improves upon the previous $O(\\log n)$-approximation in $\\tilde{O}(n^2)$ time (Lahn et al., ICML 2025). The authors also present a simpler and more practical algorithm for the special case where the input distributions are uniform (i.e., the $W_p$ matching problem), as well as a nearly matching conditional lower bound.\n\nTechnically, both algorithms build on a clustering scheme obtained via a multi-level sampling procedure inspired by Bourgain (1985). The authors use this scheme to\n\n- construct a directed spanner with $m = \\tilde{O}(n^{3/2})$ edges that preserves the $p$-power metric $d(\\cdot,\\cdot)^p$, and by combining it with the $m^{1+o(1)}$-time algorithm for graphs (Chen et al., FOCS 2022), obtain a near-quadratic $\\tilde{O}(n^2)$-time algorithm for $p$-Wasserstein distance; and\n\n- design data structures for weighted nearest-neighbor search and dynamic bichromatic closest-pair queries, which are then used to speed up the Gabow–Tarjan bipartite matching algorithm, leading to a near-quadratic $\\tilde{O}(n^2)$-time algorithm for $W_p$ matching.", "questions": "In the “Algorithm efficiency” paragraph of the experiment section, the authors report the number of operations but not the actual running time. Could the authors clarify the motivation for this choice? \n\nMoreover, in Line 465 the paper states that *“Combined with the $O(n^2)$ per-query complexity … the algorithm runs in quadratic time ...”*. However, the results seem to indicate that there are at least $n^{3/2}$ queries, which would imply a total running time of $O(n^{3.5})$ rather than $O(n^2)$. Is this a typo, or am I misunderstanding something?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- The paper addresses the efficiency challenge of the fundamental $p$-Wasserstein distance problem and achieves significant improvements (reducing the approximation ratio from $O(\\log n)$ to  $O(1)$). Moreover, their algorithms are general and work for any $p \\in [1,+\\infty]$, whereas the previous $O(\\log n)$-approximation only works for small $p$. \n- The algorithms are also practical, as demonstrated by the experimental results.\n- The results are also complete, as the authors additionally provide a (conditional) lower bound.\n- The paper is well-written and I can easily follow the presentation and understand the main idea of the algorithm and its analysis.", "weaknesses": "The only weakness in my view is that the technical contribution is not clearly explained. The paper’s main technique, a clustering scheme via multi-level sampling, seems largely based on prior work (e.g., Bourgain, 1985). The directed spanner construction and the data structures appear to follow naturally from this clustering scheme. It would be helpful if the authors could more explicitly highlight their technical novelty (for example, by clarifying how their multi-level sampling differs from prior work). \n\nIn addition, I believe the experimental section could be further improved. In particular, it would be useful to compare against the previous $O(\\log n)$-approximation algorithm of Lahn et al. (ICML 2025), and to evaluate the algorithm on some real-world datasets.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "hkkwZSppAS", "reviewer": "Reviewer_5TsD", "summary": "The paper deals with the $W_p$ optimal transport problem. In this problem we are given two node sets $A$ and $B$, a distance metric $d$ on $A \\cup B$ and two probability distributions $\\mu$ and $\\nu$ defined on $A$ and $B$ respectively. The goal is now to move the probability mass of distribution $\\mu$ from $A$ to $B$ such that $\\nu$ is formed. Alternatively $\\mu$ and $\\nu$ could also be seen as node weights on $A$ and $B$ (each summing up to $1$) and we need to match the respective weights to each other. The respective movement or matching should then minimize the $W_p$ distance with respect to the distance metric $d$.\n\nThe authors point our that it is already known that the problem can be solved exactly with running time $O\\left(n^{2+ o(1)}\\right)$ using a min-cost flow algorithm. However, this algorithm is not feasible in practice and there has been significant research over the last years both trying to find practical algorithms as well as trying to obtain a theoretical running time in $O(n^2)$, even if this means that not the optimum solution is found.\n\nThe main contribution of the submission is a randomized $2$-layer clustering scheme that allows to create a weighted directed graph of expected size $O(n^{3/2})$ such that the shortest path distance between two nodes $a \\in A$ and $b \\in B$ in this graph is at least $d(a,b)$ and at most $(4 + \\epsilon) d(a,b)$. Crucially all paths from $a$ to $b$ only contain a single weighted edge. As a result also the value $d(a,b)^p$ gets only distorted by a value of $(4+\\epsilon)^p$ by this construction (for $p \\in \\mathbb{N}$). \n\nAfterwards the authors apply the min cost flow algorithm by Chen et al. using these new distances. Given that the cost of the algorithm depends on the number of edges in the respective graph which got reduced from $O(n^{2})$ to $O(n^{3/2})$, this improves the running time accordingly. The distance graph can be calculated in $O(n^2)$, and the author end up with an $(4+ \\epsilon)$-approximation with expected running time $O(n^2)$ under reasonable assumptions. For the case that $p = \\infty$ they provide an alternative approach with a similar running time.\n\nBesides this, the submission provides an algorithm for the $W_p$ matching problem, which can be seen as the special case of the $W_p$ optimal transport problem where both $\\mu$ and $\\nu$ are uniform distributions. The advantage of this algorithm seems to be that it has a reasonable practical running time even though the theoretical running time is mostly comparable to the more general algorithm. The approximation ratio stays $(4 +\\epsilon)$ (or $(8 + \\epsilon)$ respectively). It is also shown that if one finds an $(2 - \\epsilon)$-approximation algorithm for the $W_\\infty$ matching problem in $O(n^2)$, this would imply that one could find a perfect matching in bipartite graphs in $O(n^2)$. This can be seen as a hardness result.\n\nFor their experiments the authors generated data sets using either a uniform or a truncated normal distribution in up to 10 dimensions. They do not evaluate the practical performance of their $W_p$ optimal-transport algorithm directly (probably due to the inefficient min-cost flow algorithm). Instead they only evaluated the quality of the $2$-layer clustering scheme. They showed that the size of the clustering fits the expected theoretical bounds and obtained that in their experiments the maximum distance distortion of this clustering often was between $3$ and $3.5$ while the average distortion was closer to $1.5$. They also provided some experiments for their $W_p$ matching algorithm with the approximation ratios seemingly approaching values between $2$ and $2.5$ with larger number of nodes (at most 8000). An exception seems to be the algorithm for $p = \\infty$. Here the ratio reaches the value of $3$ and it is unclear if the value would increase even more for larger point sets. The algorithm also seems to be reasonably efficient.", "questions": "- The data structures for the proximity queries are rather involved and are not directly necessary for the optimal-transport algorithm. Maybe one could improve the readability of the paper by directly providing the algorithm after the introduction of the layer clustering and present the query structures afterwards. This would also be beneficial because the need to answer these queries only gets clear once the matching algorithm gets presented. Thus presenting this algorithm directly afterwards could be helpful.\n\n- Given that the paper focuses a lot on improving the running time it could make sense to provide a summary of the analysis of the running time of the $W_p$ optimal transport algorithm somewhere in the paper or the Appendix (which could also discusses the running time for computing the clustering). \n\n- In the experimental section there is no comparison with existing algorithms which would have been very helpful to judge whether the new results also yield improvements in practice.\n\n- In Theorem 1.3 the authors say that one could calculate the perfect matching for an arbitrary graphs (if it exists) but in the appendix they only prove that one could find the perfect matching in a bipartite graph.\n\n- In line 471 the authors write that the $W_p$ matching algorithm performs 'near optimal' in practice. This seems to be a bit of a stretch for approximation ratios around $2$.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper contains non-trivial theoretical results. At first glance the improvement from an exact $O(n^{2+ o(1)})$ algorithm to an $O(n^2)$-approximation algorithm (under certain assumptions) seems to be not that impressive. However, the authors provide a lot of literature that deal with this problem and their algorithm improves the best existing result significantly.", "weaknesses": "The fact that the new $W_p$ optimal-transport algorithm was not tested in the experiments seems to indicate that it is very inefficient in practice. It is unfortunate that the authors did not provide experimental results for existing algorithms for a better comparison.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "b81DPUwyZp", "reviewer": "Reviewer_HhqV", "summary": "This paper studies the problem of computing the $W_p$-distance, for $p \\in [1, \\infty]$, between two distributions $\\mu$ and $\\nu$ supported on finite point sets of a metric space. The authors give a $(4 + \\epsilon)$-approximate algorithm for computing the optimal transport between $\\mu$ and $\\nu$ with expected runtime $O(n^2 + (n^{3/2} \\epsilon^{-1} \\log^2 \\Delta \\log U)^{1 + o(1)})$, where $\\Delta$ is the ratio between the maximum and minimum pairwise distances of the support, and $U$ is the ratio between the maximum and minimum probabilities. Also, when $\\mu$ and $\\nu$ are uniform distributions supported on two sets of the same size, the authors give a $(4 + \\epsilon)$-approximate combinatorial algorithm with expected runtime $O(n^2 \\epsilon^{-2} \\log^2 \\Delta)$ and a $(8 + \\epsilon)$-approximate combinatorial algorithm with expected runtime $O(n^2 + n^{5/3} \\epsilon^{-2} \\log^2 \\Delta)$. Finally, the authors show that given an $O(n^2)$-time algorithm for $p = \\infty$ that is either $(2 - \\epsilon)$-multiplicative approximate or $(\\Delta / 2 - \\epsilon)$-approximate with constant $\\epsilon$, then there exists an $O(n^2)$-time algorithm for computing a perfect matching in general graphs provided exists.", "questions": "- Line 62;64: Why is deciding the existence of a perfect matching in a dense graph is a simpler task?\n- Is the approximation ratio $4 + \\epsilon$ tight for the given algrorithm?\n- Line 166: Would the degree of any particular point be as large as $nt$ instead of $n$?\n- Line 178-180: Should the first term in the equation be $\\Pr[p \\in C_{w_j}[t] \\mid w_j \\in P_0 \\setminus P_1]$? Also, should $j$ in the last term be $j-1$?\n- Line 183-185: Should the term in the summation be $\\Pr[w_s \\in P_0 \\setminus P_1] * Pr[p \\in C_{w_s}[t] \\mid w_s \\in P_0 \\setminus P_1]$? I believe the range of the second summation should be $s = 0, 1, \\ldots, n - 1$. Finally, if the previous question is correct, should the last term be $O(1)$ now that $\\Pr[w_s \\in P_0 \\setminus P_1] = 1 / \\sqrt{n}$?\n- Line 228-235: This paragraph seems confusing to me. Could you further justify its correctness and runtime? In my understanding here we should do the following: A heap for each $C$ containing $A \\cap C$ is maintained, and for each $p$, a heap containing $X =  \\{(a_C, C) \\mid p \\in C\\}$ is maintained. So each update $p$ involves updating the heaps associated with $C$ that contains $p$, which might result in updates in the heaps associated with points in $C$. Please correct me if I'm wrong.\n- At the end of Section 3.1, Lemma 2.3 is applied to conclude Theorem 1.1 for $p = \\infty$. Yet it seems to me that Lemma 2.3 only works for $p < \\infty$. Do I miss anything?", "limitations": "", "rating": 2, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "This paper studies an interesting and important problem. It is generally well-written, and the proofs are easy to follow. The technical contributions, despite mainly motivated by prior work and some classic techniques, are non-trivial and lead to improvement over prior work in terms of runtime.", "weaknesses": "Overall, I find the contributions of this paper marginal. It seems that the main technical contribution comes from the data structure given in Section 2, which at a high level is similar to Bourgain's multi-level sampling. Yet, the analysis of its guarantees follows quite straightforwardly from its construction. The remaining components of the algorithms then apply ideas from some classic algorithms.\n\nBesides, the presentation of this paper requires considerable improvement. A related work section that places this paper into a broader context is missing. The proofs could be restructured following the theorem-lemma style instead of a plain description, and the guarantees of the data structure can be stated more explicitly.\n\nI have some concerns regarding the correctness of the proofs (see Questions).\n\nDetailed comments:\n- Line 45-46: I believe that $W_p$ OT being a matching problem also requires $\\mu$ and $\\nu$ to have the supports of the same size.\n- $\\Delta$ is used to denote both the aspect ratio and the diameter, which might lead to confusion.\n- Line 200: Should $C_y$ be $C_y[i]$? Same for $C_x$.\n- Using Algorithm environments to formally describe algorithms presented in the paper would greatly improve the readability.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)