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.

Optimal Transport is a popular distance metric for measuring similarity between distributions. Exact and approximate combinatorial algorithms for computing the optimal transport distance are hard to parallelize. This has motivated the development of numerical solvers (e.g. Sinkhorn method) that can exploit GPU parallelism and produce approximate solutions. We introduce the first parallel combinatorial algorithm to find an additive $\varepsilon$-approximation of the OT distance. The parallel complexity of our algorithm is $O(\log(n)/ \varepsilon^2)$ where $n$ is the total support size for the input distributions. In Massive Parallel Computation (MPC) frameworks such as Hadoop and MapReduce, our algorithm computes an $\varepsilon$-approximate transport plan in $O(\log (\log (n/\varepsilon))/\varepsilon^2)$ rounds with $O(n/\varepsilon)$ space per machine; all prior algorithms in the MPC framework take $\Omega(\log n)$ rounds. We also provide a GPU-friendly matrix-based interpretation of our algorithm where each step of the algorithm is row or column manipulation of the matrix. Experiments suggest that our combinatorial algorithm is faster than the state-of-the-art approximate solvers in the GPU, especially for higher values of $n$.

Maximum cardinality bipartite matching is an important graph optimization problem with several applications. For instance, maximum cardinality matching in a $\delta$-disc graph can be used in the computation of the bottleneck matching as well as the $\infty$-Wasserstein and the Lévy-Prokhorov distances between probability distributions. For any point sets $A, B \subset \mathbb{R}^2$, the $\delta$-disc graph is a bipartite graph formed by connecting every pair of points $(a,b) \in A\times B$ by an edge if the Euclidean distance between them is at most $\delta$. Using the classical Hopcroft-Karp algorithm, a maximum-cardinality matching on any $\delta$-disc graph can be found in $\tilde{O}(n^{3/2})$ time.~\footnote{We use $\tilde{O}(\cdot)$ to suppress poly-logarithmic terms in the complexity.} In this paper, we present a simplification of a recent algorithm (Lahn and Raghvendra, JoCG 2021) for the maximum cardinality matching problem and describe how a maximum cardinality matching in a $\delta$-disc graph can be computed asymptotically faster than $O(n^{3/2})$ time for any moderately dense point set. As applications, we show that if $A$ and $B$ are point sets drawn uniformly at random from a unit square, an exact bottleneck matching can be computed in $\tilde{O}(n^{4/3})$ time. On the other hand, experiments suggest that the Hopcroft-Karp algorithm seems to take roughly $\Theta (n^{3/2})$ time for this case. This translates to substantial improvements in execution time for larger inputs.

The 2-Wasserstein distance (or RMS distance) is a useful measure of similarity between probability distributions with exciting applications in machine learning. For discrete distributions, the problem of computing this distance can be expressed in terms of finding a minimum-cost perfect matching on a complete bipartite graph given by two multisets of points A, B ⊂ ℝ2, with |A| = |B| = n, where the ground distance between any two points is the squared Euclidean distance between them. Although there is a near-linear time relative ∊-approximation algorithm for the case where the ground distance is Euclidean (Sharathkumar and Agarwal, JACM 2020), all existing relative ∊-approximation algorithms for the RMS distance take Ω(n3/2) time. This is primarily because, unlike Euclidean distance, squared Euclidean distance is not a metric. In this paper, for the RMS distance, we present a new ∊-approximation algorithm that runs in O(n5/4 poly{log n, 1/∊}) time. Our algorithm is inspired by a recent approach for finding a minimum-cost perfect matching in bipartite planar graphs (Asathulla et al, TALG 2020). Their algorithm depends heavily on the existence of sublinear sized vertex separators as well as shortest path data structures that require planarity. Surprisingly, we are able to design a similar algorithm for a complete geometric graph that is far from planar and does not have any vertex separators. Central components of our algorithm include a quadtree-based distance that approximates the squared Euclidean distance and a data structure that supports both Hungarian search and augmentation in sublinear time.

Given a weighted planar bipartite graph G ( A ∪ B , E ) where each edge has an integer edge cost, we give an Õ( n 4/3 log nC ) time algorithm to compute minimum-cost perfect matching; here C is the maximum edge cost in the graph. The previous best-known planarity exploiting algorithm has a running time of O ( n 3/2 log n ) and is achieved by using planar separators (Lipton and Tarjan ’80). Our algorithm is based on the bit-scaling paradigm (Gabow and Tarjan ’89). For each scale, our algorithm first executes O ( n 1/3 ) iterations of Gabow and Tarjan’s algorithm in O ( n 4/3 ) time leaving only O ( n 2/3 ) vertices unmatched. Next, it constructs a compressed residual graph H with O ( n 2/3 ) vertices and O ( n ) edges. This is achieved by using an r -division of the planar graph G with r = n 2/3 . For each partition of the r -division, there is an edge between two vertices of H if and only if they are connected by a directed path inside the partition. Using existing efficient shortest-path data structures, the remaining O ( n 2/3 ) vertices are matched by iteratively computing a minimum-cost augmenting path, each taking Õ( n 2/3 ) time. Augmentation changes the residual graph, so the algorithm updates the compressed representation for each partition affected by the change in Õ( n 2/3 ) time. We bound the total number of affected partitions over all the augmenting paths by O ( n 2/3 log n ). Therefore, the total time taken by the algorithm is Õ( n 4/3 ).

We present a weighted approach to compute a maximum cardinality matching in an arbitrary bipartite graph. Our main result is a new algorithm that takes as input a weighted bipartite graph $G(A\cup B,E)$ with edge weights of $0$ or $1$. Let $w \leq n$ be an upper bound on the weight of any matching in $G$. Consider the subgraph induced by all the edges of $G$ with a weight $0$. Suppose every connected component in this subgraph has $\mathcal{O}(r)$ vertices and $\mathcal{O}(mr/n)$ edges. We present an algorithm to compute a maximum cardinality matching in $G$ in $\tilde{\mathcal{O}}( m(\sqrt{w}+ \sqrt{r}+\frac{wr}{n}))$ time. When all the edge weights are $1$ (symmetrically when all weights are $0$), our algorithm will be identical to the well-known Hopcroft-Karp (HK) algorithm, which runs in $\mathcal{O}(m\sqrt{n})$ time. However, if we can carefully assign weights of $0$ and $1$ on its edges such that both $w$ and $r$ are sub-linear in $n$ and $wr=\mathcal{O}(n^γ)$ for $γ< 3/2$, then we can compute maximum cardinality matching in $G$ in $o(m\sqrt{n})$ time. Using our algorithm, we obtain a new $\tilde{\mathcal{O}}(n^{4/3}/\varepsilon^4)$ time algorithm to compute an $\varepsilon$-approximate bottleneck matching of $A,B\subset\mathbb{R}^2$ and an $\frac{1}{\varepsilon^{\mathcal{O}(d)}}n^{1+\frac{d-1}{2d-1}}\mathrm{poly}\log n$ time algorithm for computing $\varepsilon$-approximate bottleneck matching in $d$-dimensions. All previous algorithms take $Ω(n^{3/2})$ time. Given any graph $G(A \cup B,E)$ that has an easily computable balanced vertex separator for every subgraph $G'(V',E')$ of size $|V'|^δ$, for $δ\in [1/2,1)$, we can apply our algorithm to compute a maximum matching in $\tilde{\mathcal{O}}(mn^{\fracδ{1+δ}})$ time improving upon the $\mathcal{O}(m\sqrt{n})$ time taken by the HK-Algorithm.

Transportation cost is an attractive similarity measure between probability distributions due to its many useful theoretical properties. However, solving optimal transport exactly can be prohibitively expensive. Therefore, there has been significant effort towards the design of scalable approximation algorithms. Previous combinatorial results [Sharathkumar, Agarwal STOC '12, Agarwal, Sharathkumar STOC '14] have focused primarily on the design of near-linear time multiplicative approximation algorithms. There has also been an effort to design approximate solutions with additive errors [Cuturi NIPS '13, Altschuler \etal\ NIPS '17, Dvurechensky \etal\, ICML '18, Quanrud, SOSA '19] within a time bound that is linear in the size of the cost matrix and polynomial in $C/\delta$; here $C$ is the largest value in the cost matrix and $\delta$ is the additive error. We present an adaptation of the classical graph algorithm of Gabow and Tarjan and provide a novel analysis of this algorithm that bounds its execution time by $\BigO(\frac{n^2 C}{\delta}+ \frac{nC^2}{\delta^2})$. Our algorithm is extremely simple and executes, for an arbitrarily small constant $\eps$, only $\lfloor \frac{2C}{(1-\eps)\delta}\rfloor + 1$ iterations, where each iteration consists only of a Dijkstra-type search followed by a depth-first search. We also provide empirical results that suggest our algorithm is competitive with respect to a sequential implementation of the Sinkhorn algorithm in execution time. Moreover, our algorithm quickly computes a solution for very small values of $\delta$ whereas Sinkhorn algorithm slows down due to numerical instability.

We give an Õ(n7/5 log(nC)) time1 algorithm to compute a minimum-cost maximum cardinality matching (optimal matching) in Kh-minor free graphs with h = O(1) and integer edge weights having magnitude at most C. This improves upon the Õ(n10/7 log C) algorithm of Cohen et al‥ [SODA 2017] and the O(n3/2 log(nC)) algorithm of Gabow and Tarjan [SIAM J. Comput. 1989]. For a graph with m edges and n vertices, the well-known Hungarian Algorithm computes a shortest augmenting path in each phase in O(m) time, yielding an optimal matching in O(mn) time. The Gabow-Tarjan [SIAM J. Comput. 1989] algorithm computes, in each phase, a maximal set of vertex-disjoint shortest augmenting paths (for appropriately defined costs) in O(m) time. This reduces the number of phases from n to and the total execution time to . To obtain our speed-up, we relax the conditions on the augmenting paths and iteratively compute, in each phase, a set of carefully selected augmenting paths that are not restricted to be shortest or vertex-disjoint. As a result, our algorithm computes substantially more augmenting paths in each phase, reducing the number of phases from to O(n2/5). By using small vertex separators, the execution of each phase takes Õ(m) time on average. For planar graphs, we combine our algorithm with efficient shortest path data structures to obtain a minimum-cost perfect matching in Õ(n6/5 log (nC)) time. This improves upon the recent Õ(n4/3 log (nC)) time algorithm by Asathulla et al. [SODA 2018].

Given a weighted planar bipartite graph G(A ∪ B, E) where each edge has a positive integer edge cost, we give an Õ(n4/3 log nC) time algorithm to compute minimum-cost perfect matching; here C is the maximum edge cost in the graph. The previous best known planarity exploiting algorithm has a running time of O(n3/2 log n) and is achieved by using planar separators (Lipton and Tarjan ’80). Our algorithm is based on the bit-scaling paradigm (Gabow and Tarjan ’89). For each scale, our algorithm first executes O(n1/3) iterations of Gabow and Tarjan's algorithm in O(n4/3) time leaving only O(n2/3) vertices unmatched. Next, it constructs a compressed residual graph H with O(n2/3) vertices and O(n) edges. This is achieved by using an r-division of the planar graph G with r = n2/3. For each partition of the r-division, there is an edge between two vertices of H if and only if they are connected by a directed path inside the partition. Using existing efficient shortest-path data structures, the remaining O(n2/3) vertices are matched by iteratively computing a minimum-cost augmenting path each taking Õ(n2/3) time. Augmentation changes the residual graph, so the algorithm updates the compressed representation for each affected partition in O(n2/3) time. We bound the total number of affected partitions over all the augmenting paths by O(n2/3 log n). Therefore, the total time taken by the algorithm is Õ(n4/3).