Scalable Approximation Algorithms for $p$-Wasserstein Distance and Its Variants
Authors: Nathaniel Lahn
Sharath Raghvendra
Emma Saarinen
Pouyan Shirzadian
Venue: ICML 2025 poster (2025)
Keywords: p-Wasserstein Distance
Optimal Transport
Approximation Algorithms
Primary area: optimization->discrete_and_combinatorial_optimization

Abstract: The $p$-Wasserstein distance measures the cost of optimally transporting one distribution to another, where the cost of moving a unit mass from $a$ to $b$ is the $p^{th}$ power of the ground distance $\mathrm{d}(a,b)$ between them. 
Despite its strong theoretical properties, its use in practice  -- especially for $p \ge 2$ -- is limited due to two key challenges: sensitivity to noise and a lack of scalable algorithms.
We identify noise sensitivity as a key reason why some existing approximation algorithms for $p=1$ fail to generalize to $p \ge 2$ and then present new algorithms for approximating the $p$-Wasserstein distance and its variant. First, when $\mathrm{d}(\cdot,\cdot)$ is a metric, for any constant $p \ge 2$, we present a novel relative $O(\log n)$-approximation algorithm to compute the $p$-Wasserstein distance between any two discrete distributions of size $n$. The algorithm runs in $O(n^2 \log U\log \Delta\log n)$ time, where $\log U$ is the bit-length of the input probabilities and $\Delta$ is the ratio of the largest to the smallest pairwise distance. 
We use $p$ hierarchically well-separated trees to define a distance that approximates the $p$-Wasserstein cost within a factor of $O(\log n)$ and then present a simple primal-dual algorithm to compute the $p$-Wasserstein cost with respect to this distance. 
Second, due to the noise sensitivity of the $p$-Wasserstein distance, we show that existing combinatorial approaches require $\Omega(n^2/\delta^p)$ time to approximate the $p$-Wasserstein distance within an additive error of $\delta$. In contrast, we show that, for any arbitrary distance $\mathrm{d}(\cdot,\cdot)$, a recent noise-resistant variant of the $p$-Wasserstein distance, called the $p$-RPW distance, can be approximated in $O(n^2/\delta^3)$ time.

Scalable Approximation Algorithms for p-Wasserstein Distance and Its Variants
Nathaniel Lahn't Sharath Raghvendra*
2 Emma Saarinen 2
Pouyan Shirzadian
Abstract
1. Introduction
Let μ and v be discrete probability distributions supported on sets and B, respectively, with |.|| + |B] = n. For each pair (a, b) A x B, let d(a, b) denote the distance between and b. For simplicity, assume that the diameter of AA B is i.e. (a d(u, b) = 1. A transport plan is function σ AA x B + Rz0 that assigns a mass to each pair (a, b) such that ELE aaa b) I ju(a) and ELEA (a, b) S (b)). Given a parameter p 2 1, suppose the cost of moving a unit of mass from point E to a point b B is given by d(a, b)p. The p W-sssersteinn cost of any transport plan between μ and is defined
The p- p-Wasserstein distance measures the cost of optimally transporting one distribution to another, where the cost of moving a unit mass from a to is the tth power the ground distance d( 6) between them. Despite its strong theoretical prop- erties, its use in practice especially for p 2 2 is limited due to two key challenges: sensitivity noise and lack of scalable algorithms. We iden- tify noise sensitivity as a key reason why some existing approximation algorithms for p 1 fail to generalize to p 2 2 and then present new algo- rithms for approximating the y- Wasserstein dis- tance and its variant. First, when dl:: a metric, for any constant p 2 2, present novel relative Olog n)- approximation algorithm to compute the p- Wassersteii distance between any two dis- crete distributions of size n. The algorithm runs in (n² log U log A logn) time, where log U is the bit-length of the input probabilities and A is the ratio of the largest to the smallest pairwise distance. We use p hierarchically well-separated trees to define distance that approximates the p- Wasserstein cost within a factor of Ollog n) and then present a simple primal-dual algorithm to compute the p- g-Waaserstein cost with respect to this distance. Second, due to the noise sensitivity of the p- Wasserstein distance, we show that exist- ing combinatorial approaches require Q(n² /8P) time to approximate the p- Wasserstein distance within an additive error of 6. In contrast, we show that, for any arbitrary distance d(.,,), +), recent noise-resistant variant of the p-Wasserstein dis- tance, called the p-RPW distance, can be approxi- mated in (n²² /88) time.
1/p
Up(a):= =
Σ D
o(a,b) x d(a, b)"
In the a-partial p- Wasserstein problem, we wish to compute the transport plan of that transports ar fraction of the mass and has the smallest p- p-Wasserstein cost. We refer to this cost as the a-partial p- p-Wasserstein distance and denote it by Wp,o (μ, v) = wp(ot). The p- W-sssessttiin distance, denoted by Wp(u, 2)), is equal to the cost up(ot).
A standard approach estimating the p- Wasseesstein dis- tance between μ and is to draw samples from both and construct empirical distributions pnn and Vn, by assigning mass of 1/n to each sampled point. Computing the p- Wasserstein distance between these empirical distributions corresponds to a well-known problem in combinatorial op- timization: the assignment problem. In this setting, the optimal transport plan is matching set of vertex- disjoint edges, each carrying mass of 1/n.
The p- p-aassersteinn distance is effective in capturing geo- metric similarity between distributions, especially when 2 2. When the ground distance d(.,.) is a metric, the p- Wasserstenin distance is also metric. Furthermore, when the diameter of the supports bounded, the empirical p- Wasserstein distance converges to the true distance as the sample size increases (Fournier and Guillin, 2015). These favorable theoretical properties have led to its widespread adoption in applications across machine learning (Chang et al., 2023; Chuang et al., 2022), computer vision (Backurs et al., 2020; Lai et al., 2022), and natural language process- ing (Alvarez-Melis and Jaakkola, 2018; Yurochkin al.,
"Authors are listed Alphabetically 'Radford University "North Carolina State University 'Virginia Tech. Correspon- dence to: N. Lahn <nlahn @raddoril.dd.. S. Raghvendra <skraghve @ ncsu.edu>, E. Saarinen <cesaarin @ nesuu,dduz, Shirzadian <pshirzadian@ vt.du>..
Proceedings of the nnnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Scalable Approximation Algorithms for y- Wasserstein Distance and Its Variants
2019). However, despite its strengths, the p- p-Wasserstein distance (for p 2 2) faces two significant limitations in prac- tice: sensitivity to noise and lack of scalable algorithms.
tree metric cannot approximate the 2-Wasserstein distance within bounded factor.
Subsequent work by Agarwal and Sharathkumar (2014) pre- sented an algorithm to find an ((1/89.63)) j-aprroximation of the 1- 1-Wasserstein distance in (n2++) time, which was later improved by Sherman (2017) to a (1 + s) approxima- tion algorithm in Ó /=2) time. These algorithms rely on property
Noise sensitivity: A major drawback of the p- Wassersteii distance is its sensitivity to small amounts of noise. Specif- ically, a noise mass of 8 can increase the p- p-Wasserstein distance by as much as 8 For example, a 1% noise (8 0.01) can increase the 2-Wasserstein distance by more than 10% (or 0.1) (Raghvendra et al., 2024). This makes the y- W--sssssstein distance less reliable for real-world datasets which often contain noise or outliers. To address this limita- tion, Raghvendra et al. (2024) introduced the (p, k--Robust Partial Wasserstein ((p. k)-RPW) distance. It is defined as the smallest such that Wp, Wp.(1=e)(11 2) S kg. This formulation generalizes several classical metrics including total variation, the p- p-Wasserstein, and the well-known Lèvy- Prokhorov distance (Prokhorov, 1956). Importantly, the (p, k)-RPW distance is provably more robust: modifying a 5 mass in the distribution alters the (p. k--RPW distance by at most 5, in contrast to the potentially much larger distortion in the standard p- Wasserstein case.
the assumption that the cost function is a metric that holds when p 1, but fails for p 2 since the costs are d(.,) ),, which do not satisfy the triangle inequality. As result, these methods do not extend to higher values of p.
For p = 1, significantly faster algorithms exist in geometric settings (e.g., low-dimensional Euclidean space); see, for instance, (Agarwal al., 2022; Indyk, 2007: Sharathku- mar and Agarwal, 2012a;b; Agarwal and Varadarajan, 2004; Fox and Lu, 2020). Owing to their simplicity, some of these algorithms have since been adapted to the design of data structures and algorithms for more complex tasks in the Wasserstein space; see, for instance, nearest neighbor searching (Backurs et al., 2020; Andoni et al., 2018) and barycenter computation (Agarwal et al., 2025). In contrast, for p 2 2, there are no fast approximation algorithms known that work for any metric space. The only existing fast ap- proximation algorithms are restricted to specific settings, for instance wwo-dimensionll Euclidean settings (Lahn and Raghvendra, 2019; 2021). In this paper, we address the following major open challenge:
Exact algorithms for computing the p- Wasserstein distance take Ō(n²)) time (Orlin, 1988)'. For the case where the precision of costs and probabilities are bounded, one can use interior point methods to improve the execution time to ((n2+o()) (Chen et al., 2022); how- ever, these algorithms are complex and, without significant simplifications, is unlikely that these methods will lead to a usable implementation for machine learning applica- tions. Similarly, for distances that admit dynamic weighted bichromatic closest pair (BCP) data structure with query and update time of (n), one can improve the execution time of exact algorithms to (nn2 (n)). Unfortunately, the only known BCP data structures are in the two-dimensional setting and the data structures are based on sophisticated techniques that do not have an implementation (Eppstein, 1995; Chan, 2020). Due to lack of efficient, usable exact algorithms, researchers have focused on designing relative and additive approximation algorithms.
Exact Algorithms:
For any integer p 2 and any metric d(-..), can we design an Ó(n²) time Ollog n))-elaatvve approximation algorithm to compute the p- Wasserstein distance between two discrete distributions?
Additive Approximations: A transport plan o is a 8- additive approximation, or simply 8-close, if wp(o) I w,(a" + 6. There are several simple and highly paral- lelizable quadratic algorithms for computing additive ap- proximations for any cost matrix (Cuturi, 2013; Altschuler et al., 2017: 2019; Lin et al., 2019; Dvurechensky et al., 2018; Jambulapati et al., 2019; Lahn et al., 2019; 2023). However, the best-known algorithms run in Ó(n² /8P) time, which we show is likely difficult to improve upon. Note that the execution times of these algorithms increase with p and do not converge for p = Do.
Relative Approximations: An a-relative approximation algorithm (or simply an a-approximation) returns a trans- port plan whose p- -asssersteinn cost is at most a times the optimal value. When the ground distance d(---) is a met- ric and p- 1, Charikar (2002) presented an O(n?) time Olog n log log n)) aaproximmtionn algorithm. The key idea is to approximate the original metric with a tree metric that preserves the J-Wasserstein distance within a multiplica- tive factor of (log n log l). Once embedded into the tree, the 1-Wasserstein distance under the tree metric can be computed efficiently using simple greedy algorithm. How- ever, this approach does not generalize to p 2 since the
1.1. Our Results
In this paper, we identify noise sensitivity as a key barrier preventing the extension of scalable relative and additive approximation algorithm for the 1-Wasserstein distance to the case where p 2 2, and we present new algorithms that overcome this challenge.
Relative Approximation: First, for any finite integer 2 2, we present an Ö(n²) time Ollogn n) rellavvee approxima-
O()) hides poly log n factors in the running time
2

---

Scalable Approximation Algorithms for y- p-Wasserstcin Distance and Its Variants
tion algorithm for the p- p-Wasserstein distance. We are not aware of any sub-linear approximation algorithm for the p- Wasserstein problem that runs in near-quadratic time.
can be used to determine whether an unweighted bipartite graph has perfect matching in faster than n²25 time. While there are several theoretically almost-quadratic algorithms (e.g. (Chen et al., 2022)) for finding a perfect matching, the fastest practical algorithm remains the classic Hopcroft- Karp algorithm Hopcroft and Karp (1973), which runs in (nn-5) time. This reduction suggests a barrier to design- ing efficient and implementable algorithms for computing the p- Wasserstein distance faster than On / time. This computational hardness extends to robust variants as well, including the o-partial y- Wasserstein distance (Chapel et al., 2020) and the A-ROBOT distance (Mukherjee et al., 2021).
Theorem 1.1. Given two discrete distributions a and v defined over point sets A and B. respectively, where |A| + [B] n, a metric d(.,-), and a parameter p 2 1. there exists an O(log n) -approximation algorithm compute the p- Wasserstein distance that runs in On²² log U log A log n) time, where log U is the number of bits required to represent the probability values and A is the spread of A U B.
Previous relative approximation algorithms for the 1- Wasserstein distance rely on embedding the ground met- ric d(.,-) into a tree metric by using hierarchically well- separated tree (HST) (Charikar, 2002; Kleinberg and Tardos, 2002). Embedding a metric into a tree metric introduces distortion -a a form of noise in the edge costs - which is tolerable for = 1 but problematic for p 2 2 due to the increased sensitivity of the p- p-Wasserstein distance to such perturbations. As result, the p- Wasserstein cost computed over the tree metric can deviate significantly from the true cost, yielding unbounded approximation guarantees. To overcome this limitation, we propose a new approach that reduces the impact of distortion by considering multiple tree embeddings. Specifically, we construct p different HST embeddings and, for each pair of points, use the minimum edge cost across these embeddings. We show that this strat- egy yields more stable distance, and prove that the opti- mal p- Wasserstein distance computed using these modified edge costs provides an O apprrxxixinnnnnnnn to the true p- Wasserstein distance.
Interestingly, we show that the (p, kPW distance can be approximated in O(n2/89 time, for any p 2 and fixed constant k, and this runtime is independent of p. Thus, the approximate (p, k)-RPW distance overcomes both major limitations of the p- p-Wasserstein distance: it is robust to noise and supports significantly more scalable algorithms.
Theorem 1.2. Given two discrete distributions js and defined over point sets A and B, respectively, where |A| + [B] n, a metric d(.,)) and parameters 2 1, 8 > 0, and any fixed constant k >> O. there exists a 8-additive approximation algorithm for the (p. kP-P distance that runs in 0n22/53) time.
Our improved algorithm for approximating the (p. k))-PWW distance is motivated by a key observation: augmenting- path- optimal transport algorithms tend to transport most of the mass quickly, with the majority of the run- time spent on routing a small residual fraction (Hopcroft and Karp. 1973; Lahn et al., 2019; Lahn and Raghvendra, 2019). For example, using the classical Hopcroft-Kaut Karp al- gorithm (Hopcroft and Karp, 1973) for bipartite perfect matching, a matching of size (1 - 5)n can be computed in On2 /8) time for any constant 5 > 0, whereas completing the remaining ón edges requires 2.5 time. A similar be- havior also holds for several combinatorial algorithms that approximate the optimal transport cost.
Unlike the case of the Wasserrstein distance, a simple greedy algorithm does not yield an optimal solution under our new distance. Instead, we adopt a primal-dual frame- work to compute the optimal transport plan. In general, primal-dual algorithms perform (n) graph search steps to identify so-called augmenting paths, with each search typically taking On²²) resulting in an overall time complexity of Onn)). To improve efficiency, we leverage the structure of the p tree embeddings and design a dynamic BCP data structure to accelerate the search process. Us- ing this data structure, we conduct each graph search in O(n A log time, where A is the spread, which is the ratio between the largest and smallest nonzero edge costs. This results in a total runtime of (n² log A log n) for com- puting the optimal solution.
In contrast, to approximate the (p, RPW distance e be- tween distributions and v up to an additive error of 5, it suffices to compute a transport plan that moves at least (1 -- 5) mass with a total cost at most + 5. We leverage this relaxation to design a faster approximation algorithm by adapting the LMR algorithm (Lahn al., 2019). Specifi- cally, we find a transport plan that deviates from the optimal in both cost and transported mass by at most d, thereby ap- proximating the RPW distance efficiently. We also provide an implementation of our algorithm. Experiments (included in Appendix C) suggest that our algorithm outperforms the algorithm by Raghvendra et al. (2024), especially for larger values of p and smaller values of 5.
We also provide an implementation of our algorithm and em- pirically analyze the efficiency and accuracy of the solution produced by our algorithm on stochastic data sets.
Additive Approximation: Using the sensitivity of p- Wasserstein distance to noise, we show that any algorithm that approximates the p- Wasserstein distance within an ad- ditive factor of and runs in better than O(n² /8P) time
3

---

Scalable Approximation Algorithms for y- Wassersstti Distance and Its Variants
2. A Tree Based Distance
2.2. Distance Function
In this section, we present a novel distance function based on p independently constructed hierarchically separated trees (HSTs) that approximates the ground metrie d(->). We prove that the p- p-Wasserstein distance with respect to our tree-based distance is, in expectation, an O(log n) ap- proximation of the p- -aasssesteein distance with respect to d(.-)..
Let μ and v be two discrete distributions supported on point sets A and B, with |A| + [B| n, and let d(.,-) be the ground metric. In this section, we introduce an HST-based distance function dr(+;)) .) and show that the D-Wasssrstein distance with (: as the ground distance would be an (log n) apprrxxximtion of the p- Wasserstein distance with .(..)) as the ground distance.
We begin by reviewing HSTs and their key properties as established by Fakcharoenphol et al. (2003) in Section 2.1. In Section 2.2, we formally define our new distance function and demonstrate its approximation guarantees. Then, in Section 2.3, we present dynamic weighted bichromatic closest pair (BCP) data structure for our tree-based distance. Finally, in Section 3, we use this BCP data structure to design a quar-quadratic time algorithm for computing the p- Wasserstein distance with respect to our tree- rree-based distance.
Let T (Ti...., Tp} be a set of p independently con- structed HSTs over the metric space (A U B, d). For any pair of points (a,6) E A x B and each tree index i E [1. p). let cluster Ca b) denote the least common ancestor of and b in the tree Ti and let lev; (a. b) denote the level of the cluster (a b) in the tree, Recall that the diameter of the cluster at level j is at most 2- hence, 2 (a.b) is an up- per bound on d(a, b). We define the tree distance between and b as
2.1. Hierarchically Well-Separated Trees
dr (a, == min {2 1e-levv(a,))). |E[1,)]
Let dX,d) be a metric space. Consider a partition P = (C....., C+} of X into clusters, where each point in X belongs to exactly one cluster in P. For each cluster C' E P. let (C) denote the diameter of C, and let X denote the subset of X contained in C. For each point x E X and radius 5 > 0, define a ball centered at x and radius 8, denoted by B(x, 6), as the points of X that are at a distance no more than from æ, i.c., B(x, 5) := {u E X d(z, u) s 5]. For any integer j 2 0, define k == [Bx,, 2-x] as the number of points of X that are at a distance at most 2- from x.
This definition takes the best (i.e., tightest) estimate among all the p trees an upper bound for cl(a, b). The next lemma relates d7 (a, b) to the ground distance d(a,
Lemma 2.1. For any pair of points (a, E A x B, dr b)P 2 b)a and ((d7(a,b)p) =
(((4logn)d(a,,))
Proof. j = arg pt (e2_evv(ab)),, i.c dr b) = 2 2-lovy(a.b), Recall that C2(a, b) denotes the least common ancestor of a and b in Tj. Then,
Fakcharoenphol et al. (2003) presented randomized algo- rithm that, given metric space (X, with unit diameter and spread A, constructs a hierarchically well-separated tree (HST) T over the points in X in (n²) time. The root of the tree corresponds to all points in X, the children of each cluster C form a partition of the subset Xc, and every point in X forms a leaf cluster in the tree. The level of a cluster C' defined as its distance (in number of edges) from the root: the root is level and the levels increase as we descend the tree. Their construction guarantees two properties:
dr (a, b) (2 (2-leryaabip 2 Dia 2 d(a,
where the second inequality holds from (H1) and the last inequality holds since a, b E Cy(a, b). Let h denote the maximum height of the HSTs in T. For any level j [0, A).
Pr[lev(a, b) = 51 < Pr [lvv: (a. b) S S E [1,p]]
(HI) For any cluster C level j. Diam(Xe) < 2-
II
Pr lev. (a, b) رح 1=1
(H2) For any pair of points u, v E X, the probability that the least common ancestor of the pair (u, u) is at a
d(a, b)P . -Ta,, by
2 2J T x I(u, v) where T(u, v, := Hep + Hy: H Hap here, H3 is the kth Harmonic number.
level at most j is at most d(u.v)
(2-1-1)p
d(a, b)p (211))p-1 l' b.j) --+0+1)
The tree T has a height Olog A). Moreover, every point x E X participates in each of the (log A) clusters that are ancestors of x in T.
where the last inequality holds since ky E [1, n] and I(a, b,j) < 2H, We can bound the expected distortion of the tree distance dy b) by summing over all levels from
4

---

Scalable Approximation Algorithms for Passsesstinn Distance and Its Variants
0 to h as follows.
We present dynamic data structure that maintains this pair and supports insertions and deletions in U and V in O log log A) time, For each tree T; E T and each cluster C' E T, let Uc (resp. Vc) denote the subset of U (resp. V) that participates in C. Let lev(C) denote the level of C in tree T, The following key lemma is instrumental in the design of our data structure.
h
E [|dT(a 6))1 =
Σ
Pr (rllev(a,b) 6) = 12
j=0
d(a, b) 10)..(2H))) (2H,) ".,a,,r,,,,,,)) I(a,b,j), -5
2-1(++1
1=0
Lemma 2.3. Suppose (a',b)) E V X U is the weighted bichromatic closest pair between V and U. Then,
=2- (211)) d(a, 6)"² Σ b(a, j=U
- 2F (,H,,) )01 d(a, b)p x 4Hn
(1)
2 g-lev(C)p - max w(b) mmxx w(a)
= O((4log))) n)? x d(a, b))),
min T,T,CET
BEUC VeVo
where Inequality (1) holds since
Data Structure: For each tree T, T and each cluster C E Ti, define the candidate closest pair for cluster C as the pair (ac, bc) E Vc x UC that minimizes 2lev(C), maXteUe w closest pair of C, our data structure stores two max max-heaps for each cluster C: Heape for points in Ve and Heap' for points in Uc. with weights as keys. The candidate closest pair for C' is simply given by the roots of these heaps.
Σ b(j,b, Σ Hag Haa Hag42 Hag
w(a). To maintain the candidate
j=0
= =Hug + +Hug + Hag + HH2
- -Hug.
Hahh HRg,2
< 4Ha.
We maintain global min-heap GlobalHeap that stores the candidate pair for each cluster C, with their weighted distance as key. From Lemma 2.3, the weighted bichro- matic closest pair corresponds to the minimum entry in GlobalHeap.
Let o* (resp. o7) be the transport plan between μ and v that minimizes the p- Wasserstein cost with respect to the distances d(..)) (resp. .)) For any transport plan o', let wT. (o') denote the p- Wassersteii cost of a' under the tree distances (: Using linearity of pectation, we bound the expected value E[wT,,(0ff (0) by ((44 log w, )*). Combining it with Jensen's 's inequal- ity, we have the following lemma.
To insert (resp. delete) a point b B into (resp. from) the set U, we proceed as follows: for each tree T; T and every cluster C along the path from the leaf containing b to the root of T,, we insert b into (resp. remove b from) the heap Heap'. U If the operation alters the candidate closest pair for any such cluster, we update its entry in the global min-heap GlobalHeap. An identical procedure is used to insert or delete a point in the set V by updating the corresponding heaps Heap'.
Lemma 2.2. [[WT, p (o7)) )
Next, we describe simple dynamic data structure for main- taining a weighted bichromatic closest pair between two subsets U C B and V S A. In Section 3, we use this data structure in a primal-dual algorithm for computing an approximation of the p- p-Wasserstein distance.
Analysis: Building all p HSTs takes O(ma² time. The current bichromatic closest pair is always at the root of GlobalHeap and can be queried in O(1) time. Each point v E A U B belongs in Ollog A) clusters per HST, or Opplog A) clusters overall. Inserting or deleting a point either U or V requires updating the corresponding heaps in all these Opplog A) affected clusters, recomput- ing candidate pairs, and updating GlobalHeap. This takes OOp log log A) time per update.
2.3. Dynamic Bichromatic Closest Pair Data Structure
Let w(u) denote the weight associated with point u. For any pair of points a e A and b E B. define the weighted distance between them
Lemma 2.4. We can build a data structure DD in Opm²) time that stores weighted point sets U and V, returns the bichromatic closest pair in O(1) time, and allows for in- sertion and deletion of points to U U UV in O(p log n log A) time.
(a,
b) == d)( - w(a)- -wb).
Given two sets of points U S. B and V S A, the weighted bichromatic closest pair between U and V is the pair mini- mizing the weighted distance, i.c.,
Remark: The space complexity of the data structure is Oppn log A). This can be reduced to O(pn) by optimizing
arg min (dr(u, r)p - w(()- - w()).
5

---

Scalable Approximation Algorithms for y- p-Wasserstein Distance and Its Variants
the heap storage: instead of storing all points in Ue and Vo, we store only one representative per child cluster (the point with the maximum weight) in the non-leaf clusters.
of Ga is defined as follows: For each (a, A x B, if (a, b) E M, add a directed edge from b to a; if (a, b) M, add directed edge from to b. We set the cost of any edge between a and b, regardless of direction, to be its slack, which is given by d(a, 6)d - y(a) - y(b). We also add a zero-cost edge directed from the source vertex to every vertex in BF and add a zero-cost edge directed from every vertex in Ap to
3. An O (loogn))Apppooo Approximation Algorithm for the Optimal Transport
Given a dynamic weighted bichromatic closest pair (BCP) data structure with a query/update time of ((n), the opti- mal p- Wassersteinn distance between two distributions can be computed in (n) log U) time; here U is the num- ber of bits required to store the probability associated with each point. For completeness, we briefly describe the clas- sical BCP-based algorithm for the assignment problem introduced by (Vaidya, 1989).
The algorithm starts with an empty matching M and ini- tializes all dual weights y(v) + 0 for all E A U B. It performs a sequence of Hungarian Search procedures each of which augments the matching while maintaining dual feasibility. Each Hungarian Search increases the size of M by 1, and after n such steps, we obtain minimum-cost perfect matching.
3.1. Primal-Dual Framework
Hungarian Search Procedure: The Hungarian search pro- cedure consists of the following steps:
Consider complete bipartite graph G(AU Ax B) de- fined on sets A and B. A matching M S A x B is set of vertex-disjoint edges. The size of matching is the number of edges it contains. A vertex u is called free (with respect to M) if no edges of M are incident on v. Let A and B be the set of free vertices in A and B, respectively; we drop the superscript M when the matching is clear from the context. An altermating path is path that alternates between edges in M and those that are not in M, and an augmenting path is one that starts and ends at free vertex. We augment M along such a path P by setting M + M @ P. where edges in P in M are removed and those in P M are added, increasing the size of matching M by 1.
Execute Dijkstra's 's shortest path algorithm with s as the source on GM. For any vertex v, let P. denote the shortest path from s u, and let &, be the total cost of Pe. using slacks as edge costs.
For any b E B, l. < le, set y(b) + y(b) - Lo + br
For any a E A, < &,, set y(a) + y(a) - le + la
Let P be a path obtained by removing and t from P.- Augment M along P.
Note that all steps in the Hungarian search except the execu- tion of Dijkstra's 's algorithm take (n) time. We present an efficient ((n((n)) implementation of Dijkstra's shortest path algorithm on GM-
We associate a dual weight v(v) with each vertex V E AUB, and say that the matching M and dual weights v()) is dual feasible if:
y(a) + v(b) S c(aa 6) v(a) + (b) = d(a, 6)P
for all (a, b) $ M,
(2)
Efficient Dijkstra's shortest path algorithm: Recollect, that Dijkstra's 's algorithm incrementally builds a shortest path free rooted at the source s. Let U C B be the points already added to the shortest path tree, and let V C A be those not yet added. Initially, U = By and V = A the points in Bp. For each b E U, la 0 and w(b) = y(b) l.. For each E V, assign w(a) y(a). Build a BCP structure (U, V) using these weights. In each iteration of the algorithm, select
for all (a, b) E M.
(3)
These are the classical feasibility conditions for the assign- ment problem.
Lemma 3.1. Ifa a perfect matching M along with set of dual weights is dual feasible, then M is a minimum-cost perfect matching.
We now describe how a dynamic bichromatic closest pair data structure can be used to implement the Hungarian algo-
(a, =
arg ei ((dd, + ). (4)
rithm in (n)) time and how it applies to computing the p- Wasserstein distance.
This is equivalent to
s(a', + ty = d(a'' u(a') - (v(b)- - kx)
3.2. Hungarian Algorithm
= d(d''b''p w(a') - -wbb).
To find augmenting paths, the Hungarian algorithm con- structs a residual graph Gy- The vertex set for GM is
Thus, the pair minimizing equation 4 is exactly the BCP between U and V and can be found in 0(1) time.
A U B, source vertex 8 and a sink vertex The edge set
6

---

Scalable Approximation Algorithms for y- Wasserstein Distance and Its Variants
Once the pair (a, b) is selected, the algorithm will remove from V, add it to the shortest path tree, and set la ls + s(a, b). If is free, we have found augmenting path. Otherwise, let b' be the vertex matched with a in M. Set by = fa w(V/) = y(b(- - ly,, and add 1/ to U.
point methods. Thus, the best-known simple, implementable algorithm to determine if there is a perfect matching in G' is by Hopcroft and Karp (1973) and takes ( time. Obtaining a simple, implementable algorithm with sub- n 2.5 execution time remains a major open question. In Section 4.1, we present a reduction (similar to the one pre- sented by (Blanchet et al., 2018)) demonstrating that any algorithm & that produces a ò-additive approximation of the p- Wasserstein distance in sub-n? /8" time can be adapted to achieve a sub-n2.5 algorithm for finding perfect matching, given clever selection for the error parameter d.
Each iteration of Dijkstra's ' deletes one vertex from V and adds at most one vertex to U. The algorithm terminates with an augmenting path in at most n iterations.
Time Complexity: Each iteration of Dijkstra's 's algorithm during the Hungarian Search adds at most one vertex to U and removes one vertex from V. Since there are at most n such iterations, the total number of insertions and deletions to the BCP data structure is (n). Given that each insertion and deletion takes O(plog A log n) time, the total execu- tion time for the Hungarian search is O(pn log A log n.). As the assignment problem requires such searches to compute a perfect matching, the overall complexity is Opnn log A n).
4.1. Reduction
Given an instance of the perfect matching problem on the dense graph G(AU B, E), with |A| [ = n, we convert G to an instance of the p- Wasseesstein problem as follows: Create distributions μ and v by assigning a mass of 1/n for each point of A and B. For each pair (a, b) E Ax
x B, assign distance d(a, b) = 0f (a, b) E E and a distance d(a, b) = 1 otherwise. Now, given this instance of p- p-Wasserstein problem, use any 5-additive approximation algorithm A with running time expressed as T(n, 6, p) to compute a transport plan a and p- p-Wasserstein cost up((r such that Wp(o) I "op(a") + 6, where o" is an optimal transport plan. The optimal p- Wassersteii cost is equal to ((n [M1))/))// where M* is the max-cardinality matching from our original graph G. Therefore, A returns a plan with p- Wasserstein cost most ((n D/n)))))/ '/ + 8. Note that if perfect matching exists in G, then our total transport cost must be (wp(o))² < 8
For the optimal transport problem, when the probabilities as- sociated with each node can be represented using Olog U) bits, (Atkinson and Vaidya, 1995) showed that On log U) Hungarian searches are sufficient to find the optimal trans- port. Therefore the p- Wasserstein distance between two distributions μ and V under the distance (-.,)) -) can be found in O(m= log U log A logn) time, where log U is the num- ber of bits required to represent the probability at every point in the support of the distributions.
Sharathkumar and Agarwal (2012a) demonstrated that the optimal solution for the assignment problem can be computed efficiently by combining the BCP data structure with the cost-scaling framework. Their algorithm performs Ón2322 log(n BCP queries, where log C is the number of bits needed to represent the edge costs. This approach can be integrated with our BCP data structure a straight- forward way to compute an Olog n) approxxxxiannnnnn to the p- Wasserstein distance. The overall execution time of this algorithm is dominated by the construction of HSTs and is Ó(pn² + pn3/2 log" A).
Remark:
We convert the transport plan a into a matching by first scaling all vertex masses and transported masses by a fac- tor of and then ensuring that the plan is integral using the cycle canceling method introduced by Kang and Payor (2015) in linear time. This produces perfect match- ing M' of points in A and B. We then remove all edges from the matching with d(a, b) = 1 (those not originally present in G), which removes at most 52 n edges (i.e., [M'] 2 |M|- - 88- - n)), as all removed edges have a cost of 1. We then match the remaining 8" n unmatched vertices using an augmenting path per unmatched vertex, with each path being found in (n2)) time. Thus, the final running time for the reduction is the sum of the time for our algorithm A and the time for remaining augmentations. Explicitly, we have a running time of T(n, 5, p) + n".
4. Additive Approximation for Wasserstein Distance and Its Variants
In this section, we show that designing fast additive approx- imation algorithms for the p- p-aasserstein distance is chal- lenging due to its sensitivity to noise. In contrast, we show that a robust variant of the p- Wasserstein distance called the (p, k)) RPW distance can be approximated efficiently.
For any additive approximate solver with T(n, 6,p) = n² 8((1=)), we can select 5 = 1/n and obtain an n n,25-(22/2) /2) time algorithm for finding a perfect matching. Our reduction shows that is unlikely that the existing tech- niques can be adapted to obtain 8-additive approximation algorithm for p- Wasserstein distance with an execution time of O(n² ((22//(1-=)) for any constant & > 0.
Consider the problem of determining if an unweighted bi- partite graph G(A U B, E) has perfect matching. All 2.b-n.26 2.5 time algorithms for this problem are either matrix- multiplication-based approaches or sophisticated interior
7

---

Scalable Approximation Algorithms for y- p-aasserstein Distance and Its Variants
4.2. Algorithm for RPW distance
Note that the rounding of masses might create extra mass on the points of A. Thus, when projecting back to the original space, our algorithm removes such excess mass as described in (Lahn et al., 2019).
In this section, we present an additive approximation al- gorithm for the (p, AP-RPW distance, for any p 2 1. For simplicity in presentation, we describe our algorithm assum- ing k = 1 and refer to the distance as the p-RPW distance, noting that similar approach easily extends to any constant > O. For two distributions jl and v, let w" be the p-RPW distance between them, i.c., w" is the smallest value such that Wp.1--- (H. p) I w". Any value w is 6-close to the p-RPW distance if [w- w"I S 6.
Our algorithm initializes ē to an empty transport plan and executes j :- / phases of the GT algorithm on the scaled space, which builds a partial transport plan incrementally in each phase using a primal-dual Hungarian- approach. In each phase, the GT algorithm computes set of aug- menting paths and augments ē along those paths. For each phase 1 S S j", let oj denote the partial transport plan â maintained after the execution of the jth phase. Define ú == = define aj as the amount of mass transported by شj, and let rú = max 1 % first term is an estimate of the untransported mass when we project oj back to the original space, and the second term is an estimate of projected cost. Note that as our algorithm executes the phases (i.e., j increases), the value mi will decrease until the estimated projected cost surpasses the estimated projected untransported mass.
At a high level, given any guess value g of w², our algorithm computes partial transport plan 0g and uses the cost and the amount of transported mass in og to obtain an estimate Wg. which is an upper bound (up to an additive error 8) on the p-RPW distance. We show that, interestingly, when g E [w", - + 2). the computed transport plan oy transports at least 1 wh -d 5 mass with a cost most w" + 8. We use this observation to show that the returned estimate wa would id be a 8-close p- p-RPW distance. Consider the series 9: = T for = L..... : 1. most i away from w'. We execute our algorithm for all gi 's and return the smallest upper bound as a 8-close p-RPW distance.
( )1pp) J: here, the
Note that at least one of the 9, values is at
If ( (" ) < 1 2i for all phases j S J", then our algorithm projects the transport plan āj. to a transport plan og on the original space and returns mng := m.. as an estimate of the p-RPW distance.
Our algorithm executes a modified version of the LMR al- gorithm (Lahn et al., 2019). Given parameter e > 0, the LMR algorithm computes, for each value a E [0, 1]. an E-close a-partial transport plan by first scaling the problem instance to integer supplies, demands, and costs, and then executing a number of phases of Gabow and Tarjan's 's in- leger transportation algorithm (Gabow and Tarjan, 1989) (from now on referred to as the GT algorithm). We present an algorithm that, in similar fashion, first transforms the problem instance to a scaled space, where all supplies, de- mands, and costs are integral, then executes several phases of the GT algorithm with modified stopping condition, and finally projects the transport plan back to the original problem space. Next, we describe our algorithm in detail.
Otherwise, suppose at an intermediate phase j S j", (1)==2)1// exceeds the mass 1 - #1 for the first the cost ( 81 (")-2))// 2 1 a: and, for all k < j. time, i.e., ( < <1 os For any a E (a)=1, aj), let 6th denote the a- a-partial transport plan computed during the execution of the jth phase of the GT algorithm, i.e., the transport plan obtained after augmenting 0j-1 by a - aj-1 mass. Let
a' be the value in (a)-1, such that (0(0"")=)))
a", As discussed by Phatak et al. (2022) and Raghven- 51 dra et al. (2024), the value a" can be computed by tak- ing the intersection point of a line segment and with the curve = (1 - ) in constant time. Our algorithm then projects the transport plan ga to a transport plan dg in the original space and returns 0g as well as the value my = max {1 {1- 2- "eld""""ee)///" estimate of the p-RPW distance given the guess g. This completes the description of our algorithm.
Given a guess value g, our algorithm lifts the problem in- stance to scaled space as follows: Let $1 2pm and 82 = 0²5 For any point a E A (resp. b E B), define the scaled mass ji(a)= := [s| uuaa] (resp. ((b)= := (bb))- For any pair of points (a, 6) E A x B, we set their distance to be à((abb) b)- = [d(a,b)² For any transport plan à in the lifted space, let ( denote the 1 1-Wasserstein cost of ò under d(.,-) distances.
2p
= 1 - a as an
Efficiency: For each guess value g, our algorithm executes at most j = 5 phases of the GT algorithm, each in (n?)) time. Projecting the input instance to the lifted space and converting the computed transport plan back to the origi- nal space also takes O(n?) time. Consequently, our algo- rithm processes each guess value in 0 (%) time; summing across all O(1/8) guess values, the running time of our al- gorithm would be 0( (5))
At any point in the algorithm, we can also project our cur- rent transport plan back to our original problem space as follows: given a lifted transport plan o, we project to the original problem space by defining a transport plan dg. where (y(a,b) = (a, (a,b)/s: for all pairs (a,b) A x B.
8

---

Scalable Approximation Algorithms for y- paassesstenn Distance and Its Variants
10'
2,2
Data Structure Operations
p=2
2.75
n=20000
--30000
Approxination Facter
p=4
2.50
D-00000
2,0
P-5
10f-
an (3/2)
Apgrroximaatio
2.5
D-40000
n=400
2.00
1.8
1.75
1.6
10a-
150
1:25
10000
20000
30000
40000
50000
10
1.00
2
3
4
(a) Approx. factor (Uniform)
(b) Data structure operations (Uniform)
(c) Effecc the number HSTs (Uniform)
10a
2.2
Data Seructure Operations
p
Approxinnaton Fector
Average Dierortion
0-4
2.0
5
(3/2)
10f
an
1.8
1.6
p-3
102
p
1.4
10000
30000
50000
10a
2 10
10
10
10
(d) Approx. factor (Normal)
(e) Data structure operations (Normal)
(f) Expected tree distance distortion
Figure 1. Our experimental results, (a) and (d): the approximation factor of our computed matching as function of n. (b), (e): the number of calls to the data structure operations by the Hungarian algorithm, in log-log scale, (c): the approximation factor of our computed matching for p 4 as varying the number of HSTs, and (f) the average distortion of the tree distances as function n, in log-log scale.
For each guess gr, we show that our algorithm returns a value that is guaranteed be least w" 6/ 6/16; here w" is the optimal (p, 1) RPW distance. Furthermore, show that for any guess g [w" w- + 5/4), our algo- rithm returns a value that is at most w² 5. See Appendix B for a detailed analysis.
Correctness:
change as a function of n, while increases exponentially as p increases, which is in line with our result in Lemma 2.1.
Algorithmm's 's Accuracy: For given value of and p and both datasets, we executed our algorithm on n samples A and B from the dataset and computed the ratio of our matching cost to the optimal matching cost. For each value of n and we report the average ratio over 15 executions. As shown in Figure (a) and (d), as n grows, the approximation factors slightly increase; however, the approximation factors do change function of p. We also measured the accuracy of our algorithm by varying the number of HSTs used in our construction. For 4, we varied the number of HSTs m from to 6 and measured the approximation ratio of the computed matching. As shown in Figure I(c), introducing more than HSTs seems to have an insignificant effect on the algorithm's accuracy.
5. Experiments
In this section, we present the empirical analysis of both our tree-based distance from Section 2.2 as well as the algorithm from (Section 3). All computations puter with a 2.6 GHz 6-Core Intel Core i7 CPU and 16 GB RAM, using a single calculation thread. The code is avail- able at https://githuub.com www/poooaass///ssaaooooo Wasserstenn.
performed
com-
Datasets: We perform experiments on two datasets, namely (i) samples drawn from the uniform distribution inside the 10-dimensional unit hypereube (Uniform dataset), and (ii) samples drawn from truncated normal distribution inside the O0-dimensional unit hypercube (Normal dataset).
Algorithm's Efficiency: We analyze the complexity of algorithm by counting the number of calls made to the data structure. As we can see in Figure ((b) and (e), our experi- ments suggest that the number of data structure operations grows faster than linear but sub-n 1.5 which is significantly better than the worst-case upper bound of O(n²). Inter- estingly, this empirical behavior is reminiscent of the per- formance guarantees established for the scaling algorithm of (Sharathkumar and Agarwal, 2012a).
Results: We present our results in three parts:
Tree Distances Accuracy: For two sets of points from the Uniform dataset, we compute the average distortion of the tree distance across all pairs of points. As depicted in Figure 1(f,, the average distortion does not significantly
9

---

Scalable Approximation Algorithms for p-Wasserstein Distance and Its Variants
Acknowledgment
In Proceedings of the International Congress of Mathe- maticians: Rio de Janeiro 2018, pages 3287-3318. World Scientific, 2018.
This research was supported by NSF grants CCF 2314753 and CCF 2223871. We would also like to acknowledge Ad- vanced Research Computing (ARC) Virginia Tech well as A Root Cluster (ARC) at North Carolina State University, which provided us with the computational resources used to run the experiments.
David S Atkinson and Pravin M. Vaidya. Using geometry to solve the transportation problem in the plane. Algorith- mica, 13(5):442 461, 1995.
Arturs Backurs, Yihe Dong, Piotr Indyk, Ilya Razenshteyn, and Tal Wagner. Scalable nearest neighbor search for optimal transport. In Proc. 37th International Conference on Machine Learning, pages 497-506, 2020.
Impact Statement
This paper presents work whose goal advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Jose Blanchet, Arun Jambulapati, Carson Kent, and Aaron Sidford. Towards optimal running times for optimal trans- port. arXiv preprint 1r11:1181000771,,2 2018.
References
Timothy M Chan. Dynamic generalized closest pair: Revis- iting eppstein's 's technique. In Symposium on Simplicity Algorithms, pages 33-37. SIAM, 2020.
Pankaj K. Agarwal and R. Sharathkumar. Approximation algorithms for bipartite matching with metric and geo- metric costs. In Sympos. on Theory of Computing, pages 555-564, 2014.
Wanxing Chang, Shi, and Jingya Wang. Csot: Cur- riculum and structure-awand optimal transport for learn- ing with noisy labels. ArXiv preprint arXiv:2312.06221,, 2023.
Pankaj K. Agarwal and Kasturi R. Varadarajan. A near- linear fonttant-factor approximation for Euclidean bipar- tite matching? In 20th Internat. Sympos. on Computa- tional Geometry, pages 247-252, 2004.
Laetitia Chapel, Mokhtar Z. Alaya, and Gilles Gasso. Partial optimal tranport with applications on positive-unlabeled learning. Advances in Neural Information Processing Systems, 33:2903- 2913, 2020.
Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirza- dian, and Rachita Sowle. An improved e-approximation algorithm for geometric bipartite matching. In 18th Scan- dinavian Sympos. and Workshops on Algorithum Theory, 2022.
Moses S Charikar. Similarity estimation techniques from rounding algorithms. In Proceedings of the thiry-fourth annual ACM symposium on Theory of computing, pages 380-388, 2002.
Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Keegan Yao. Efficient approximation algorithm for computing wasserstein barycenter under euclidean metric. In Proceedings of the 2025 Annual ACM-SIAM Sympo- sium on Discrete Algorithms (SODA), pages 1809-4826, SIAM, 2025.
Li Chen, Rasmus Kyng, Yang P Liu, Richard Peng, Maxim- ilian Probst Gutenberg, and Sushant Sachdeva. Maximum flow and minimum-cost flow in almost-linear time. In Proc. 63rd IEEE Annual Sympos. Foundations of Com- puter Science, pages 612 623.. IEEE, 2022.
Jason Altschuler, Jonathan Niles Weed, and Philippe Rigol- let. Near-linear time approximation algorithms for opti- mal transport via Sinkhorn iteration. Advances in neural information processing systems, 30, 2017.
Ching Yao Chuang. Devon Hjelm, Xin Wang. Vibhav Vineet, Neel Joshi, Antonio Torralba, Stefanie Jegelka, and Yale Song. Robust contrastive learning against noisy views. In Proc. IEEE/CVF Conference on Computer Vi- sion and Pattern Recognition, pages 16670-16681, 2022.
Jason Altschuler, Francis Bach, Alessandro Rudi, and Jonathan Neees-Weed. Massively scalable sinkhorn dis- tances via the nyström method. Advances in neural infor- mation processing systems, 32, 2019,
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, pages 2292-2300, 2013.
David Alvarez-Melis and Tommi Jaakkola. Gromov- Wasserstein alignment of word embedding spaces. arXiv preprint arXiv:
Frederik Michel Dekking, Cornelis Kraaikamp, Hen- drik Paul Lopuhaä, and Ludolf Erwin Meester. A Modern Introduction to Probability and Statistics: Understanding why and how. Springer Science & Business Media, 2006.
2018
Alexandr Andoni, Piotr Indyk, and Ilya Razenshteyn. Ap- proximate nearest neighbor search in high dimensions.
10

---

Scalable Approximation Algorithms for paassersttin Distance and Its Variants
Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complex- ity by accelerated gradient descent is better than by sinkhorn's ss algorithm. In International conference on machine learning, pages 1367-1376. PMLR, 2018.
International Symposium on Computational Geomelry, 2019.
Nathaniel Lahn and Sharath Raghvendra. An O time approximation algorithm for RMS matching a plane, In Proc. 2021 ACM- Sympos. on Discrete Algorithms, pages 869-888. SIAM, 2021.
David Eppstein. Dynamic Euclidean minimum spanning trees and extrema of binary functions. Discrete & Com- putational Geometry, 13:111-1221 1995.
Nathamiel Lahn, Deepika Mulchandami, and Sharath Raghvendra. A graph theoretic additive approximation of optimal transport. Advances Neural Information Processing Systems, 32, 2019.
Jittat Fakcharoenphol, Satish Rao, and Kunal Talwar. A tight bound on approximating arbitrary metries by tree metries. In Proc. 35th Annual ACM Symposium on Theory of Computing, pages 448-455, 2003.
Nathaniel Lahn, Sharath Raghvendra, and Kaiyi Zhang. A combinatorial algorithm for approximating the opti- mal transport in the parallel and mpe settings. Advances in Neural Information Processing Systems, 36:21675- 21686, 2023.
Nicolas Fournier and Arnaud Guillin. On the rate of conver- gence in Wasserstein distance of the empirical measure. Probability Theory and Related Fields, 162(3):707-738 2015.
Zhengfeng Lai. Chao Wang, Sen-ching Cheung, and Chen- Nee Chuah. Sar: Self-adaptive refinement on pseudo labels for multiclass-imbalaneed semi-supervised learn- ing. In Proc. IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 4091-4100, 2022.
Kyle Fox and Jiashuai Lu. A near-linear time approxima- tion scheme for geometric transportation with arbitrary supplies and spread. In Proc. 36th Anmual Sympos. on Computational Geometry, pages 45:1-45:18, 2020.
Harold N. Gabow and Robert F. Tarjan. Faster scaling algo- rithms for network problems. SIAM Journal of Comput- ing, 1013:1003)) October 1989. ISSN 0097-5397.
Tianyi Lin, Nhat Ho, and Michael Jordan. On the effi- ciency of the sinkhorn and greenkhorn algorithms and their acceleration for optimal transport. arXiv preprint arXiv: 1:X.v19906.013311 2019.
Harold N. Gabow and Robert E. Tarjan. Faster scaling algorithms for general graph-matching problems. Journal ACM, 38(4):815-853, 1991.
Debarghya Mukherjee, Aritra Guha, Justin Solomon, Yuekai Sun, and Mikhail Yurochkin. Outlier-robust opti- mal transport. In Internat. Conference Machine Learn- ing. pages 7850-7860. PMLR, 2021.
John Hopcroft and Richard M Karp. An 5/2 algorithm for maximum matchings in bipartite graphs. SIAM Jour- nal of computing, 2(4):225-231, 1973.
James Orlin. A faster strongly polynomial minimum cost flow algorithm. In Proc. 20th Annual ACM Sympos. The- ory of Computing, pages 377-387, 1988.
Piotr Indyk. A near linear time constant factor approxi- mation for euclidean bichromatic matching (cost). In Proc. 18th Annual ACM-SIAM Sympos. on Discrete Al- gorithmss, pages 39-42. 2007. URL http: acm. 1283388.
Abhijeet Phatak, Sharath Raghvendra, Chittaranjan Tripathy, and Kaiyi Zhang. Computing all optimal partial trans- ports. Proc. 11th Internat. Conference on Learning Representations, 2022.
org/ stg/eitation... wwwmiidniiiimfe
Arun Jambulapati, Aaron Sidford, and Kevin Tian. A direct ((1/c) iteration parallel algorithm for optimal transport. arXiv preprint arXiv: 006600018, 2019.
Y. V. Prokhorov. Convergence of random processes and limit theorems in probability theory. Teor. Veroyatnost. i Primenen., 1:177-238, 1956. [English translation: The- ory Probab. Appl. (1956), 157-214].
Donggu Kang and James Payor. Flow rounding. arXiv preprint arXiv: 100.::1157..81111 0899, 2015.
Sharath Raghvendra, Pouyan Shirzadian, and Kaiyi Zhang A new robust partial p Wassersteein based metric for com- paring distributions. In 4/st Internat. Conference on Machine Learning, 2024.
Jon Kleinberg and Eva Tardos. Approximation algorithms for classification problems with pairwise relationships: Metric labeling and markov random fields. Journal of the ACM (JACM), 49(5):616-63 639, 2002.
R Sharathkumar and Pankaj K Agarwal. Algorithms for the transportation problem in geometric settings. In Proc. 23rd Annual SSAM Sympos. on Discrete Algorithms, pages 306-317, 2012a.
Nathaniel Lahn and Sharath Raghvendra. A weighted ap- proach to the maximum cardinality bipartite matching problem with applications in geometric settings. In 35th

---

Scalable Approximation Algorithms for Paassersttin Distance and Its Variants
R. Sharathkumar and Pankaj K. Agarwal. A near-linear time approximation algorithm for geometric bipartite match- ing. In Proc. 44th ACM Annual Sympos. Theory of Computing, pages 385-394, 2012b.
Jonah Sherman. Generalized preconditioning and undi- rected minimum-cost flow. In Proc. 28th Annual ACM- SIAM Sympos. Discrete Algo., pages 772-780, 2017.
Pravin Vaidya. Geometry helps matching. SIAM
Journal of Computing, 18(6):1001-1255, 1989.
Mikhail Yurochkin, Sebastian Claici, Edward Chien, Farzaneh Mirzazadeh, and Justin M Solomon. Hierarchi- optimal transport for document representation. Ad- vances in Neural Information Processing Systems, 32,
2019.
12

---

Scalable Approximation Algorithms for paasersstiin Distance and Its Variants
A. Missing Proofs and Details of Section 2 Lemma 2.2. ((7}) O((log n));)("))
Proof. From the definition of the p- Wasserstein distance and since of is an optimal transport plan with respect to the tree distances,
1/p
1/p
E UT
( ((
=F
Σ (a,))AAB
Ha
SR
Σ
(a
(a,b)EAxB
1/p
E
o'(a, b)P (a.))
Σ
(5)
where the last inequality holds by Jensen's inequality (Dekking et al., 2006). Using the linearity of expectation and Lemma 2.1, from Equation (5) and for some constants c > 0,
d'(a,b)dr(aab) 1)
1/v
1/2
WT, x {o7)
Σ
Σ
b)
(a.b)EAXB
(o,))
1/P
1/p
Σ
ba( (4 d(a b)p
(cglogg)
Σ (a.b)EAxB
o d' d(a, 6)p
(a.))EAXB
(c2 log n))p((")
Lemma 2.3. Suppose (a" .b*) E V x is the weighted bichromatic closest pair between V and U. Then,
dr dr,n(a",
2 lev(C) max w(b) max w(a)
min CET
bEllC
aEVo
Proof. For any cluster C' at aevel of tree T, E T, define d(C) 2 2- 2 Define the weighted C-distance of
and b as dc. the weighted C-distance. Consider any tree T; E T, any cluster C' E T, a level j, and any pair (a, b) Ve Uc. From the definition of the tree distances,
d(C)) w(a) w(b) and the weighted C-closest pair the pair (a", (a",b") bb)) Vo x Uc minimizing
(a, b) (a, (() w(b) < 0(2=(((0)))))) (2-lovs(a.b)) w(a) w(b) s d(C)² w(a) w(b).
Furthermore, if C' is the cluster determining the tree distance of and b (i.e., the tree distance of and is determined by the tree T; and C' is the least common ancestor of and b in T.). then
(a, b) d(C)p w(a) w(b)
Therefore, for any pair of points (a, 6) V x U, if C2,b denotes the set of all clusters of T3 containing both a and 6,
dt d b) min
{d(C)- w(a) w(b)).
(6)
CEC,
Consequently, if (a*, b*) V x U is the weighted bichromatic closest pair, then
dr (a", b²) min (o.b)eV xU (a, =
min (U,))EVX TIET. CEC,, de b)
min T,ET, CET, (au
da min
min T.T. CET; (a. XUe
de
(7)
13

---

Scalable Approximation Algorithms for p-aasserstein Distance and Its Variants
In other words, to find the weighted bichromatic closest pair, one can compute the weighted closest pair inside each cluster in all trees in T and find the closest one among such pairs. Finally, from Equation (6),
min dc,w(a,b) b) (ox))eeexUe
min (o.
{d(C)p ww() w(6)
d(C)²
-mxx w(a) max w(b). BEUl
(8)
.VVo
Therefore, the weighted bichromatic closest pair for the cluster C' would be the pair (a, b) E Vc x UC where a (resp. b) has the largest dual weight among all points in Ve (resp. Uc). Combining Equations (7) and (8),
{
{
faccy
d b²)
min
(
min
de dc,w(a,b)
min T,T,,CE,,
max w(a) max w(b) j. beUc
T,ET, OeT: (al
aeVc
as claimed.
B. Missing Proofs and Details of Section 4.2
In this section, we show that the value w returned by our algorithm from Section 4.2 is a ô-additive approximation of the p p-RPW distance. Recall that our algorithm computed an estimate myg for any guess value gi i6 for all i E [1,5/4]. In Zn where w" is the p-RPW distance (Lemma B. We then show that for guess value g" [w" 222 + 5/4), the computed estimate m,- is at most w" + 5 (Lemma B.3). Combining the two lemmas, by picking the minimum value among all values in {mi.... ../.)) our algorithm correctly returns a 6-additive approximation of the p-RPW distance.
our proof, we first show that for any guess value g. the computed estimate my is at least w 8
Lemma B.1. For any value E (0, 1). let m, denote the value returned by our algorithm given the guess value Then, i 2 w²
& 2p
Proof of Lemma B.1. We begin by bounding the change in the cost and the transported mass of a transport plan ā when projecting to the original space.
Lemma B.2. For any value g & (0, 1), lel â denote an o-partial transport plan in the scaled space, and let a denote its projection to the original space. Then, a transports at least a
2p mass and has a cost wp(o) S [bb(a)s22 1/p
Proof. Recall that our algorithm constructs the scaled space by defining j(aa) = [[s [aa)] (resp. ((b) = 81 v(b)) ) for each point a E (resp. b E B), where 81 = _. Therefore, for any point a E A, the transport plan transports at most ji(a)= S $1 p(a) 1 mass from a to the points of B, i.e., for any point a A,
b(a
((a)
1 81
< u(a) +
81
81
where the first inequality holds since a partial transport plan. Recall that our algorithm projects back to the original
space by first defining transport plan transporting (a, 6) = (a, b)/ si mass between each pair (a, 6) E A x B. From Equation (9), for any point a E A, the transport plan transports up s1 2pm excess mass, which our algorithm then adjusts by arbitrarily removing masses the edges incident on Therefore,
ΣΣ ((a,b) b) BEAbEB
Σ a
o(a,b)
1 $1
Σ ((.)))
(10)
81
81
81
81
2p
BEB
Next, we show that the cost of a is at most
2(d)
Recall that our algorithm defines the scaled distance of each pair
81
(a, b) E A x B as b(a, = d)a for 82
266
Since d(a.)?
S d(a, b) and (a, b) S (a, b)/s1.
2p 82
o(a. b) syd(a, b) w()
(wp(o))p=
Σ
oaa b) b)
Σ
81
(a.))
(a.0)EAx
14

---

Scalable Approximation Algorithms for p-Wasserstein Distance and Its Variants
For any transport plan (in the original space), let a:()) denote the amount of mass transported by As shown in Lemma 2.2 in (Raghvendra al., 2024), if we denotes the p-RPW distance between μ and then w I max a()), wp(a)) Given a guess value [0, 1J. recall that our algorithm computes an á-partial transport plan & in the scaled space and returns a value mg = max 81
) and projected transport plan og. From Lemma B.2, a(rr) 2
2p
and up(o,) & (2) 1/P Hence.
m max 1 â
( ((()) [0(0))82) 1/p)
2 max {1 a(0g)
2p" u,loo)) 2 w²
8
$1 $1
2p
This completes the proof of Lemma B.I.
Lemma B.3. For any value g" with g" E [w", w² + 4). suppose thn- denotes the returned estimate of the p-RPW for the guess g" Then, mg- "0 8.
Proof of Lemma B.3. Let 6g denote the transport plan in the scaled space computed by our algorithm for the guess value g", and let 0g- denote the projection of āg- to the original space. To prove Lemma B.3, we first show in Lemma B.5 that the transport plan og. transports at least $1 (1 wt 5) mass, i.c., if â denotes the amount of mass transported by o g' then
S ~" + 6.
(11)
We then show in Lemma B.6 that in the scaled space,
(20r(ōr)))
1/P
S w + 8.
(12)
$1
Combining Equations (11) and (12), mg. S w + 6, claimed.
The following lemma helps in proving Lemmas B.5 and B.6.
Lemma B.4. Let a denote optimal (1 w²) >-partial transport plan the original space, and let 5* denote its projection
to the scaled space. Then, w(+")= < s1(w") + $1. Furthermore, the total amount of mass transported by O* is at least s1 w" 2p
Proof. From the definition of the scaled distances as well the projection of transport plans to the scaled space,
( (cla,,))) b)P
81 (wmp(a"))
= Σ
)'(a, b) å(a,b) s
Σ
$1a" (a,
+ +1)
+ 81 (1 w
$2
82
(a
(xx
81 (w")) p
81.
$2
For any point b E B, the transport plan a transports at most v(b) $ 0())+1 mass from to the points of A, i.e., for any point
b E B,
Σ (1 b) S slu(b) S ((b) +
(13)
DEA
After defining scaled transport plan & transporting ba) s b) mass between each pair (a, b) E A x B, from Equation (13), for any point E B, the transport plan ē transports up to unit of excess mass, which our algorithm then adjusts by arbitrarily removing masses on the edges incident on b. Therefore,
8 2p
ΣΣ o(a,b) 2
Σ DEB
Σ aEA
(1) -1 =
b(a,b) - n = =81(1= w²) )- -n 2 s1 w Ax
Σ
BEB aEA
15

---

Scalable Approximation Algorithms for yo g-Wasserstein Distance and Its Variants
Lemma B.5. For the guess 9" with 0 S g² -w"S let az. denote the amount of mass transported by the transport plan og. at the end of phase jr. Then, 1 áz S wt + 8. $1
g
Proof. We begin by introducing the notations necessary to prove this lemma. For any augmenting path P, let the net-cost of P. denoted by o(P), be the rate of change in the cost of the transport plan when augmented along P. i.e. augmenting the transport plan o along the path P by 3 mass would increase the cost of a by 8 x Φ(P). Each phase j of Gabow and Tarjan's 's algorithm is associated with a net-cost Φy and in phase j. it finds and augments the maintained transport plan along a set of augmenting paths, each with a net-cost oj. is shown that
(P1) in each phase j. the net-cost value 6y increases by least one, i.e., Φ3 2 Φ + 1 (Gabow and Tarjan, 1991).
For transport plan in the original space, we project to the scaled space follows. We first set ò (a. b) s1 for all pairs (a, b) E A x B, and, for each point à E B such that the amount of mass transported out of according to ā is more than p(6), we arbitrarily remove mass transportation on the edges incident to it.
Recall that our algorithm stops after executing j = Sp phases of the Gabow and Tarjan's 's algorithm. We show that the transport plan o.- computed at the end of phase j" transports least 81 (1 =-- -)) mass to conclude the lemma statement.
Suppose o-- transports s1 =² 2Φ 3 mass, for some > (otherwise, the lemma statement holds trivially). Suppose, instead of stopping carly, our algorithm would have continued executing the phases until transporting a total of (1 w² & mass, and let P = (Pi, .2..... ......PP, Pk)) denote the sequence of augmenting paths computed by the phases to transport the remaining 3 mass. Note that by property (P1), all augmenting paths in P have net-costs more than j z=. Therefore, if B; denotes the amount of mass augmented along each augmenting path P; and o.. denotes the $1 (1 w" 2p -partial transport plan obtained after augmenting oz- along the paths in P, then
8p
8p
8p3
w w6ww)) ) (6,
oPP.)..
(14)
Combined with Lemma B.4 and the property that each intermediate transport plan is 51 -close as shown in (Phatak et al., 2022),
8p3 8²
81 (1w""
S (6 - - S (6 S w(0") + S
+ 231
(15)
82
where the first inequality holds from Equation (14), the third inequality holds since ā is a 81 close (1(1 - .² transport plan and o* is a transport plan that transports at least 81 - w² 4 2p ) mass, and the last inequality holds from Lemma B.4. Recall that g" 2 w". Plugging in the value $2
5 ppartill
2p
"(g")" into Equation (15), 2p
8 < $16² (2p(w")) +2) s $18 x
(1
< 51 2
(16)
8p 8(g^)p
4
4p
In other words, the transport plan ãz- computed at the end of phase j transports at least aj+ 2 s/ - w²
6
616
2p
2
mass, and 0y" < w² + 6.
Lemma B.6. For the guess g" with 0 S 9" w² &. any transport plan â computed by our algorithm in the scaled space that transports a mass a E [st w² 8), 81 (1 w² 2p satisfies
(22((())
1/p
w 8.
51
Proof. Recall that, by the properties of Gabow and Tarjan's 's algorithm, the transport plan would be a si sclose a-partial transport plan. If 0% denotes the a- e-partial OT plan in the scaled space, then using Lemma B.4,
s1(w")
((6) < (6x2) + +1 S 6(60)) + s1 -
+ 251.
82
16

---

Scalable Approximation Algorithms for p-Wasserstein Distance and Its Variants
n= 1000
2000
5000
our daparttes
digritim UMA 11W
UVBt bused 1/W
UMB NFW
122
1200
82-
Time (becord21
5000
Time (oacord41
4C00
13
203
620
2000
22
Figure 2. Comparing the running time of RPW algorithm (blue lines) with the algorithm of (Raghvendra al., 2024) (orange lines)
increases.
Therefore,
$2
(6)
5
$2
(s(W")²
+ 281
S (w=)p= + 282 = -w=)²++ + 8(gg)²
S (w + ō(w* + 88/4)²
S (w" + 8)y,
81
81
where the last inequality holds since
(w+ 55 = (w*)p++ +
P 1op0-110
8
P--22 282 2
+... 5p
2 (w") } +
w
P
w
+
1
P
w
2p-282
+ 8"
p
4p
p
2 (w")²
+
p
1
wP8 +
w
P
w
....
+
4
5(w" + 8/4)
2 (w")p +
p
C. Additional Experimental Results
In this section, we present an empirical comparison between our RPW algorithm presented in Section 4.2 and the current state-of-the-art RPW algorithm presented in (Raghvendra et al., 2024)(L We demonstrate that the algorithm is both implementable and efficient. All experiments were performed North Carolina State University's 's computing cluster (ARC), on a single calculation CPU thread. Implementation can be found at https::////. //tpps:/ggt /ttph/gtttt/tt///// com/ saarinenemma/ FasterRPR.
Results: We perform experiments on samples drawn from the uniform distribution inside the unit square. For two sets of n points, one sampled from the left half of the unit square and the second drawn from the right, we execute both LMR-RPW and our RPW algorithms to find the p-RPW distance between the datasets. We measure efficiency by comparing runtime in seconds. We present two strengths of our algorithm:
Dependence on p: For 8 0.1 and any n E (1000, 2000, 5000), we compute the p-RPW distance between samples from the left half and samples from the right half of the unit square for increasing values of As shown in Figure 2, for each sample size, the running time of our algorithm grows much slower as increases compared to L.RRRPWW, and we observe that for 4, our algorithm computes an approximate p-RPW significantly faster than LMR-RPW.
Dependence on d: For n = 1000 and p = 3, we compute the p-RPW distance between two sets of n samples drawn from each half of the unit square for 8 E (0.1, 0.08, 0.06, 0.04). As shown in Figure 3, as the additive error decreases (i.e., the accuracy increases), the runtime of the LMR-RPW algorithm increases at a greater rate, and our algorithm outperforms the LMR-RPW for all values of 5.
17

---

Scalable Approximation Algorithms for p-aaseerttiin Distance and Its Variants
10000, p
Nggoithm
Time (soconds)
20
delta
Figure 3. Comparing the running time of our RPW algorithm (blue lines) with the algorithm of (Raghvendra et al., 2024) (orange lines) as 8 inereases.
18


Reviews:
[{"review_id": "XJkqv6k6P2", "reviewer": "Reviewer_wEzT", "summary": "This paper is concerned with developing a new algorithm for estimating p-Wasserstein distances. Notably, this algorithm enables approximating  the $p$-Wasserstein distance between distributions supported on $n$ atoms up to a multiplicative factor which scales as $O(log(n)$ (in expectation), i.e. the approximation satisfies $0\\leq \\mathbb E[A(\\mu,\\nu)]\\leq C\\log(n) W_p(\\mu,\\nu)$. The approximation can be computed in time scaling as $O(n^2\\log n \\log U\\log \\Delta)$, where $U,\\Delta$ are problem dependent parameters. The approach is based on using a number of hierarchically well-separated trees which are constructed independently and computing a \"tree-based distance\" for each instance. This distance can be computed efficiently using a dynamic bichromatic closest pair (BCP) data structure which can be used to compute Wasserstein distances exactly in $O(n^2\\Phi(n)\\log(U))$ time, where $\\Phi(n)$ is the query/update time for the underlying data structure (in this implementation, $O(\\log(\\Delta)\\log(n))$). \n\nIt is then argued that a $\\delta$-additive approximation of the $p$-Wasserstein distance with an execution time of $O(n^2/\\delta^{p(1-\\epsilon)})$ is unlikely to exist. On the other hand, by adapting a version of the LMR algorithm algorithm, it is shown that a $\\delta$-additive approximation of the so-called RPW problem can be obtained in $O(n^2/\\delta^3)$ time. The paper concludes with some experimental validations of these findings.\n\n## update after rebuttal\n\nThe author's rebuttal addressed my primary question with this work, I have updated my score in consequence.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough review and constructive feedback. We answer the main concern raised by the reviewer below.\n\n>The importance of $O(\\log n)$ approximation factor.\n\n**Response:** A tree-embedding-based $O(\\log n)$-approximation for the $1$-Wasserstein distance originally developed by [1, 2] has been effectively utilized in various applications. These applications include the design of the FlowTree algorithm [3], $k$-NN and LSH data structures [4, 5], $1$-Wasserstein barycenter computation [6, 7], streaming algorithms for $1$-Wasserstein [8], design of $(1+\\varepsilon)$-relative approximation algorithms [9, 10] for the $1$-Wasserstein distance, improving the accuracy of additive approximations such as Sinkhorn for the $1$-Wasserstein distance [11], and defining related metrics such as the tree-sliced Wasserstein distance as a proxy for the $1$-Wasserstein distance [12]. For many of these applications, a bound of $O(\\log n)$ on the approximation factor is critical.\n\n\nDespite considerable efforts including some lower bounds [13], there are very few known approximation algorithms for the $p$-Wasserstein distance for $p > 1$ [14, 15]. In this paper, we introduce the **first $O(\\log n)$-approximation algorithm** for computing the $p$-Wasserstein distance for any fixed $p \\ge 1$. Extending our algorithm to applications such as $k$-NN under the $2$-Wasserstein distance or boosting the accuracy to $(1+\\varepsilon)$-approximation remains an important future direction of research.\n\n----\n\n[1] M. S. Charikar. \"Similarity estimation techniques from rounding algorithms\". STOC, 2002.\n\n\n[2] J. Kleinberg and E. Tardos. \"Approximation algorithms for classification problems with pairwise relationships: Metric labeling and Markov random fields.\" JACM, 2002.\n\n[3] A. Backurs, Y. Dong, P. Indyk, I. Razenshteyn, and T. Wagner. \"Scalable nearest neighbor search for optimal transport.\" ICML, 2020.\n\n[4] P. Indyk and N. Thaper. \"Fast image retrieval via embeddings\". International Workshop on Statistical and Computational Theories of Vision, 2003.\n\n[5] T. Liu, A. Moore, K. Yang, and A. Gray. \"An investigation of practical approximate nearest neighbor algorithms.\" NeurIPS 2004.\n\n[6] T. Le, V. Huynh, N. Ho, D. Phung, and M. Yamada. \"Tree-Wasserstein barycenter for large-scale multilevel clustering and scalable Bayes.\" ArXiv:1910.04483, 2019.\n\n[7] P. K. Agarwal, S. Raghvendra, P. Shirzadian, and K. Yao. \"Efficient Approximation Algorithm for Computing Wasserstein Barycenter under Euclidean Metric.\" SODA, 2025.\n\n[8] X. Chen, R. Jayaram, A. Levi, and E. Waingarten. \"New streaming algorithms for high dimensional EMD and MST.\" STOC, 2022.\n\n[9] J. Sherman. \"Generalized preconditioning and undirected minimum-cost flow.\" SODA, 2017.\n\n[10] P. K. Agarwal, S. Raghvendra, P. Shirzadian, and K. Yao. \"Fast and accurate approximations of the optimal transport in semi-discrete and discrete settings.\" SODA, 2024.\n\n[11] P. K. Agarwal, S. Raghvendra, P. Shirzadian, and R. Sowle. \"A higher precision algorithm for computing the 1-Wasserstein distance.\" ICLR, 2023.\n\n\n[12] T. Le, M. Yamada, K. Fukumizu, and M. Cuturi. \"Tree-sliced variants of Wasserstein distances.\" NeurIPS 2019.\n\n[13] A. Andoni, A. Naor, and O. Neiman. \"Impossibility of sketching of the 3d transportation metric with quadratic cost.\" ICALP, 2016.\n\n[14] N. Lahn and S. Raghvendra. \"An $O(n^{5/4})$ Time $\\varepsilon$-Approximation Algorithm for RMS Matching in a Plane.\" SODA, 2021.\n\n[15] P. K. Agarwal and J. M. Phillips. \"On bipartite matching under the RMS distance.\" CCCG 2006.", "claims_and_evidence": "The paper claims to provide an improved algorithm for approximating the $p$-Wasserstein distance between distributions supported on $n$ atoms up to a multiplicative factor which scales as $O(log(n)$ (in expectation), i.e. the approximation satisfies $0\\leq \\mathbb E[A(\\mu,\\nu)]\\leq C\\log(n) W_p(\\mu,\\nu)$ and that this approximation can be computed in time scaling as $O(n^2\\log n \\log U\\log \\Delta)$. The article carefully explains the ideas underlying this approximation and the mathematical results are supported by proofs provided in the supplement. As such, I have no concerns regarding the accuracy of the claims made in the paper. \n\nThe claims regarding the RPW problem and $\\delta$-additive approximations are also well supported and appear reasonable to me.\n\nThese claims are also empirically validated in numerical experiments which I find to be sufficient. Notably, the multiplicative scaling factor in the approximation is seen to scale sublinearly.", "methods_and_evaluation_criteria": "While the experiments are not extensive, they adequately demonstrate the claimed results.", "theoretical_claims": "I did not verify the proofs in the supplement.", "experimental_designs_or_analyses": "I believe the experiments are sound.", "relation_to_broader_scientific_literature": "The paper advances our understanding of algorithms for numerical resolutions of optimal transport and furnishes an efficient algorithm for estimating the RPW problem (based on a modification of a known algorithm).", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "I believe the paper is quite well-written and provides clear explanations. The algorithmic aspects of the work are particularly well explained. The algorithm for the RPW distance is also a nice contribution and its scaling appears favorable. My concern lies with the quality of the approximation provided for the Wasserstein distance. In effect, a  $O(\\log n)$ multiplicative factor in the quality of the approximation appears very undesirable even if it is more efficient to compute. Although the experiments demonstrate that in two simple examples the dependence on $n$ is not bad, the approximation error is on the order of 1.5, which implies that this value is not a very good estimate for the actual Wasserstein distance.", "other_comments_or_suggestions": "Line 188 left column: distace -> distance", "questions_for_authors": "My only question pertains to the $O(\\log n)$ multiplicative factor in the approximation of the Wasserstein distance. Perhaps I have misunderstood something, but this appears to be be a very serious limitation of the approach and implies that the quality of the approximation is effectively unknown.", "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": "3ywFrtJOIN", "reviewer": "Reviewer_NUZa", "summary": "This paper develops efficient algorithms for approximating the $p$-Wasserstein distance $W_p$ and a robust variant of $W_p$. In particular, when the input measures are uniform over discrete point sets of size $n$, they provide an approximation algorithm for $W_p$ with relative error $O(\\log n)$ that runs in time $O(n^2 \\log \\Delta)$, where $\\Delta$ is the ratio between largest and largest edge cost. Previous results of this form only worked with $p=1$.\n\nFor additive error, they provide a reduction suggesting that getting error $\\delta$ requires time $O(n^2/\\delta^p)$ unless one resorts to a class of currently impractical algorithms. Then, they show that a robust variant of the $W_p$ can be estimated in time $O(n^2/\\delta^3)$ for all $p$.\n\n## Update after Rebuttal\nI maintain my positive evaluation.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough review and constructive feedback. We will update the paper to include your suggestions. We assume that distances can be queried in $O(1)$ time and will state it explicitly. We will also include additional dependencies of execution time on $p$ that arise due to taking the $p$th power of the costs.\n\n>Can the authors elaborate on the $n^{3/2}$ scaling they observe - in particular, do they think there is room for a tighter computational complexity bound?\n\n**Response:**\nYes, it is possible to achieve a bound of $\\tilde{O}(n^{3/2})$ for the special case of the minimum-cost bipartite matching problem provided the HSTs are precomputed. Sharathkumar and Agarwal [1] showed that Gabow and Tarjan's algorithm [2] for min-cost bipartite matching can be implemented using $\\tilde{O}(n^{3/2})$ queries to the bichromatic closest pair data structure. Combining this with the structure of Section 2, we can bound the execution time by $\\tilde{O}(n^{3/2})$ time algorithm. We will include this discussion in the next version of our paper.\n\n[1] R. Sharathkumar and P. K. Agarwal. \"Algorithms for the transportation problem in geometric settings.\" SODA, 2012.\n\n[2] H. Gabow and R. E. Tarjan. \"Faster scaling algorithms for network problems.\" SIAM Journal on Computing 1989.\n\n----\n----\n\n>Since I am not as familiar with their robust variant, I am less sure of the significance of that result. E.g. it could be the case that it is easier to calculate for reasons that make it less useful in practice, though I am not claiming that to be the case.\n\n**Response:** The $p$-Wasserstein distance's sensitivity to noise is a key factor limiting its practical utility and making the design of an additive approximation particularly challenging.\n\nA small noisy mass of $\\delta^p$ can affect the $p$-Wasserstein distance by $\\delta$. Consequently, any $\\delta$-additive approximation algorithm must transport all but $\\delta^p$ mass in an approximately minimum-cost way to ensure a bound of $\\delta$ on the additive approximation, which, under reasonable assumptions takes $\\Omega(n^2/\\delta^p)$ time. In the extreme case, when $p =\\infty$, any approximate solver must transport all the mass in an approximately minimum-cost way which requires $\\Omega(n^{2.5})$ time. In essence, the sensitivity to noise is precisely why additive approximations for the $p$-Wasserstein distance require $\\Omega(n^2/\\delta^p)$ time. \n\nIn contrast, $p$-RPW is more robust, as a noise of $\\delta$ impacts it by at most $\\delta$. Our algorithm essentially transports all but $\\delta$ mass in an approximate minimum-cost way and has an execution time of $O(n^2/\\delta^3)$. Therefore, the $p$-RPW distance is not only robust to noise but also allows for faster algorithms, a direct consequence of this robustness.  \n\n----\n----\n\n>Can the authors please confirm that they will publish their code if the paper is accepted?\n\n**Response:** Yes, we will make our GitHub repository public.\n\n----\n----\n\n\n>When $p=1$, how close is the robust W1 variant you consider to Dudley's bounded Lipschitz distance (IPM w.r.t Lipschitz and bounded functions)?\n\n**Response:** For $p = 1$ and in a metric space with a unit diameter, the $1$-RPW can be written as $\\sup_{Lip(f) \\le 1-||f||_\\infty} |\\int f d\\mu - \\int f d\\nu|$, which is upper-bounded by Dudley's bounded Lipschitz distance.", "claims_and_evidence": "Yes. This is primarily a theoretical paper and all theorems are accompanied by proofs which appear sound. They supported their guarantees with some basic experiments, with results that match or beat those predicted by the theory (which I think are sufficient given theoretical nature of the paper).", "methods_and_evaluation_criteria": "Yes, I think the evaluation settings are fine. I assume they are from the setting where both sample sets are drawn from the same distribution (could be clarified). This makes sense because when the true distance is small the relative error is more tolerable.", "theoretical_claims": "Yes, I read through all of the proofs for the relative approximation result and the reduction, and I read through the main parts of the additive approximation result (though I just skimmed over some auxiliary lemmas in the appendix). Everything looks good to me, assuming that they accurately describe the data structure results from previous work, which I am not very familiar with and did not verify.\n\nThere is a bit of ambiguity on the computational model. For example, reading the input points / computing the cost matrix takes $n^2 d$ time in $d$ dimensional Euclidean space. I think it's fine to primarily assume the distance matrix can be queried in constant time, but there should probably be a remark on this. Also, I am used to seeing additive approximation results which scale with the largest value of the cost matrix. I'm just confirming that there is not such a dependence here (assuming that $\\Delta$ is small and that we can do arithmetic with entries of the distance matrix in constant time).", "experimental_designs_or_analyses": "The experiment setup seem fine as above. The results are in line (or even a bit better) than their theory predicts. By the way, can the authors elaborate on the $n^{3/2}$ scaling they observe - in particular, do they think there is room for a tighter computational complexity bound?\n\nCan the authors please confirm that they will publish their code if the paper is accepted?", "relation_to_broader_scientific_literature": "There is lots of work on statistics and computation of OT and its robust variants, which they discuss. I think they appropriately cite the relevant work on OT computation, though I work more on the statistical side so I cannot be certain. I am only loosely familiar with the data structure results that they use.", "essential_references_not_discussed": "Not that I know of.", "other_strengths_and_weaknesses": "This paper is well-written, and I think that the first result and its proof ideas are quite nice - honestly enough to justify acceptance on their own. I could imagine that their tree metric might have other applications in computational geometry.\n\nSince I am not as familiar with their robust variant, I am less sure of the significance of that result. E.g. it could be the case that it is easier to calculate for reasons that make it less useful in practice, though I am not claiming that to be the case.", "other_comments_or_suggestions": "Could you include the dependence of your complexities on p in the theorem statements? It seems to pretty well accounted for in the proofs, although I suppose there should be another multiplicative overhead of p when you perform arithmetic on entries of the exponentiated cost matrix (depending on the computational model).\n\nDefine LMR initialism\n\nIt seems that each cluster C is defined as a subset of points, so I am not sure what the notation X_C is needed for.\n\nThere seems to be a typo in the definition of $k_j^x$\n\nI think the quantity $d(C)$ used in the statement of Lemma 2.3 is only defined in the proof.", "questions_for_authors": "When p=1, how close is the robust W1 variant you consider to Dudley's bounded Lipschitz distance (IPM wrt Lipschitz and bounded functions)? That distance is also upper bounded by W1 and TV.", "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": "VwwyBD6v79", "reviewer": "Reviewer_GvkL", "summary": "This paper aims to provide an $O(\\log n)$ approximation algorithm for p-Wasserstein distance that runs in $O(n^2 \\log n \\log U \\log \\Delta)$ time. This is done with a collection of $p$ HST trees and with a dynamic BCP data structure to efficiently find augmenting paths in the dual framework for the OT problem. The authors also present an algorithm for an additive approximation of a noise-resistant variant of $p$-Wasserstein called $(p,k)-RPW$. Their approach to this borrows ideas from the work of [Lahn et al. 19], scaling the problem to an integer one and using the [Gabow and Tarjan 89] algorithm. Some simple experiments demonstrate the scaling behavior that is expected, with also evidence that the average-case performance for the first method is between $n$ and $n^{3/2}$.\n\n## Update after rebuttal:\n\nI felt sufficiently satisfied with the explanations for my concerns, so I've updated my score to accept.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough review and constructive feedback. We will add a short discussion and a reference about the convergence rate of the Wasserstein distance, as well as the formal definition of spread and the corrected inequality in Appendix A in our next version. We answer the main concern raised by the reviewer below.\n\n\n>My main reason for not providing a higher score is that the presentation seems a bit disjointed, with two related approximation problems solved with completely different techniques. Much of the intro describes noise sensitivity as an issue, but then the majority of text is spent on the first problem.\n\n\n**Response:** Noise sensitivity of $p$-Wasserstein distance is the primary reason why both relative and additive approximations are challenging to design. We elaborate on this below and edit the paper to include this discussion.\n\n*Relative Approximations:* Embedding a metric into a tree metric introduces noise, or distortion, to the edge costs. Due to the high sensitivity, the optimal $p$-Wasserstein distance with respect to these noisy edge costs fails to serve as a relative approximation to the $p$-Wasserstein distance with respect to the original costs, i.e., the approximation quality is unbounded. To address this, we reduce the distortion in edge costs by selecting the smallest distortion across $p$ different tree embeddings. We demonstrate that the optimal $p$-Wasserstein distance computed with respect to these edge costs is in fact a $O(\\log n)$-approximation of the true $p$-Wasserstein distance.  \n\n\n*Additive Approximations:*  A small noisy mass of $\\delta^p$ can affect the $p$-Wasserstein distance by $\\delta$. Consequently, any $\\delta$-additive approximation algorithm must transport all but $\\delta^p$ mass in an approximately minimum-cost way to ensure a bound of $\\delta$ on the additive approximation, which, under reasonable assumptions takes $\\Omega(n^2/\\delta^p)$ time. In the extreme case, when $p =\\infty$, any approximate solver must transport all the mass in an approximately minimum-cost way which requires $\\Omega(n^{2.5})$ time. In essence, the sensitivity to noise is precisely why additive approximations for the $p$-Wasserstein distance require $\\Omega(n^2/\\delta^p)$ time. \n\nIn contrast, $p$-RPW is more robust, as a noise of $\\delta$ impacts it by at most $\\delta$. Our algorithm essentially transports all but $\\delta$ mass in an approximate minimum-cost way and has an execution time of $O(n^2/\\delta^3)$. Thus, it admits faster algorithms precisely because of its robustness to noise. \n\n----\n----\n\n>In line 207, lev(a,b) is never explicitly defined, I don't think. It is easy enough to infer, but please clarify.\n\n**Response:** The notation is defined in line 186.\n\n\n----\n----\n\n>I don't think I understood the second inequality in the derivation 181-189 (second column, for inequality 1). It was my impression that the telescoping sum would eliminate $H_{k^a_0}$ and $H_{k^b_0}$ if $h \\geq 2$. I believe the result still holds, but would have expected different terms. Please clarify.\n\n\n**Response:** Thanks for raising this point. While the result of the equation in lines 181-189 is correct, the subscript of H in the first line of the equation has to be changed from $j-2$ to $j+2$. Note that in our notation, the root node (the largest cell) is at level 0. The slight typo resulted since the original paper presenting the HST construction [1] had the root at level $h$. We will change $j-2$ to $j+2$ in our next version.\n\n[1] J. Fakcharoenphol, S. Rao, and K. Talwar. \"A tight bound on approximating arbitrary metrics by tree metrics.\" STOC, 2003.", "claims_and_evidence": "This is mostly a theoretical paper, so its evidence is mostly in the proofs. As noted in the section below, I was able to check Section 2 and Appendix A carefully, which included:\n* a O(log n) in expectation approximation of the true p-Wasserstein distance with p HST trees\n* a description of a data structure that allows for efficient construction, retrieval, and insertion/deletion \nThe proofs were correct. I was not able to check Section 3 (which used these tools to calculate their tree approximation distance in near-quadratic time) and Section 4/Appendix B (which proposed their method for additive approximation of $(p,k)$-RPW) in as much detail, but the techniques used were sensible.\n\nThere was also some empirical support for the theoretically established complexities via some simple experiments.", "methods_and_evaluation_criteria": "Yes, the experiments were reasonably done.", "theoretical_claims": "I checked the arguments of Section 2 (including appendix A) in detail and am convinced of their correctness. Some questions on minor details in latter sections. For sections 3 and 4, I did not have time to check them carefully, but the parts in the main text seem fine upon a more cursory inspection.", "experimental_designs_or_analyses": "The experiments were relatively simple, but reasonable, in my opinion.", "relation_to_broader_scientific_literature": "The work does a good job of referencing and discussing prior work. I did not notice any glaring omissions.\n\nOne note though, is that I would have preferred that the authors mention the poor rate of convergence for empirical measures to underlying continuous measures in $p$-Wasserstein distance (a reference is [Weed & Bach 19]: *Sharp asymptotic ... empirical measures in Wasserstein distance*). This is a major practical limitation of the Wasserstein distance that felt conspicuously missing from the sentence in lines 054-055. I understand that it doesn't fit the story well, as I don't think this approach addresses this problem at all, but it felt wrong to omit it entirely.", "essential_references_not_discussed": "See above.", "other_strengths_and_weaknesses": "See below.", "other_comments_or_suggestions": "Rating Explanation: I appreciated the work and the technical depth of the approach to an important problem. My main reason for not providing a higher score is that the presentation seems a bit disjointed, with two related approximation problems solved with completely different techniques. Much of the intro describes noise sensitivity as an issue, but then the majority of text is spent on the first problem which is not looking at the $p-RPW$ distance.\n\nOne clarity nit: At the start of 2.1, I would remind people that one is considering a discrete, finite metric space. I found myself confused at $k^u_j$ being an integer for a moment.", "questions_for_authors": "1. Is the *spread* equal to the ratio of largest to smallest edge costs? These are both denoted with the same $\\Delta$ notation, but are never explicitly connected.\n2. In line 207, lev(a,b) is never explicitly defined, I don't think. It is easy enough to infer, but please clarify.\n3. I don't think I understood the second inequality in the derivation 181-189 (second column, for inequality 1). It was my impression that the telescoping sum would eliminate $H_{k^a_0}$ and $H_{k^b_0}$ if $h \\geq 2$. I believe the result still holds, but would have expected different terms. Please clarify.\n4. In appendix A, for the proof of lemma 2.2, I expected the second equality to be a $\\leq$, as the optimal transport plan for the base metric will be suboptimal for the tree metric, no?", "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": "f1yxgiq9Mk", "reviewer": "Reviewer_ZSVF", "summary": "This work introduces a method to compute a $O(\\log n)$-approximation of the p-Wasserstein distance in $O(n^2 \\log n)$ time for $p\\ge 2$. The method is based on the construction of Hierarchically well Separated Trees (HSTs), and the Hungarian algorithm with Bichromatic Closest Pairs. Furthermore, the authors provide a method to approximate a variant of the p-Wasserstein distance which is robust to noise, named $(p,k)$-RPW, which can be computed in $O(n^2/\\delta^2)$ for an additive error of $\\delta$. This is an improvement compared to other methods such as entropic regularized OT computed with Sinkhorn which are computed in $O(n^2/\\delta^p)$.\n\n\n## Update after rebuttal\n\nI maintain my score as I believe that additional experiments (especially for RPW) and comparisons would strengthen this work.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough review and constructive feedback. The two main contributions of our paper are as follows:\n\n\n* A $O(\\log n)$ relative approximation algorithm for $p=1$ has been known for almost three decades and has had a significant impact. Despite significant effort, no such algorithm was known for $p=2$. We present the first $O(\\log n)$ relative approximation algorithm for the $p$-Wasserstein distance, for any fixed value $p > 1$.\n\n* Under reasonable assumptions, we show that any algorithm that additively approximates $p$-Wasserstein distance requires $\\Omega(n^2/\\delta^p)$ time. This hardness also extends to other robust variants including $\\lambda$-ROBOT and partial $p$-Wasserstein distance. In contrast, the $p$-RPW distance, due to its robustness, admits an $O(n^2/\\delta^3)$ time algorithm for any value $p \\ge 1$. \n\nWe note that our contributions are significant from a theoretical standpoint. The primary reason to implement and test our relative approximation algorithm is to understand the gap between its worst-case theoretical guarantees and its practical performance. Having said that, we can include a comparison of our algorithm with the standard Hungarian algorithm as well as an implementation of our approximation algorithm of $p$-RPW with that of [1].\n\n**Comparison with OT.emd:** There are two reasons for not making a direct comparison between our algorithm and OT.emd.\n\n* OT.emd requires $O(n^2)$ space which means we can only execute them on small instances. This restricts our ability to compare our algorithm with OT.emd and Sinkhorn on larger inputs and understand how their performance scales with input size.\n\n* Our prototype implementation is written in Python where as OT.emd is based on a highly optimized C++ implementation. Comparing the performance of codes written in different languages, especially only in small instances, can mislead us which is why we have avoided a direct comparison. \n\n----\n----\n\n>Maybe I missed it, but I am not sure what \"$an$\" and \"$an^{(3/2)}$\" refer to in Figure 1.b and 1.d.\n\n**Response:** The plot uses a $\\log$–$\\log$ scale, where any polynomial function appears as a straight line, with the slope indicating the polynomial’s exponent. In Plot 1.b, we included the functions $f(n) = an$ and $f(n) = an^{3/2}$ (for some constant $a$) to illustrate that, in our experiments, the running time of the Hungarian algorithm with our data structure is bounded by $O(n^{3/2})$.\n\n----\n----\n\n>Did you implement the algorithm to compute RPW?\n\n**Response:** We show that additive approximation requires $\\Omega(n^2/\\delta^p)$ time due to their sensitivity to noise and the robustness of $p$-RPW allows for significantly faster additive approximation with an execution time of $O(n^2/\\delta^3)$ (for all values of $p$). Our intention was to highlight this gap, which can be accomplished without the need for an implementation. \nNonetheless, we are in the process of implementing our algorithm and commit to comparing it with the algorithm of [1] in the next version of our paper.\n\n----\n----\n\n>How does your algorithm introduced to compute RPW differ from the ones in [1]?\n\n**Response:**\nThe algorithm introduced in [1] computes the OT profile and uses it to approximate the $p$-RPW, which takes $O(n^2/\\delta^p+ n/\\delta^{2p})$ time. We use a guessing procedure along with an early stopping criteria to reduce the execution time to $O(n^2/\\delta^3)$. More precisely, our algorithm picks a guess $g$ from one of the $O(1/\\delta)$ values as an approximation of the true $p$-RPW distance. For this value, it executes only $O(p/\\delta^2)$ iterations of the LMR algorithm [2]. It then returns the minimum estimate achieved across all the $O(1/\\delta)$ guesses. In contrast, the algorithm in [1] runs the LMR algorithm [2] for $O(1/\\delta^p)$ iterations. \n\n\n[1] S. Raghvendra, P. Shirzadian, and K. Zhang. \"A New Robust Partial $p$-Wasserstein-Based Metric for Comparing Distributions.\" ICML 2024.\n\n[2] N. Lahn, D. Mulchandani, and S. Raghvendra. \"A graph theoretic additive approximation of optimal transport.\" NeurIPS 2019.", "claims_and_evidence": "The authors claim that they provide a method to compute a $O(\\log n)$-approximation of the p-Wasserstein distance in $O(n^2 \\log n)$ time for $p\\ge 2$ (Theorem 1.1), and that they can approximate the $(p,k)$-RPW distance in $O(n^2/\\delta^3)$ time for an additive error of $\\delta$ (Theorem 1.2). Both claims are supported by proofs. More precisely, in Section 2, a distance is derived, and is shown to be a $O(\\log n)$ approximation of the Wasserstein distance, and in Section 3, an algorithm with a complexity of $O(n^2 \\log n)$ to compute the proposed distance is described. Finally, the algorithm to approximate the RPW distance is derived in Section 4.", "methods_and_evaluation_criteria": "The methods proposed make sense. The evaluation criteria allow to verify the claims in practice. Nonetheless, I believe it would have been interesting to provide comparisons with known baselines, e.g. comparing the runtime of the method with the implementation of the OT problem in the library Python Optimal Transport. And likewise, it would have been nice to compare the RPW algorithm with e.g. Sinkhorn. Also, I am not sure the results of RPW are presented in Figure 1?", "theoretical_claims": "The theoretical claims seem good. The claims (see Claims and Evidence section) are supported through proofs in Section 2, 3, 4 and in Appendix. However, I did not check all the details.", "experimental_designs_or_analyses": "The experimental design seems good. They demonstrate the complexity of the first algorithm in Figure 1.b and 1.e, and the approximation in Figure 1.f. They also verify the impact of several parameters on the results. However, there are no comparisons with baselines which are usually used to compute the $p$-Wasserstein distance. Moreover, a practical study of the algorithm to approximate RPW seems to be lacking.", "relation_to_broader_scientific_literature": "One of the key contribution is to use a set of $p$ Hierarchically well Separated Tree Metrics instead of one, which was done for $p=1$ in [1, 2]. Another key contribution is to show a method to provide an approximation of the $p$-RPW distance, introduced in [3].\n\n[1] Moses S Charikar. Similarity estimation techniques from rounding algorithms. In Proceedings of the thiry-fourth annual ACM symposium on Theory of computing, pages 380–388, 2002.\n\n[2] Jon Kleinberg and Eva Tardos. Approximation algorithms for classification problems with pairwise relationships: Metric labeling and markov random fields. Journal of the ACM (JACM), 49(5):616–639, 2002.\n\n[3] Sharath Raghvendra, Pouyan Shirzadian, and Kaiyi Zhang. A new robust partial p-Wasserstein-based metric for comparing distributions. In 41st International Conference on Machine Learning, 2024.", "essential_references_not_discussed": "All the essential references seem to be discussed.", "other_strengths_and_weaknesses": "**Strengths**\n- Interesting method to compute an approximation of the $p$-Wasserstein distance in $O(n^2 \\log n)$ time.\n- Also provide an interesting method to compute an approximation of the $p$-RPW distance.\n- Verify the theoretical results in practice.\n\n\n**Weaknesses**\n- The procedures proposed can be a bit hard to understand, as they are described in the text. Algorithms/Figures could maybe help the readers to better understand the methods.\n- The experimental verifications are done only with the procedure to approximate the $p$-Wasserstein distance (if I understand well). Moreover, there are no comparisons with baselines which are usually used to compute the $p$-Wasserstein distance.", "other_comments_or_suggestions": "I suggest improving readibility by adding algorithms or figures to help the reader understand how to implement the algorithm.\n\nIt would also be nice to have experimental results for the computation of RPW, and comparisons with benchmark methods. For instance, the algorithm should be faster than the hungarian algorithm, but it does not appear on the experimental sections.\n\nMaybe I missed it, but I am not sure what \"$an$\" and \"$an^{(3/2)}$\" refer to in Figure 1.b and 1.d.\n\nTypos:\n- Line 189, 1st column: \"distace\"", "questions_for_authors": "1. Is your approximation algorithm faster than ot.emd in the library POT?\n2. Did you implement the algorithm to compute RPW?\n3. How does your algorithm introduced to compute RPW differ from the ones in [1]?\n\n[1] Sharath Raghvendra, Pouyan Shirzadian, and Kaiyi Zhang. A new robust partial p-Wasserstein-based metric for comparing distributions. In 41st International Conference on Machine Learning, 2024.", "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:

Meta-review: The authors propose an $\mathcal{O}(\log{n})$-approximation algorithm for $p$-Wasserstein ($p \ge 2$) for discrete measures in time complexity $\mathcal{O}(n^2 \log{n})$, where $n$ is the number of supports of input measures. Intuitively, the proposed method is based on the construction of Hierarchically Well-Separated Trees (Fakcharoenphol et al., 2003) and the Hungarian algorithm with Dynamic Bichromatic Closest Pairs data structure. Additionally, the authors propose to approximate a robust variant of $p$-Wasserstein (i.e., $p$-PRW) within an additive error $\delta$ in time complexity $\mathcal{O}(n^2 / \delta^3)$, in contrast to $\Omega(n^2 / \delta^p)$ for existing combinatorial approaches in the literature. All Reviewers agree that it is a valuable contribution. The Reviewers raised concerns on empirical comparison to other popular approaches in the literature to illustrate the practical benefits of the proposed approach, i.e., whether one can realize an efficient implementation of the proposed approach to match with its theoretical complexity. We suggest the authors to incorporate the Reviewers' comments and discussions in the rebuttal into the updated version, which helps to strengthen the work.
Decision: Accept (poster)