A Robust Exact Algorithm for the Euclidean Bipartite Matching Problem
Authors: Akshaykumar G Gattani
Sharath Raghvendra
Pouyan Shirzadian
Venue: NeurIPS 2023 poster (2023)
Keywords: Euclidean bipartite matching
exact algorithms
primal dual method

Abstract: Algorithms for the minimum-cost bipartite matching can be used to estimate Wasserstein distance between two distributions.
Given two sets $A$ and $B$ of $n$ points in a $2$-dimensional Euclidean space, one can use a fast implementation of the Hungarian method to compute a minimum-cost bipartite matching of $A$ and $B$ in $\tilde{O}(n^2)$ time. Let $\Delta$ be the spread, i.e., the ratio of the distance of the farthest to the closest pair of points in $A\cup B$. In this paper, we present a new algorithm to compute a minimum-cost bipartite matching of $A$ and $B$ with a similar worst-case execution time of $\tilde{O}(n^2 \log \Delta)$. However, when $A$ and $B$ are drawn independently and identically from a fixed distribution that is not known to the algorithm, the execution time of our algorithm is, in expectation, $\tilde{O}(n^{7/4}\log \Delta)$.

To the best of our knowledge, our algorithm is the first one to achieve a sub-quadratic execution time even for stochastic point sets with real-valued coordinates.
Our algorithm extends to any dimension $d$, where it runs in $\tilde{O}(n^{2-\frac{1}{2d}}\Phi(n))$ time for stochastic point sets $A$ and $B$; here $\Phi(n)$ is the query/update time of a dynamic weighted nearest neighbor data structure. 
Our algorithm can be seen as a careful adaptation of the Hungarian method in the geometric divide-and-conquer framework.

A Robust Exact Algorithm for the Euclidean Bipartite Matching Problem
Akshaykumar G. Gattani', Sharath Raghvendra', and Pouyan Shirzadian'
Department of Computer Science, Virginia Tech
Abstract
Algorithms for the minimum-coe bipartite matching can be used to estimate Wasserstein distance between two distributions. Given two sets and B n points in a 2-dimensional Euclidean space, one can use fast implementation of the Hungarian method to compute a minimum-cost bipartite matching A and in Ō(n²) time. Let be the spread, i.c., the ratio the distance of the farthest to the closest pair of points in A U B. In this paper, we present new algorithm to compute minimum-com bipartite matching of A and B with a similar worst-case execution time of Ón² log A). However, when A and are drawn independently and identically from fixed distribution that is not known to the algorithm, the execution time of our algorithm is, in expectation, (n/7/ log A).
To the best of knowledge, quadratic execution time even for stochastic point sets with realevalued coordinates. Our algorithm extends to any dimension d, where runs in Ō z ((n)) time for stochastic point sets A and B; here ((n) is the query/update time of a dy- namic weighted nearest neighbor data structure. Our algorithm can be seen as a careful adaptation of the Hungarian method in the geometric divide-and-conquer framework.
algorithm the first
achieve sub-
Introduction
Given two distributions jr and v defined on the sets Sp, Sy S Rd, respectively, and an integer p 2 1, the p- p-aassrsteein distance between l and is the minimum cost required to transport mass from one to the other. More formally, let [(, v) be the set of all probability measures on Sn x S, with marginal distributions f and v. The p- W-ssserstein distance between p and v is defined as
1/p
Wp(nvv)
inf (ET(H.P) Sp XS.
||sI S2 5211ddy(s1, (2)
where ||sı s2ll denotes the Euclidean distance of $1 and $2.
In many applications, the distribution y (resp. v)) may not be known, but one may have access to n i.idd samples A (resp. B) from y (resp. v). The empirical distribution Pa (resp. vn,) is discrete distribution with A (resp. B) as support, where each sample has a mass of 1/n. The empirical p- Wassersttii distance is simply the p- p-Wasserstein distance between the empirical distributions Nn and P.. As n + co, the empirical >-assersteein distance between and v converges to the real p- passersstnin distance between them 77 [5] 16 46 with several results showing sharp upper
and lower bounds on the empirical p- Wasserstein distances
[10 17 19
45
491
See also [36]
"Following convention from Theoretical Computer Science, all authors ordered alphahetically.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

for a survey on such results. Due to these properties, empirical p- Wasserstein distance has found
applications in training generative adversarial networks 8 9
24
33 B5], image retrieval [3
42
graph predictions 27 37.. clustering stability validation [8 31 and two-sample tests [3 14 23
25
41].
The convergence of the empirical distribution to the real distribution, however, exhibits the
"curse of dimensionality": that is, the convergence rate decreases as the dimension increases [6] 49
The empirical p- Wasserstein distance therefore is most useful in low-dimensional settings.
Computing the empirical p- p-Wasserstein distance can be done by solving an instance of the minimum- cost bipartite matching problem as follows. Consider the complete bipartite graph on A and B where each edge from E A to b E B has cost of aa bll. A bipartite matching (or simply a matching) M on A B is a set of vertex- vertex-disjoint edges of (A, B). The matching M is said to be perfect matching if [ = n. We define the cost of a matching M, denoted by wp (M)). as "P((M (aa ||a blll_ The pth power Euclidean bipartite matching is a perfect matching M; with the minimum-com cost under (,).). The empirical p- Wasserstein distance between and v is Wp(m:: vn)) = (* up (M,))) 1/p
The minimum-cost matching M; and therefore, the empirical p- Wassersteiin distance be computed in (nr)) time using the Hungarian algorithm [28]. In this paper, we adapt the Hungarian algorithm in a geometric wwiid-anddconque framework. The execution time of our algorithm is similar to that of the Hungarian algorithm (Remark 3.2). However, when the input point sets A and B are samples drawn from the same distribution, the algorithm is asymptotically faster than the Hungarian algorithm.
Related Work. Computing an optimal minimum-cost bipartite matching can be done in On time using the classical Hungarian algorithm (28). In geometric settings, the efficiency of the Hungarian algorithm can be improved to (nn² where ((n) is the query and update time of dynamic weighted nearest neighbor data structure therefore, we get a (nn²) time exact algorithm. For higher dimensions, however, this leads to only slightly sub-cubic execution time.
48 For two dimensions, & (n) log O(1), n and
For graphs with n vertices and m edges, there are faster algorithms that compute an optimal minimum- provided the edge costs are integers that are bounded by c However, in our setting, the coordinates of the input points are real-valued and the Euclidean distances, due the presence of a square-root, can be irrational. As a result, these algorithms only provide an approximate solution and the Hungarian algorithm and its fast implementations remain the only known exact algorithm for the problem.
cost matching [II 21 22
For 2-dimensional points with integer coordinates bounded by I, Sharathkumar 4 presented weakly polynomial Onn3/2 log I') time algorithm to compute the exact 1-Wasserstein distance. The main idea was to use the approximation algorithm by Sharathkumar and Agarwal 44 to initially find a planar sub-graph that traps all the edges of the optimal matching. Then, they use the planar separator based algorithm by Lipton and Tarjan (32) to find a minimum-cost matching in this sub- graph. This algorithm, however, restricted only to two dimensions and for p 1. The design of a faster algorithm to compute the pth power Euclidean bipartite matching for d-dimensional point sets remains major open problem. We would also note the extensive work on the design of
approximation algorithms for the I-Wasserstein distance (2 4 the 2- Wasserstein distance [
20 26
30
38
40
as well as
Our Results. For any point set P in the Euclidean space, the spread of P is the ratio of the distance of the farthest pair to the closest pair in P. The main result of this paper is presented in Theorem I.I Theorem 1.1. There exists randomized algorithm that, given any two point sets A and B sampled independently and identically from a distribution j inside the unit d-dimensional hypercube, where |A| = [B| = n and p is not known to the algorithm, computes an exact Euclidean bipartite matching between A and B and has expected running time of Ó(n² 2a (()) log A); here, A is the spread of the points in AUA. B.
When A and B are i.i.d samples from a distribution 1s in the unit square (d = 2) and l is not known to the algorithm, our algorithm achieves weakly polynomial sub-quadratic execution time of 7/ log A). Our algorithm easily extends to any constant p > 1 and any fixed d-dimensional
'The execution time of these algorithms have dependence log C. making them weakly polynomial.

---

v(h
vla)
y(b)
y(b)
y
y(a
(a)
(b)
(c)
Figure 1: (a) An infeasible edge, (b) a feasible matching edge satisfying constraint (2) (c) a feasible
non-matching edge satisfying constraint (T.
space (we present our result for any arbitrary well as the proof of all our claims in the Appendix).
For instance, for p and 2 as well for p = and d 3, the execution time of our algorithm can be bounded by O( (n) log A). However, for simplicity in exposition of the ideas, we restrict the presentation of our algorithm to two-dimensional (d = 2) and Fuclidean settings (n = 1). To the best of our knowledge, this is the first sub-quadratic weakly polynomial exact algorithm for computing the pth power Euclidean bipartite matching for 2-dimensional stochastic point sets. Furthermore, for many distributions, such as the uniform distribution on the unit square, one can show that, with high probability, the spread of the point set is bounded by a polynomial in For all such distributions, we achieve a strongly polynomial sub-quadratic exact algorithm.
Our algorithm can be seen an elegant adaptation of the classical Hungarian algorithm in quadtree based divide-and-conquer framework. The previous speed-ups of the Hungarian algorithm in geometric settings are based on sophisticated data structures that are hard to implement. In contrast, our geometric improvements are relatively simple, which allows us to implement and compare our algorithm with the standard implementation of the Hungarian algorithm. Experiments suggest that our algorithm outperforms the Hungarian algorithm for synthetic data samples drawn from two (possibly different) distributions and long as the I-Wasserstein distance between μ and is small. Our algorithm also outperforms the Hungarian algorithm on data samples drawn from the New York Taxi dataset, where one sample is drawn from the set of request pick-up locations (μ) and the other sample is drawn from the request drop-off locations (v). Experiments also suggest that our analysis may not be tight, at least for the where μ is the uniform distribution inside the unit square. Obtaining a tighter analysis of our algorithm remains an important open question. Next, we present overview of our hwvide-and-conuee algorithm.
Our Approach. The classical Hungarian algorithm is based on the primal-dus dual framework where, along with a matching M, for every point v E A U B, the algorithm maintains non-negative dual weight y(v). The matching M and dual weights y(-) are feasible if, for every edge (a, 6) E A x B,
ybb) y(a) s || bl].
(a, b) $ M, (a, M.
(1)
((6) v(a) ||a bll.
(2)
The Hungarian algorithm initializes M to be empty matching and the dual weights of every vertex v E AUB to be O. It then incrementally builds a perfect matching M while maintaining its feasibility and returns the constructed feasible perfect matching. is well-known that a perfect matching that is feasible is also a minimum-cost matching.
In the Euclidean setting, we interpret the dual weight y(v) of any point the radius of disc centered around v. We refer to these discs as the dual discs. For any edge (a, 6). if the dual disc of a is in the interior of the dual disc of (Figure a)). then y(b) y(a) > l|e 6l1. making the
dual assignment infeasible. Therefore, for any feasible dual weight assignment, the dual disc
cannot be completely inside the dual disc of b. The condition 2 corresponds to Figure Ibb where
the boundary of the disc of a touches the boundary of the disc of b from inside. Figure [c) is an example of an edge that only satisfies II]
In order to implement the Hungarian algorithm in the geometric divide-and-conque framework, we can use quadtree Q with the unit square as the root node. Each node of this tree represents square O and if has more than one input point, then it has four children that are obtained by splitting D into four equal squares. Squares with no more than one point become leaf node in this tree.

---

b
a
b
(
(b)
Figure 2: (a) Feasible matchings for children that is infeasible when combined at the parent (the edge (b, a') is infeasible). (b) Restricted feasible matchings inside the children, which is also restricted feasible when combined at parent.
Consider the following naive implementation of the Hungarian algorithm in quadtree based geomet- ric divide-and-conquee framework. Given square D of the quadtree Q, recursively find a feasible matching in each of the four children. Doing so, however, may lead to a dual weight assignment that violates the feasibility condition for bridge edges, i.c., edges that cross the boundary of the children squares. See Figure a) that illustrates a dual feasible assignment at two children squares, which when combined causes the edge (a', b) to violate the feasibility condition (I. Thus, implementing the conquer step becomes challenging since we cannot simply combine the feasible dual assignments the children to achieve a feasible dual assignment the parent.
To overcome this difficulty, we add an additional condition that restricts the dual disc of a point inside any child square C' to always stay inside D' i.e., y(v) I d'v,,'), where d(v, is the Euclidean distance of u to the boundary of ['. Figure 2 shows a restricted dual assignment at the two children squares, which remains feasible when combined at the parent. Such restricted feasible matching inside each the four children squares of II trivially combine to form a feasible matching framework to compute such restricted feasible
O. Thus, can apply the divide-and-come matching. Unfortunately, however, the restricted feasible matchings computed at the children may be of extremely low cardinality. This is particularly true when most of the input points are close to the boundaries of the children squares. See, for instance, in Figure Baa, the restricted feasible matching an empty one since most of the blue points are close the boundary. As result, the divide step may not make any progress in computing a perfect matching and the conquer step may have do significant amount of work to find perfect matching, causing the execution time of the conquer step to be the same as that of the standard implementation of the Hungarian algorithm.
We observe that if the quadtree Q is constructed using random-shift, then, in expectation, only sub-linearly many points are close to the boundary of any square of Q. By combining this with the fact that the average length of matching edge for n points drawn from an arbitrary distribution μ inside unit square is x 1/vn,, we are able to bound the expected number of unmatched points in the restricted matchings by 3(43/),, leading to an execution time of 7/74)) for the conquer steps across all cells at each level of the quadtree and an overall execution time of 7n274 /4 log A).
2 Preliminaries
We begin by introducing the notations necessary to describe our algorithm. Given A U B and any square D. let Ap A A O, Bn B n D. and ni An U B=|. Let tn denote the side-length of [L For any point v E An U Bg inside D, let d D) denote the Fuclidean distance of o to its closest point on the boundary of I. We say that v is 8-close to O if d(v D) I 86g. Let n denote the number of points of An U B3 that are 8-close to D.
Randomly Shifted Quadtree. Given the input points A U B inside a unit square, a randomly shifted quadtree on the input can be constructed as follows. Let ę be a point chosen uniformly at random from the unit square [o. 1|2. Define the square 0 : [-4, 4 + È to be the root of the quadtree Q. Recursively construct Q by decomposing any square O with na > 1 into four equal squares, each of which become the children of D. Any square with exactly one point becomes a leaf square of the quadtree. Given that the spread of the points in A U B is A, the height of Q is O(log A).

---

(a)
(b)
(c)
Figure 3: (a) D-MCM where points of B (blue points) are unmatched, (b) the boundaries of the cells of the randomly shifted quadtree (solid line) and the 8-close points to the cells of the quadtree (points in the shaded area), (c) the matchings constrained to the cells of the quadtree.
For any square of the quadtree Q, in the next lemma, we show that due to the random shift, there are not many points that are very close to the boundary of C.
Lemma 2.1. For any square D of a randomly shifted quadtree and any E (0,1/2), E [m O [nn])).
Constrained Matching, For any square O. we introduce variant of the minimum-cost matching problem, which refer by the minimum-cost D-comstrained matching problem. Given square D. consider a matching Mn of the points in An U Bn. Any point E An U Bn is afree point with respect to Mn if v is not matched in Mn. Let A5 (resp. B denote the set of free points of An (resp. Bn) with respect to Mn. We define the D-contrrained cost of the matching Mn, denoted by (M()) as
(n) = Σ aa-
Σ d(b,
(3)
Ma
BEBE
For any square C. the minimum-cost D-constrained matching (D-MCM), denoted by Mn, is simply a matching with the minimum -cconstrained cost. In Lemma 2.2 we show that for the root square O of Q, the minimum-cost D -constrained matching is also minimum-cost matching on A AU B.
For any arbitrary square C of Q, MC might be of very small cardinality. See, for instance in Figure aa) where B is the set of blue points, Equation (3) is minimized when all blue points are free points; I.e., (M
E B d(b, [) and MD is an empty matching. Nonetheless, by using the bound on the number of 8-close points to any square D of randomly- quadtree, we can bound the expected number of free points with respect any --MCM MC. the example of FigureBJb) the solid lines show the boundary of the squares of a randomly-shife quadtree, and as shown in FigureB only a few points will be free with respect the optimal matchings constrained the squares of the quadtree.
Lemma 2.2. For any square O of a randomly shifted quadtree and any minimum-cost E-constrained matching MD (i) the expected number of free points of Bu with respect to Må Ô (n²/4) in 2 2-dimensions (and Ó(n--2i 2a) ) in d dimensions), and (ii) if is the root square, then Ma is minimum-cost perfect matching on AU B.
For the rest of the paper, we design a primal-dual method to compute MMMMCM, which by Lemma 2.2 also minimum cost matching on U B.
Constrained Feasibility. Similar to the Hungarian algorithm, we devise a primal-dual method to compute a D-MCM. For any square O of Q. say that matching Mu on Ag B3 along with set of non-negative dual weights y(-) for the points in A U Bu is D-feasible if,
v(b) v(a) < ||a bl1.
V(a, b) A B, V(a,b) Mo.
(4)
(6) - y(a) = ||a bl].
(5)
(66 < d(6, O),
vb E B3.
(6)
s(a) 0,
Va E A5.
(7)
5

---

b2
b
b2
b3
a
(b)
(c)
Figure 4: (a) Cl-eesible matching Mg where b1 is a f-rre point and by is a free point that is not D-free (b) an admissible path P (dashed line) from bi to point b3 whose vertex slack is zero, and (c) the O-optimal matching obtained from augmenting MD along P.
Note that (4) and are conditions identical to the ones maintained by the Hungarian algorithm. Condition ensures that the dual discs stay inside the square and condition (7) ensures that the C-MCM is of minimum cost.
Given a D-feasible matching Mo,v(-), we say that any point b Bu is D-free point with respect to Mn if & is not matched in Mc and y(b) < d(b, D). Note that a free point may not be D-free. For instance, in Figure 4 a). the blue point 62 is free but not a L-free point since its dual weight is equal to its distance to L. Let B denote the set of all D-free points in Bo. We say that --fcasible matching Mn. y-)) is D-optimal if it does not have any D-free points, i.e., B = 0. The following lemma shows that any O-optimal matching Mn. y-) D-MCM. Thus, we focus on computing a C-optimal matching.
Lemma 2.3. Let Ma-y( (-) be D-optimal matching on Ag U Bg. Then, Mg is a minimun-cost D-constrained matching.
Given square D of Q and its four children 11 D22 Dg and 44. we can obtain a L-feasible matching O by simply combining the D.-optimal matchings from each child C. We use this property to design our divide-an and- and-conquer algorithm for computing D-optimal matching.
O Lemma 2.4. For any square C. let D, [1,4] be the of all children of O and let Mi, v()) denote --poiimall matching. Then, the matching Ui_1 M,, v(-) is I-feasible matching.
NexI, we define slack and an augmenting path with respect to a D-feasible matching.
Slacks. For any square D. any f-feasible matching Mu. y(),, and any pair of points (a, b) E An x Br, we define the edge slack on the feasibility conditions for edge (a. b). denoted by s(a. b). as s(a, b) = ||a bll - y(b) + y(a). We say that (a. b) is admissible if s(a, b) = 0. Furthermore, for any point b E Bo, we define the vertex slack of with respect to condition (6) to be (b) ==
d(b, D) - y(b). We refer to the edge slack (resp. vertex slack) as slack when the edge (resp. vertex) obvious from the context.
While the definition of edge slacks are identical to what is used in the Hungarian algorithm, the definition of vertex slacks is new. Since the dual weights associated with any L-feasible matching both edge slacks and vertex slacks will be non-negative.
satisfies
Admissible Augmenting Path. For any square D. suppose Mn, v(-) is D-feasible matching on An U Bn. An alternating path with respect to Mn is a simple path on An U Bn whose edges alternate between those in M= and those not in Mn. An admissible augmenting path (or simply admissible path) is any alternating path consisting only of zero slack edges that starts at a --free point b E BE and ends either (i) free point Añ (ii) point D/ E BD with slack of 0. We can augment a D-feasible matching Mg along an admissible path P by setting Mp + Mp @ P. i.e., we remove all matching edges of from Mp and add all non-matching edges of P to Mo (see Figure The following lemma shows that augmenting a matching along an admissible path does not violate C-feasibility conditions and also reduces the number of O-free points of Bu by one.

---

Lemma 2.3. Suppose Mn, v()) is a C-feasible matching and P is admissible path. After augmenting Mn along P, the matching Mn. y(.) remains D-feasible. Furthermore, the augmentation reduces the number of --rree points of Bu with respect to Ma by one.
Residual Network. In order to assist in computing admissible paths with respect to a D-feasible matching Mo. v(),, define residual network of An U B3 with respect to Mn follows. The vertex set of the residual network is the set of points Ag U Bn and a source vertex 8. For any pair of points (a, b) Ao B3. if (a, b) E MD, add an edge directed from to b with weight s(a, b) to the residual network. Otherwise, (a, b) & Mn. we add an edge directed from to a with weight s(a b). In addition, we add zero-weight edges from the source BB Note that any zero-weight directed path from the source vertex & to a free point a E A or zero-slack point b Bu in this residual network is an admissible path.
every D-free point in
Algorithm
In this section, we describe algorithm for computing optimal matching on point set AU B. Our algorithm builds randomly shifted quadtree Q. For any square D of the quadtree Q, we denote the set of children of O in Q by C[U].
We describe our divide-and-conquer algorithm with respect to an arbitrary square D in Q. Our algorithm computes D-optimal matching Mo, y(-).
Base case. If O is a leaf of Q, let v be the only point in D. If E B, set y(u) + d(u, D] and
MD < 0. If o E A, set y(v) 0 and Mp + 0 and return (see lines 1-4 of Algorithm
Divide step. If is not leaf of Q, for each child O' C[D), recursively compute the C' O'pptiiall matching v-) on Ag U B (see lines 6 and 7 of Algorithm D
Conquer step. For any child E' E C[J]. let MC', 3-) denote the C'-optimal matching returned by the algorithm. Set Mo = Uareqp Mo as the union of all matchings computed inside the children of O. Let B denote the -ffree points of Bn with respect Mn. To obtain C-optimal matching, our algorithm iteratively executes the CONSTR ANSTRAINEDHUNGE ARAAAARIIAAARCC procedure that adjusts the dual weights to find admissible path P with respect to Mo. Then, the AUGMENT procedure updates Mc by augmenting the matching along P. The algorithm continues to search and augment until there no U-free points of B_, i.e., until Mo v(-) is a O-optimal matching. See lines 8-11 11 and AUGMENT
of Algorithm|] | We describe the details of the CONSTRANDDHUNGGRRNNEENNN procedures below.
CONSTRAINED SONSTTAINEDHCNNGARANNSERRRE procedure: Compute the residual network with respect to Mn and execute Dijkstra's 's shortest path algorithm starting from the source s. Let K, denote the shortest path distance from s to as computed by the Dijkstra's algorithm. Let
K minf min Kas min Rib + (b))) aE4f
be B
and let denote the point that achieves this minimum. Note that E AD U Bn. Let Pa denote the shortest path from s to and denote the path obtained by removing from P.. For any E An U Bn, if K2 < K, update its dual weight to y(u) + y(v) + K..
2- AUGMENT procedure: Augment Mg along P.
Our algorithm returns the matching Mn- computed at the root D* of Q as the minimum-com matching between A and B. This completes the description of our algorithm. The pseudo-code of our divide-and-conquer algorithm is provided in Algorithm
3.1 Proof of correctness
Recollect that for the root square D of Q, from Lemma 2.2 D* -optimal matching is minimum- cost matching of A and B. Therefore, it suffices to show that, at each square D Q, our algorithm computes C-optimal matching.
We begin by showing this for any leaf square D. which by definition contains exactly one point. Since there are no edges inside D, conditions 44) and 5 hold trivially. Let v be the only point inside O.
7

---

Algorithm D&CHUNGARIAN An, Bu)
Input: A square D of quadtree Q and two point sets Ag and By inside
Output: A minimum-cost C-constrained matching Mn, 3())
1: if is leaf of Q then
b Base case
2:
Mu 0
3: Let v be the only point An U B
4: if' E An then y(v) < 0 else y(0) + d(v, [)
5:
else
6:
// Let 11, D2, 33 and D1 be the children of O in Q
D Divide step
7:
Mo. y f- D&CHUNGARIAN(D An, B,.), Vi E [1,4]
8:
Mn UkEt Ma.
Conquer step
9:
while B5 + do
10:
P. CONSTRAINEDEIINGRRINSSSACCCCCCC Ma + AUGMENT(Mg,P)
11:
12:
return Mn.y(.)
If B=, then y(u) d(v, [) and therefore. condition (6) holds. Otherwise, An and our algorithm sets y(v) = 0: therefore, condition (7) holds. In both cases, it is easy to see that B is empty and therefore the empty matching Mg and dual weight y(v) is a D-optimal matching.
For any non-leaf square of Q, at the beginning of the conquer step, the algorithm simply combines the constrained optimal matchings computed at its non-empty children. From Lemma 2.4 the resulting matching will be O-feasible matching. From Lemma bellow after the execution of the HONNTTRANNDD
procedure, the matching Mg vv-) remains a C-feasible matching and the path returned by the procedure is an admissible path. From Lemma 2.5 the execution of AUGMENT procedure then reduces the number of Ciffree points by one while maintaining the D-fcasibility conditions. Therefore, during the execution of the conquer step, the matching Mn. v-)) remains --feasble and the number of -ffree points reduces by one at each iteration; thus, the algorithm terminates with --optimal matching, as desired.
Lemma 3.1. Given any square of Q and any D-feasible matching Mn, v(),, after executing the procedure on D. the updated dual weights remain D-feasible. Furthermore, the returned path P is an admissible path.
3.2 Proof of Efficiency
For any square O. let To denote the execution time of the conquer step of our algorithm when executed on D. Additionally, for any i, let L(i) denote set of all squares of the quadtree that are processed at depth i of recursion and let T2 denote the total execution time of our algorithm across all such squares, i..., T, = E Th.
Recall that the conquer step of our algorithm consists of iterations, where each iteration, it executes the CONSTRAINEDI SANNTTAINEDHUNGAAIAASSLLL SEARCH procedure followed by the AUGMENT procedure. For any square D of Q, from Lemma 2.5) each iteration of the conquer step reduces the number Q-free points of Bu by one. Therefore, the number of iterations is bounded by the number of O-free points of BD, which at most total number free points with respect to the -MMM computed for all children C' of C. By invoking Lemma 2.2 on all children of [, we bound the expected number of iterations of the conquer step O by (n for dimensions (and ( 1- + for d dimensions).
2.2
Using any dynamic weighted nearest neighbor data structure with an update/query time of & (nm). one can execute the CONSTRANNEDHUNGGAIAAASSS For planar point sets, (( poly (log Ra) and therefore, each execution of the CONSTRAINED- HUNGARIANSEARCH procedure takes Ó(nm) time. Furthermore, the AUGMENT procedure simply augments the matching along an admissible path, which can be done in O(nn) time. Therefore, the conquer step of our algorithm executes ((n3/4) iterations for 2 dimensions (and O -a4) iterations for d dimensions), in expectation, where each iteration takes ((m)) time for 2 dimen- sions (and ((z=)(na)) for d dimensions). Since Eoer() no I E [T] Σoec E [Tb]
procedure in (ng time

---

200
Uniform vs. Uniform
(b) Uniform vs. Gaussian
(c) Pickup vs. Drop-off
Figure 5: The running time of algorithm (D&C Hungarian Algorithm) and the Hungarian algorithm for computing the 1-Wasserstein distance between μ and v, where (a) both μ and v are Uniform, (b) j is Uniform and is Gaussian, and (c) is Pickup and is Drop-off.
δ (EDecto) n /3/""m)) = (17/44 for 2 dimensions (and E [7.] = ,(24 ((n)) for d dimen-
sions). Summing all levels of Q, the expected execution time of our algorithm is On log A) for 2 dimensions (and Ō(2%n2 da (n) log A) for d dimensions). Using slightly sophisticated argument which is presented in the Appendix, we remove the 2 from the execution time for the d-dimensional settings, leading to Theorem 1.1
Remark 3.2. When the input points A J are chosen arbitrarily, the number of iterations of the conquer step can be ((n)), leading to an execution time of ((n))) per square and O(n² log A) in total.
4 Experiments
In this section, we present the results of our experiments comparing the execution time our algorithm to that of the Hungarian algorithm. Both algorithms are implemented in Java and share the same data structures For both algorithms, we use the classical implementation path algorithm used in the Hungarian search procedure and do not use any weighted nearest neighbor data structures. All computations are performed using a single calculation thread on computer with 2.6 GHz 6-Core Intel Core i7 CPU and 16 GB of 2667 MHz DDR4 RAM.
Dijkstra's 's shortest
Datasets. We test our algorithm on 2-dimensional synthetic and world datasets. For the synthetic dataset, we use two distributions, namely (i) a uniform distribution defined on the unit square (Uniform), and (ii) a Gaussian distribution constrained the unit square with a randomly chosen mean inside the unit square and a standard deviation of 0.25 (Gaussian). For real-world dataset, we employ the New York Taxi dataset [47] and obtain two distributions, namely (i) the distribution of pickup locations (Pickup) and (ii) the distribution of drop-off locations (Drop-off) of passengers. We
filtered the datasets by considering trips in seven dates in 2014 with (i) trip duration of at least minutes, and (ii) trip velocity of at most 1 110mph.
Tests. In each test, we conducted experiments using two sets of n i.i.d samples from distributions μe and within the unit square. We executed both our divide-and-conque algorithm and the Hungarian algorithm to compute the empirical I-sssserttein distance between the same set of samples and compared the execution times. In our first experiment, we used Uniform distribution for both μ and v to compare the performance of algorithms when the distributions are the same. For the second experiment synthetic datasets, we compared running times with different distributions Uniform and v being Gaussian, and in the third experiment, which is on the real-world datasets, j was set to the Pickup distribution while was the Drop-off distribution. In our experiments on the synthetic datasets, we recorded the number of iterations our algorithm executes at each square of the quadtree to find out how tight our analysis are. The results of our experiments shown in Figures and 6
being
Results. As depicted in Figure 5 experiments suggest that our algorithm outperforms the Hungarian algorithm, and the improvement is most significant when the distributions are the same (Figure 5))). More generally, experiments suggest that our algorithm performs significantly better than the Hungar-
Our implementations available http://ggghbbbb.oommgggaatttt/lt Euuersaeant-elaaeannaaaaaaaa Euclidean-Bipartieee Bipartite-

---

2003
1750
Sy10.53
Senb.ss
Number Swentaree
1392
0982
500
2:0
0
10000
20000
20000
00030
70000
10000
20000
30000
10000
S2000
(a) Uniform vs. Uniform
(b) Gaussian vs. Gaussian
Figure 6: The number of iterations of the conquer step for square with n points when executed on samples from the same distributions ((e) Uniform and (f) Gaussian).
ian algorithm when the matching cost is low (Figures 5 a) and (c)). The reason for this improvement lies in the fact that the Hungarian algorithm performs global Hungarian search' for matching each point, whereas our algorithm isolates shorter edges of the optimal matching within smaller quadtree sub-problems and execute 'local Hungarian search' for them. Therefore, while the asymptotic im- provements are shown for samples from the same distribution, we expect our algorithm to perform better when the optimal matching has many edges with small cost, such when points are drawn from two similar but not identical distributions. This is evident in our third experiment (Figure (J))), where we sample pick-up locations and n drop-off locations from the New York Taxi dataset. Notably, pick-up and drop-off locations tend to follow different distributions, with pick-ups having higher density around Manhattan.
Additionally, our algorithm handles larger point sets efficiently; for example, it computes the optimal solution for 50 000 points from the Uniform distribution approximately 700 seconds. Finally, in Figure O our experiments on samples from the same synthetic dataset suggest that for a square with n points, the number of iterations of the conquer step of our algorithm might be bounded by O(n") for a 0.55. This perhaps suggests that our upper-bound of (nz3/4) may be an overestimate, at least for the uniform and Gaussian distributions.
Conclusion
In this paper, we adapted the classical primal-dual approaches for computing minimum-coo cost matching within the divide-and-conquer framework. This approach led to sub-quadratic weakly-polynomia exact algorithm for computing minimum-cost matching stochastic point sets. Notably, our algorithm incorporates randomly shifted quadtree, a structure previously used only in approximation algorithms for p-rp-Ma WaWssserseeinn distance. We conclude by highlighting some open problems for further investigation.
An interesting question is whether it is possible further exploit the geometry of the dual weights and give tighter bound for the number of iterations of our algorithm at each square; thus, improving the analysis of our algorithm.
Our analysis only requires the fact that the average length of the matching edges in large cells is small. Stochastic point sets have this property and therefore, algorithms for them. Is there a sub-quadratic exact algorithm for the case where the matching edges long? Understanding this may shed light into the design of sub-quadratic exact algorithm for arbitrary inputs.
achieve sub-quadratic
Finally, our algorithm is weakly polynomial algorithm as its running time is dependent on the spread of the point set. Can we remove the dependence on spread while maintaining the simplicity of our algorithm?
10

---

Acknowledgement
We would like acknowledge Advanced Research Computing (ARC) at Virginia Tech, which provided us with the computational resources used to run the experiments. Research presented in this paper was funded by NSF CCF-1909171 and NSF CCF-2223871. We would like thank the anonymous reviewers for their useful feedback.
References
[1] Pankaj K Agarwal and Jeff M Phillips. On bipartite matching under the rms distance. In Proceedings of the Eighteenth Canadian Conference on Computational Geometry', Kingston, Canada, pages 143-146, 2006.
[2] Pankaj K. Agarwal and Kasturi R. Varadarajan. A near-linear constant-factor approximation for Euclidean bipartite matching? In 20th International Symposium on Computational Geometry, pages 247-252 2004.
[3] Pankaj K Agarwal, Alon Efrat, and Micha Sharir. Vertical decomposition of shallow levels in 3-dimensional arrangements and its applications. In Proceedings of the eleventh annual symposium Computational geometry, pages 39-50, 1995.
[4] Pankaj K Agarwal, Hsien-Chih Chang, Sharath Raghvendra, and Allen Xiao. Deterministic, near-linear e-approximation algorithm for geometric bipartite matching. In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 1052-1065, 2022.
[5] Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Rachita Sowle. A higher precision algorithm for computing the 1-wasserstein distance. In The Eleventh International Conference on Learning Representations, 2022.
[6] Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Rachita Sowle. An improved e-approximation algorithm for geometric bipartite matching. In /8th Scandinavian Symposium and Workshops on Algorithm Theory (SWAT 2022). Schloss Dagsuul-LLeibnii-ZZenrrn für Informatik, 2022.
[7] Miklós Ajtai, János Komlós, and Gábor Tusnády. On optimal matchings. Combinatorica, 4: 259-264,11 1984.
[8] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference machine learning, pages
PMLR, 2017.
[9] Jiezhang Cao, Langyuan Mo, Yifan Zhang, Kui Jia, Chunhua Shen, and Mingkui Tan. Multi- marginal wasserstein gan. Advances in Neural Information Processing Systems, 32, 2019.
[10] Sergio Caracciolo, Carlo Lucibello, Giorgio Parisi, and Gabriele Sicuro. Scaling hypothesis for the euclidean bipartite matching problem. Physical Review E, 9011):012118, 2014.
[11] Li Chen, Rasmus Kyng, Yang P Liu, Richard Peng, Maximilian Probst Gutenberg, and Sushant Sachdeva. Maximum flow and minimum-cost flow in almost-linear time. In 2022 IEEE 63rd Anmual Symposium on Foundations of Computer Science (FOCS), pages 612-623 623. IEEE, 2022.
[12] Marco Cuturi. Sinkhorn distances: Lightspeed computation optimal transport. In Advances in Neural Information Processing Systems 26, pages 2292-23000 230000 2013.
[13] Nabarun Deb and Bodhisattva Sen. Multivariate rank-based distribution nonparametric testing using measure transportation. Journal of the American Statistical Association, 118(541): 192-207, 2023.
[14] Nabarun Deb, Bhaswar Bhattacharya, and Bodhisattva Sen. Efficiency lower bounds for distribution-fie free hotelling-type two-sample tests based on optimal transport. arXiv preprint arXiv:2104.0 01 01966, 2021.
[15] Eustasio Del Barrio and Jean-Michel Loubes. Central limit theorems for empirical transportation cost in general dimension. The Annals of Probability, 47(2):926 951, 2019.
11

---

[16] Richard Mansfield Dudley. The speed of mean cdivenkkocantell convergence. The Annals of Mathematical Statistics, 40(1):40-50, 1969.
[17] RM Dudley. Real analysis and probability, 2002.
[18] Nghia Buong-Trunge Minh-Hoang Nguyen, and Hanh TH Nguyen. Clustering stability via concept-based nonnegative matrix factorization. In Proceedings of the 3rd International Confer- ence Machine Learning and Soft Computing, pages 49 54, 2019.
[19] Nicolas Fournier and Arnaud Guillin. On the rate convergence wasserstein distance the empirical measure. Probability theory and related fields, 162(3-4):707-733, 2015.
[20] Kyle Fox and Jiashuai Lu. near- linear time approximation scheme for geometric transporta- tion with arbitrary supplies and spread. In Proc. 36th Annual Symposium on Computational Geometry, pages 45:1-45:18 2020.
[21] H. N. Gabow and R.E. Tarjan. Faster scaling algorithms for network problems. SIAM J. Comput., 18:1013-1036, October 1989. ISSN 0097-5397.
[22] Harold N. Gabow and Robert Endre Tarjan. Faster scaling algorithms for general graph-matching problems. J. ACM, 38(4):815-853, 1991.
[23] Javier González-Delgado Alberto González-Sanz. Juan Cortés, and Pierre Neuvial. Two- sample goodness-of-fit tests the flat torus based wasserstein distance and their relevance to structural biology. arXiv preprint (eXiv:2108.0016: 2021.
[24] Xin Guo. Johnny Hong, Tianyi Lin, and Nan Yang. Relaxed wasserstein with applications to gans. ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3325-3329. IEEE, 2021.
[25] Masaaki Imaizumi, Hirofumi Ota, and Takuo Hamaguchi. Hypothesis test and confidence analysis with wasserstein distance on general dimension. Neural Computation, 34(6):1448 1487, 2022.
[26] Andrey Boris Khesin, Aleksandar Nikolov, and Dmitry Paramonov. Preconditioning for the geometric transportation problem. arXiv preprint arXiv: (rXiv:1902.000 083844 2019.
[27] Soheil Kolouri, Navid Naderializadeh, Gustavo K Rohde, and Heiko Hoffmann. Wasserstein embedding for graph learning. arXiv preprint arXiv:2006.0944 2020.
[28] Harold W Kuhn. The hungarian method for the assignment problem. Naval research logistics quarterly, 2(1-2):83-97, 97, 1955.
[29] Nathaniel Lahn and Sharath Raghvendra. on25/1) time e-approximation algorithm for rms matching in plane. In Proceedings of the 2021 ACM-SIAM Symposium Discrete Algorithms (SODA), pages 869-888. SIAM, 2021.
[30] Nathaniel Lahn, Sharath Raghvendra, and Kaiyi Zhang. A push-relabel based additive approxi- mation for optimal transport. arXiv preprint arXiv:2203.03732. 2022.
[31] Tilman Lange. Volker Roth, Mikio L Braun, and Joachim M Buhmann. Stability- validation of clustering solutions. Neural computation, 1696):1299-1323 2004.
[32] Richard J Lipton and Robert Endre Tarjan. Applications a planar separator theorem. SIAM journal computing, ((3):615-627, 1980.
[33] Huidong Liu, GU Xianfeng, and Dimitris Samaras. two-step computation of the exact gan wasserstein distance. In International Conference on Machine Learning, pages 3159-3168. PMLR, 2018
[34] Huidong Liu, Xianfeng Gu, and Dimitris Samaras. Wasserstein gan with quadratic transport cost. Proceedings of the IEEF/CVF international conference on computer vision, pages 4832-488 4841, 2019.
12

---

[35] Youssef Mrouch, Tom Sercu, and Vaibhava Goel. Megan: Mean and covariance feature matching gan. In International conference on machine learning, pages 2527-2535. PMLR. 2017.
[36] Victor M Panaretos and Yoav Zemel. Statistical aspects of wasserstein distances. Annual review of statistics and its application, 6:405 431, 2019.
[37] Hermina Petric Maretic, Mireille Gheche, Giovanni Chierchia, and Pascal Frossard. Got: an optimal transport framework for graph comparison. Advances in Neural Information Processing Systems, 32, 2019.
[38] Abhijeet Phatak, Sharath Raghvendra, Chittaranjan Tripathy, and Kaiyi Zhang. Computing all optimal partial transports. In The Eleventh International Conference on Learning Representa- tions, 2022.
[39] Julien Rabin, Julie Delon, and Yann Gousseau. Circular earth mover's distance for the compari- son of local features. In 2008 19th International Conference on Pattern Recognition, pages 1-4. IEEE, 2008.
[40] Sharath Raghvendra and Pankaj K. Agarwal. A near-linear time e-approximation algorithm for geometric bipartite matching. Journal of the ACM (JACM), 67(3):119, 2020.
[41] Aaditya Ramdas, Nicolás Gareía Trillos, and Marco Cuturi. On wasserstein two-sample testing and related families of nonparametric tests. Entropy, 19(2):47, 2017.
[42] Yossi Rubner, Carlo Tomasi, and Leonidas J Guibas. The earth mover's distance metric for image retrieval. International journal of computer vision, 40(2):99, 2000.
[43] R. Sharathkumar. A sub-quadratic algorithm for bipartite matching of planar points with bounded integer coordinates. In 29th International Symposium on Computational Geometry, pages 9-16, 2013.
[44] R Sharathkumar and Pankaj K Agarwal. Algorithms for the transportation problem in geometric settings. In Proceedings of the 23rd annual ACM-SIAM symposium on Discrete Algorithms, pages 306-317. SIAM, 2012.
[45] Peter Williston Shor. Random planar matching and bin packing. PhD) thesis, Massachusetts Institute of Technology, 1985.
[46] Max Sommerfeld, Jörn Schrieber, Yoav Zemel, and Axel Munk. Optimal transport: Fast probabilistic approximation with exact solvers. J. Mach. Learn. Res., 20< 2019.
[47] NYC Taxi and Limousine Commission (TLC). Trip record data. httpp://wwat http///wwwwwwg...../.. gov/ site/tlcc http:///aaaboottttttri---eeeeeee
(age 2023. Accessed: 2023-03-01.
[48] Pravin M Vaidya. Geometry helps in matching. SIAM Journal Computing, 1816):120-1122 1989.
[49] Jonathan Weed and Francis Bach. Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance. Bernoulli, 25(4A):2620 2648, 2019.
13


Reviews:
[{"review_id": "VkUQJuF042", "reviewer": "Reviewer_NeKH", "summary": "The paper considers matching two sets $A$ and $B$ of $n$ points in the Euclidean space so as to minimize the sum of distances of matched points, when both pointsets are drawn independently and identically from the same (unknown to the algorithm) distribution.  The authors extend the well-known Hungarian method with the shifted quad-tree decomposition technique, and improve over the best known result for the worst-case pointsets. The paper is complemented with experimental results.", "questions": "- I was convinced by your motivation for studying the problem, but I was wondering if you know of any real-world matching problem where it is natural to have the points drawn by the same distribution?", "limitations": "No limitations applicable and no potential negative societal impact.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is an interesting problem to consider and I found the extension of the Hungarian method with the shifted-quad tree technique quite interesting and -- at least to me -- novel. The paper is well written.", "weaknesses": "The result of the paper is novel but I am not sure if I would call a better performance for only a special class of inputs an improvement over (slightly worse) results but that hold for general inputs.\n\nI would have preferred if the paper explicitly stated it every time a technique from the literature is used/adapted.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thoughtful review and feedback. \n\n>  The result of the paper is novel but I am not sure if I would call a better performance for only a special class of inputs an improvement over (slightly worse) results but that hold for general inputs. \n\n**Response.** We do not claim a better performance but a more *robust* performance in comparison to the Hungarian algorithm. Our algorithm has a similar worst-case performance to the Hungarian method but a faster performance for stochastic point sets (see our general response for a discussion on this). \n\n> I was convinced by your motivation for studying the problem, but I was wondering if you know of any real-world matching problem where it is natural to have the points drawn by the same distribution?\n\n**Response.** There are many ML problems that require testing if two sample sets are from the same distribution.\n\n* Distributional shifts: Does the real data set represent the same distribution as the training data from which ML models are built? [1]\n* Benchmarks: Do models built using different ML techniques on the same training data represent the same distribution? [2, 3]\n* Mutual-independence: Given a multi-variate distribution, are the two marginals of the distribution independent? [4]\n\nThe questions above reduce to the *two-sample test problem*: Given two sets of $n$ samples, determine if they are drawn from the same multivariate distribution or different ones. In the *Wasserstein two sample test*, one checks if the optimal matching cost between the samples is below a threshold to determine if they are from the same distribution [5]. The computation of the optimal matching cost can be done using our algorithm.\n\n---\n---\n\n**References.**\n\n[1] S. Rabanser, S. Günnemann, and Z. Lipton. \"Failing loudly: An empirical study of methods for detecting dataset shift.\" Advances in Neural Information Processing Systems, 2019.\n\n[2] A. Borji. \"Pros and cons of gan evaluation measures.\" Computer vision and image understanding, 2019.\n\n[3] D. Lopez-Paz, and M. Oquab. \"Revisiting classifier two-sample tests.\" arXiv preprint, 2016.\n\n[4] N. Deb, and B. Sen. \"Multivariate rank-based distribution-free nonparametric testing using measure transportation.\" Journal of the American Statistical Association, 2023.\n\n[5] M. Imaizumi, H. Ota, and T. Hamaguchi. \"Hypothesis Test and Confidence Analysis With Wasserstein Distance on General Dimension.\" Neural Computation, 2022."}, {"review_id": "GdX25nOLvi", "reviewer": "Reviewer_FnYg", "summary": "This paper proposes a new, exact algorithm for solving the Euclidean weighted bipartite matching problem. Here, we have data sets $A,B \\subset \\mathbb{R}^d,$ each of cardinality $n,$ and the weight of an edge $ab$ is defined to be $\\lVert a - b \\rVert^p$ for any integer $p\\ge 1.$ This formulation is motivated by its direct application to computing empirical $p$-Wasserstein distances.\n\nFor data drawn i.i.d. from the same unknown distribution on the unit hypercube, the expected runtime of the algorithm is shown to be weakly-polynomial time, asymptotically in the number of points $n$ and with respect to an additional \\emph{spread} parameter $\\Delta $, defined to be the ratio of largest and smallest distances between any two points. As $d \\to \\infty $, this expected runtime approaches that of the classical Hungarian algorithm, up to polylog factors and assuming an efficient data structure for weighted nearest neighbors.\n\nMost of the analysis and all of the experiments focus on the special case where $d=2$, where the improvements over the Hungarian algorithm are more pronounced. The runtime analysis rests on the key observation of Lemma 2.1: if we cut a randomly-shifted square containing the data into four equal pieces, an asymptotically small number of data points will lie near the boundary where these pieces meet. A maximum-weight matching may be constructed by repeated subdivision: optimal matchings in each of the four pieces can be combined to get a feasible solution, and the algorithm searches for \"admissible\" augmenting paths until the feasible solution can no longer be improved. For a bipartite graph on $2n$ vertices drawn iid unformly from the unit square, the expected runtime of this algorithm is $\\tilde O (n^{7/4} \\log \\Delta )$. For comparison, the Hungarian algorithm in this setting is $\\tilde O (n^2)$. Experiments on real and synthetic data show that the new algorithm can substantially outperform the Hungarian algorithm.\n\nPost-rebuttal edit: The authors have responded in a satisfactory manner to my queries. Their algorithm is an interesting addition to the arsenal of matching algorithms, but nevertheless has certain limitations such as dependence on the parameter $\\Delta $. For this reason I maintain my high rating of accept.", "questions": "line 46: \" The execution time of our algorithm is similar to that...\" I assume Remark 3.2 is what justifies this remark? It would be helpful to include a forward reference here.\n\nline 59: \"and their fast implementations\" Don't you mean \"its fast implementations\"?\n\nline 322: \"The dataset consists of the locations\" This means latitude and longitude, correct?\n\nline 323: \"We filtered the datasets by considering trips\" Can you comment more on the choices behind your data filtering and other aspects of experimental design? Specifically, could you report some summary statistics (min, max, median, mean) for the spread of these data sets, both with and without the filtering? It is understandable that some filtering would be needed to see a performance increase over the Hungarian algorithm, and I think the thresholds are reasonable. Still, it would be helpful to understand the effect of filtering on both the runtime and the parameters used to analyze runtime. Moreover, given the filtering step, I think that the suggested improvement $\\tilde O( n^{.55})$ is very speculative.", "limitations": "For scientific limitations, refer to my comments under \"Weaknesses\". I do not see any potential for these results to have direct adverse effects on society.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Matching problems are foundational to computer science and optimization. This paper proposes a solution that demonstrates both significant theoretical and empirical advantages in the particular domain of empirical $p$-Wasserstein distance where these problems have been successfully applied. In the class of work proposing improvements to the Hungarian algorithm, the authors have overcome limitations in this domain have included focus on the unweighted case, distributional assumptions, approximate vs exact algorithms, and special cases of the geometric setting studied here. \n\nThe runtime analysis is presented very clearly, and theorems are stated carefully. Although the algorithm is ultimately only weakly polynomial, the authors theoretical assertions that the algorithm outperforms the Hungarian method for a range of parameters are backed up convincingly by empirical results.\n\nAnother noteworthy feature of the algorithm is that it does not rely on sophisticated data structures like prior work. Indeed, hardly any background at all is needed to understand the proofs and implement the algorithm. This is also reflected in the simplicity of the source code which the authors have shared to aid in reproducibility.", "weaknesses": "My main criticism, with a view towards the proposed application of $p$-Wasserstein distance computation, is that the case of _unequal_ distributions is treated mostly as an afterthought. Arguably, this is a more important case in practice since it involves less restrictive assumptions on the data. One finds a few comments addressing this scattered throughout the paper, eg. in the abstract and Remark 3.2. However, to justify the claim that the algorithm performs similarly to the Hungarian algorithm in such cases, it would be beneficial to give experimental results in addition to the asymptotic (as the authors have already done in the equal-distributions case.)\n\nI suggest generally that the authors be somewhat more explicit about the limitations of their algorithm. As already noted in the introduction, computing the empirical Wasserstein distance is mainly tractable in low dimensions. Nevertheless, other readers might be interested in the high-dimensional Euclidean matching problem. I find it difficult to believe that the performance of this method would be better than the purely-combinatorial Hungarian for $d$ large, particularly because of the high cost of dividing each sub-hypercube into $2^d$ pieces.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thoughtful review and constructive feedback. We address your concerns below. \n\n> My main criticism, with a view towards the proposed application of $p$-Wasserstein distance computation, is that the case of unequal distributions is treated mostly as an afterthought. Arguably, this is a more important case in practice since it involves less restrictive assumptions on the data. To justify the claim that the algorithm performs similarly to the Hungarian algorithm in such cases, it would be beneficial to give experimental results in addition to the asymptotic (as the authors have already done in the equal-distributions case.)\n\n**Response.** Thank you for this question. The case of *equal*  distribution also has significant applications in ML. We have highlighted some of these applications in our common response.\n\nWe have also conducted experiments for the unequal case. \nIn our experiment on the NY-Taxi data set (already included in our submission), we sampled $n$ drop-off locations and $n$ pick-up locations and match them. Based on our observation, the pick-up and drop-off locations tend to follow different distributions, where pick-ups seem to have a higher density around Manhattan than the drop-offs. \nFurthermore, we conducted additional experiments, where one set is drawn from a Gaussian distribution and the other set is chosen uniformly at random from a unit square. We also included the results of a similar experiment, where one set of points are drawn from the uniform distribution over the unit square and the other set are samples from a Guassian mixture model consisting of $5$ clusters in the $2$-dimensional space. See Figure 1 in the pdf document attached to our general response. We notice an improvement in the efficiency even in these cases.\n\n> I suggest generally that the authors be somewhat more explicit about the limitations of their algorithm. As already noted in the introduction, computing the empirical Wasserstein distance is mainly tractable in low dimensions. Nevertheless, other readers might be interested in the high-dimensional Euclidean matching problem. I find it difficult to believe that the performance of this method would be better than the purely-combinatorial Hungarian for large $d$, particularly because of the high cost of dividing each sub-hypercube into $2^d$ pieces.\n\n**Response.** The divide step of our algorithm takes only $O(dn)$ time.  Indeed, a cell has $2^d$ children, which may be much higher than $n$. However, we note that the only sub-problems that we care about are the ones created by the non-empty children, which are at most $n$, i.e., $O(n)$. One can create these sub-problems by simply scanning the points and placing them in the appropriate sub-problem in $O(dn)$ time.\n\n\n> \"We filtered the datasets by considering trips\" Can you comment more on the choices behind your data filtering and other aspects of experimental design? Specifically, could you report some summary statistics (min, max, median, mean) for the spread of these data sets, both with and without the filtering? It is understandable that some filtering would be needed to see a performance increase over the Hungarian algorithm, and I think the thresholds are reasonable. Still, it would be helpful to understand the effect of filtering on both the runtime and the parameters used to analyze runtime. Moreover, given the filtering step, I think that the suggested improvement $\\tilde{O}(n^{0.55})$ is very speculative.\n\n**Response.** The objective of applying filters is to eliminate erroneous entries in the data, such as entries of trips with negative duration or implausible velocity. To show that the effect of data cleaning step on the result of the experiment was insignificant, we re-executed our algorithm on the NY Taxi dataset, this time applying only two basic filters: (1) the trip duration had to be at least $3$ minutes, and (2) the trip velocity could not exceed $112 mph$. You can see the results in Figure 2 of the PDF file uploaded as part of our general response. "}, {"review_id": "zREGlRj0Sr", "reviewer": "Reviewer_SwWP", "summary": "This paper studies the Euclidean bipartite matching problem. In the problem, there is a complete bipartite graph on parts A and B, and the cost of edge (alb) is ||a-b||^p. The goal is to compute the minimum cost perfect matching as quickly as possible. This setting is most motivated in the paper by computing the empirical p-Wasserstein distance (though I think the problem is also interesting in its own right). It is important to note that these edge costs are non-integral, in particular they can even be irrational. Many of the sub-quadratic algorithms that work when the edge costs are integral do not apply here.\n\nIt was known that the minimum cost matching in this setting can be computed with the Hungarian algorithm in time O(n^3), and in time O(n^2) in geometric settings. The main contribution of this work is showing that in geometric setting (i.e., the vertices of A,B are points in R^d), there exists an algorithm that runs in weakly polynomial sub-quadratic time. The “weakly polynomial” comes from a term in the runtime of log(Delta), where Delta gives the spread (i.e., ratio of the distance between the farthest and closest pair of points in A and B) of the point set, specifically the runtime is \\tilde{O}(n^{2-1/(2d)}*phi(n)*log(Delta)), where the phi(n) is some runtime dealing with a nearest neighbor data structure (polylog(n) in d=2). Note that when A and B are sampled from some distribution mu in the unit square, even when mu is unknown, the runtime is \\tilde{O}(n^{7/4}*log(Delta)). \n\nThe main strategy is to deploy a geometric Hungarian algorithm with a divide and conquer method using the quadtree. However, the divide portion is only helpful when points in the quad tree are not close to the boundary of their children (in some sense that I’m stating very informally). Therefore the worst-case is still just that of the Hungarian algorithm, but we get improvement when these bad boundary cases do not happen.\n", "questions": "Other comments to authors:\n- (see question in weaknesses)\n- Great figures through out\n- Nice discussion of high level technical overview at the end of Section 1\n- Assuming you’re given another page in your camera ready, I’d add the conclusion section from your Appendix to the main body. The future directions/ open problems you give are nice.", "limitations": "na", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- Improved runtime on an interesting problem.\n- Interesting technical insights in using quadtrees with a geometric Hungarian algorithm.  Even after this main idea, there are a couple hurdles the authors need to overcome. \n- Techniques are relatively practical, and this is exemplified with experiments on both synthetic and real data. This is particularly interesting because a lot of the past work apparently uses complicated, hard to implement data structures.\n- I thought the paper was extremely well-presented. The technical difficulties of the paper are really well laid out. ", "weaknesses": "- It was not totally clear what techniques from the authors are totally new and what is building off of prior work. I would recommend adding a brief discussion of the work by Sharathkumar, at least. My understanding is that the big technical idea of, and then modifying it so the divide and conquer parts actually give progress, is a new perspective. Can you comment more on this?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thorough review and feedback.\n\n> It was not totally clear what techniques from the authors are totally new and what is building off of prior work. I would recommend adding a brief discussion of the work by Sharathkumar, at least. My understanding is that the big technical idea of, and then modifying it so the divide and conquer parts actually give progress, is a new perspective. Can you comment more on this?\n\n**Response.** Your understanding of the main technical idea is correct. The novelty of our algorithm is in placing the classical Hungarian algorithm within a geometric divide-and-conquer framework.\nOur algorithm does not rely on any prior work except for a weighted nearest-neighbor based implementation of the Hungarian search step (which was proposed in [1, 2] and used in [3, 4]). Our algorithm also uses a randomly shifted quadtree, which is a popular data structure for the design of geometric approximation algorithms. To the best of our knowledge, our algorithm is the first to use them in the design of an exact Euclidean bipartite matching algorithm.\n\n**Comparison with Sharathkumar [4].** We have included a brief comparison of our result with the result of Sharathkumar [4] in lines 61--66 of our initial submission. We will extend it to also include a comparison of techniques. \n\nFrom a technical standpoint, apart from using a combinatorial primal-dual approach, the algorithm by Sharathkumar [4] and our algorithm are quite different. \nThe algorithm in [4] uses the cost scaling framework of Gabow and Tarjan [5, 6] to find an approximate solution. Using the properties of this approximate solution, they 'trap' the edges of the optimal matching in a *planar* graph, i.e., a graph that can be drawn on a plane without any overlapping edges. They then use an algorithm by Lipton and Tarjan [7] to find an optimal matching inside this planar graph.\n\nThe proof of correctness of [4] relies on (a) the $2$-dimensional geometry of the input, (b) the edge-costs being square-roots of integers (owing to integer coordinates) and (c) the triangle inequality of Euclidean distances. Thus, their algorithm does not extend to (i) higher dimensions (due to (a) not being satisfied), (ii) to points with real-valued coordinates (due to (b) not being satisfied) or (iii) to the case where edge costs are squared-Euclidean (due to (c) not being satisfied). In contrast, our algorithm extends to all these cases.\n\n---\n---\n**References.**\n\n[1] P. Vaidya. \"Geometry helps in matching.\" In Proceedings of the twentieth annual ACM symposium on theory of computing, 1988.\n\n[2] R. Sharathkumar, and P. K. Agarwal. \"Algorithms for the transportation problem in geometric settings.\" In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, 2012.\n\n[3] P. K. Agarwal, A. Efrat, and M. Sharir. \"Vertical decomposition of shallow levels in 3-dimensional arrangements and its applications.\" In Proceedings of the eleventh annual symposium on Computational geometry, 1995.\n\n[4] R. Sharathkumar. \"A sub-quadratic algorithm for bipartite matching of planar points with bounded integer coordinates.\" In Proceedings of the twenty-ninth annual symposium on Computational geometry, 2013.\n\n[5] Gabow, Harold N., and Robert E. Tarjan. \"Faster scaling algorithms for network problems.\" SIAM Journal on Computing 18, no. 5 (1989): 1013-1036.\n\n[6] Gabow, Harold N., and Robert E. Tarjan. \"Faster scaling algorithms for general graph matching problems.\" Journal of the ACM (JACM) 38, no. 4 (1991): 815-853.\n\n[7] Lipton, Richard J., and Robert Endre Tarjan. \"Applications of a planar separator theorem.\" In 18th Annual Symposium on Foundations of Computer Science (sfcs 1977), pp. 162-170. IEEE, 1977.\n"}, {"review_id": "wJdlOxwbjc", "reviewer": "Reviewer_L63E", "summary": "The paper presents a randomized algorithm for computing a minimum-cost matching of a bipartite graph induced by two point-sets, A and B, in the Euclidean plane. The best running time for the problem under consideration is n^2polylog(n). When A and B are drawn independently and identically from a fixed probability distribution, the presented algorithm achieves a running time of n^{7/4}polylog(n+D), where D is the ratio between the distance of the farthest pair of points to that of the closest pair of points. The presented algorithm can potentially improve the quasi-quadratic bound under the probability assumption and the assumption that D is not too large. The algorithm generalizes to higher dimensions.\n\n", "questions": "Could you please explain the relevance of your results to NeurIPS, beyond the importance of the matching problem as a generic optimization problem?", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The minimum-cost matching is an important combinatorial problem, even on bipartite geometric graphs.", "weaknesses": "\n* The improvement over existing bounds is only under certain probabilistic assumptions, and assumptions about the relative positions of the points in the plane. Moreover, the presented algorithm is randomized.\n* I don’t believe that the results will appeal to the broad audience of NeurIPS. The paper does not motivate the problem well with respect to the scope of NeurIPS. The problem is very restricted.\n* The techniques are not very novel; they are adaptations of the Hungarian method plus a clever use of data structures.  \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank you for reviewing our submission.\n\n> Could you please explain the relevance of your results to NeurIPS, beyond the importance of the matching problem as a generic optimization problem?\n\nMinimum-cost bipartite matching is extensively used in many applications in Machine Learning, Computer Vision, and Statistical Inference. We enumerate a few of these below. \n\n**Machine learning applications:** Optimal bipartite matching and its cost is used (i) as a loss function in GANs [1, 2], (ii) within autoencoders [3, 4], (iii) in domain adaptation [5, 6], (iv) clustering [7], and (v) self-supervised learning [8].\n\n**Computer vision applications:** Multi-object tracking [9],  object-centric learning [10, 11], object detection [12, 13], image retrieval [14], instance segmentation [15, 16], and vector graphics [17].\n\n**Statistical inference applications:** Two-sample test [18, 19], mutual-independence test [20], and distributional shifts [21].\n\nComputing exact matchings is too expensive in many of these applications. This has motivated machine learning researchers to design approximation algorithms [22--28]. Our contribution is to obtain asymptotically faster exact algorithm for stochastic inputs such as those that arise in some of the applications above.\n\n---\n---\n\n**References.**\n\n[1] H. Liu, G. U. Xianfeng, and D. Samaras. \"A two-step computation of the exact gan wasserstein distance.\" In ICML, 2018.\n\n[2] J. Cao, L. Mo, Y. Zhang, K. Jia, C. Shen, and M. Tan. \"Multi-marginal wasserstein gan.\" NeurIPS, 2019.\n\n[3] A. Kosiorek, S. Sabour, Y. W. Teh, and G. E. Hinton. \"Stacked capsule autoencoders.\" NeurIPS, 2019.\n\n[4] S. Kolouri, P. E. Pope, C. E. Martin, and G. K. Rohde. \"Sliced Wasserstein auto-encoders.\" In ICLR, 2018.\n\n[5] Y. Balaji, R. Chellappa, and S. Feizi. \"Robust optimal transport with applications in generative modeling and domain adaptation.\" NeurIPS, 2020.\n\n\n[6] C. Lee, T. Batra, M. Haris Baig, and D. Ulbricht. \"Sliced wasserstein discrepancy for unsupervised domain adaptation.\" In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2019.\n\n[7] X. Yang, C. Deng, K. Wei, J. Yan, and W. Liu. \"Adversarial learning for robust deep clustering.\" NeurIPS, 2020.\n\n[8] X. Wen, B. Zhao, A. Zheng, X. Zhang, and X. Qi. \"Self-supervised visual representation learning with semantic grouping.\" NeurIPS, 2022.\n\n[9] Y. Zhang, P. Sun, Y. Jiang, D. Yu, F. Weng, Z. Yuan, P. Luo, W. Liu, and X. Wang. \"Bytetrack: Multi-object tracking by associating every detection box.\" In ECCV, 2022.\n\n[10] F. Locatello, D. Weissenborn, T. Unterthiner, A. Mahendran, G. Heigold, J. Uszkoreit, A. Dosovitskiy, and T. Kipf. \"Object-centric learning with slot attention.\" NeurIPS, 2020.\n\n[11] J. Brady, R. S. Zimmermann, Y. Sharma, B. Schölkopf, J. von Kügelgen, and W. Brendel. \"Provably Learning Object-Centric Representations.\" arXiv preprint, 2023.\n\n[12] Y. Wang, and J. M. Solomon. \"Object dgcnn: 3d object detection using dynamic graphs.\" NeurIPS, 2021.\n\n[13] Y. Li, Y. Chen, X. Qi, Z. Li, J. Sun, and J. Jia. \"Unifying voxel-based representation with transformer for 3d object detection.\" NeurIPS, 2022.\n\n[14] Y. Rubner, C. Tomasi, and L. J. Guibas. \"The earth mover's distance as a metric for image retrieval.\" International journal of computer vision, 2000.\n\n[15] B. Dong, F. Zeng, T. Wang, X. Zhang, and Y. Wei. \"Solq: Segmenting objects by learning queries.\" NeurIPS, 2021.\n\n[16] S. Hwang, M. Heo, S. W. Oh, and S. J. Kim. \"Video instance segmentation using inter-frame communication transformers.\" NeurIPS, 2021.\n\n[17] A. Carlier, M. Danelljan, A. Alahi, and R. Timofte. \"Deepsvg: A hierarchical generative network for vector graphics animation.\" NeurIPS, 2020.\n\n[18] M. Imaizumi, H. Ota, and T. Hamaguchi. \"Hypothesis Test and Confidence Analysis With Wasserstein Distance on General Dimension.\" Neural Computation, 2022.\n\n[19] N. Deb, B. B. Bhattacharya, and B. Sen. \"Efficiency lower bounds for distribution-free Hotelling-type two-sample tests based on optimal transport.\" arXiv preprint, 2021.\n\n[20] N. Deb, and B. Sen. \"Multivariate rank-based distribution-free nonparametric testing using measure transportation.\" Journal of the American Statistical Association, 2023.\n\n[21] S. Rabanser, S. Günnemann, and Z. Lipton. \"Failing loudly: An empirical study of methods for detecting dataset shift.\" NeurIPS, 2019.\n\n\n[22] M. Cuturi. \"Sinkhorn distances: Lightspeed computation of optimal transport.\" NeurIPS, 2013.\n\n[23] J. Altschuler, J. Niles-Weed, and P. Rigollet. \"Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.\" NeurIPS, 2017.\n\n[24] P. Dvurechensky, A. Gasnikov, and A. Kroshnin. \"Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn’s algorithm.\" In ICML, 2018.\n\n[25] G. Luise, A. Rudi, M. Pontil, and C. Ciliberto. \"Differential properties of sinkhorn approximation for learning with wasserstein distance.\" NeurIPS, 2018.\n\n[26] N. Lahn, D. Mulchandani, and S. Raghvendra. \"A graph theoretic additive approximation of optimal transport.\" NeurIPS, 2019.\n\n\n[27] J. Altschuler, F. Bach, A. Rudi, and J. Niles-Weed. \"Massively scalable Sinkhorn distances via the Nyström method.\" NeurIPS, 2019.\n\n[28] P. K. Agarwal, S. Raghvendra, P. Shirzadian, and R. Sowle. \"A Higher Precision Algorithm for Computing the $1 $-Wasserstein Distance.\" In ICLR, 2022."}]

Author rebuttal:
Thank you for the very positive feedback on our work. We want to emphasize a few important points that were also presented in our manuscript and hope that these points may help address some of the reviewers' criticisms.

**Novelty.** The novelty of our algorithm is in placing the classical Hungarian algorithm within a geometric divide-and-conquer framework. Hungarian algorithm conducts a "global Hungarian search" to match each point. In contrast, our algorithm obtains speed-up by trapping shorter edges of the optimal matching in smaller sub-problems (squares) of the quadtree and matching these points using a Hungarian search that is local to the sub-problem. While our asymptotic improvements are shown for stochastic points that are drawn from the same distribution, we do expect our algorithm to perform better when the optimal matching has many edges with small cost, for instance, when points are drawn from two similar but unequal distributions. This is best exemplified by our experiment on NY-Taxi data set (already included in our submission), where we sample $n$ drop-off locations and $n$ pick-up locations and match them. Based on our observation, the pick-up and drop-off locations tend to follow different distributions, where pick-ups seem to have a higher density around Manhattan than the drop-offs. We also conducted two additional experiments on samples drawn from two different distributions (the first using Gaussian and uniform distributions, and the second using a Gaussian mixture model with $5$ clusters and a uniform distribution) and include the results in the one-page pdf submitted as part of our response.

In terms of techniques, our algorithm does not rely on any prior work except for a weighted nearest-neighbor based implementation of the Hungarian search step (which was proposed in [1, 2] and used in [3, 4]). Our algorithm also uses a randomly shifted quadtree, which is a popular data structure for the design of geometric approximation algorithms. To the best of our knowledge, our algorithm is the first to use them in the design of an exact Euclidean bipartite matching algorithm.

**Applications.** There are many ML problems that require testing if two sample sets are from the same distribution.

* Distributional shifts: Does the real data set represent the same distribution as the training data from which ML models are built? [5]
* Benchmarks: Do models built using different ML techniques on the same training data represent the same distribution? [6, 7]
* Mutual-independence: Given a multi-variate distribution, are the two marginals of the distribution independent? [8]

The questions above reduce to the *two-sample test problem*: Given two sets of $n$ samples, determine if they are drawn from the same multivariate distribution or different ones. In the *Wasserstein two sample test*, one checks if the optimal matching cost between the samples is below a threshold to determine if they are from the same distribution [9]. The computation of the optimal matching cost can be done using our algorithm.
Two-sample testing dates back to Leighton and Shor [10], who used it to evaluate the quality of pseudo-random number generators.

Several ML problems generate samples from similar (but not identical) distributions, such as the domain adaptation problem [11] and the training of a Wasserstein GAN [12]. In both instances, our algorithm can be faster than Hungarian when the generated samples are low-dimensional.

---
---
**References.**

[1] P. Vaidya. "Geometry helps in matching." In Proceedings of the twentieth annual ACM symposium on theory of computing, 1988.

[2] R. Sharathkumar, and P. K. Agarwal. "Algorithms for the transportation problem in geometric settings." In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, 2012.

[3] P. K. Agarwal, A. Efrat, and M. Sharir. "Vertical decomposition of shallow levels in 3-dimensional arrangements and its applications." In Proceedings of the eleventh annual symposium on Computational geometry, 1995.

[4] R. Sharathkumar. "A sub-quadratic algorithm for bipartite matching of planar points with bounded integer coordinates." In Proceedings of the twenty-ninth annual symposium on Computational geometry, 2013.

[5] S. Rabanser, S. Günnemann, and Z. Lipton. "Failing loudly: An empirical study of methods for detecting dataset shift." Advances in Neural Information Processing Systems, 2019.

[6] A. Borji. "Pros and cons of gan evaluation measures." Computer vision and image understanding, 2019.

[7] D. Lopez-Paz, and M. Oquab. "Revisiting classifier two-sample tests." arXiv preprint, 2016.

[8] N. Deb, and B. Sen. "Multivariate rank-based distribution-free nonparametric testing using measure transportation." Journal of the American Statistical Association, 2023.

[9] M. Imaizumi, H. Ota, and T. Hamaguchi. "Hypothesis Test and Confidence Analysis With Wasserstein Distance on General Dimension." Neural Computation, 2022.

[10] F. T. Leighton, and P. Shor. "Tight bounds for minimax grid matching, with applications to the average case analysis of algorithms." In Proceedings of the eighteenth Annual ACM symposium on theory of computing, 1986.

[11] Y. Balaji, R. Chellappa, and S. Feizi. "Robust optimal transport with applications in generative modeling and domain adaptation." NeurIPS, 2020.

[12] H. Liu, G. U. Xianfeng, and D. Samaras. "A two-step computation of the exact gan wasserstein distance." In ICML, 2018.

Meta-review: This paper reduces the runtime complexity of bipartite matching, when the points in the two partitions are drawn independently and random from unknown distribution. The authors discuss some applications like NY-Taxi data set (already included in our submission), where they sample drop-off locations and pick-up locations and match them.
Decision: Accept (poster)