Fully Dynamic Euclidean Bi-Chromatic Matching in Sublinear Update Time
Authors: Gramoz Goranci
Peter Kiss
Neel Patel
Martin P. Seybold
Eva Szilagyi
Da Wei Zheng
Venue: ICML 2025 oral (2025)
Keywords: Euclidean bi-chromatic matching
dynamic algorithm
1-Wasserstein distance
Primary area: general_machine_learning->everything_else

Abstract: We consider the Euclidean bi-chromatic matching problem in the dynamic setting, where the goal is to efficiently process point insertions and deletions while maintaining a high-quality solution. Computing the minimum cost bi-chromatic matching is one of the core problems in geometric optimization that has found many applications, most notably in estimating Wasserstein distance between two distributions. In this work, we present the first fully dynamic algorithm for Euclidean bi-chromatic matching with sublinear update time. For any fixed $\varepsilon > 0$, our algorithm achieves $O(1/\varepsilon)$-approximation and handles updates in $O(n^{\varepsilon})$ time.  Our experiments show that our algorithm enables effective monitoring of the distributional drift in the Wasserstein distance on real and synthetic data sets, while outperforming the runtime of baseline approximations by orders of magnitudes.

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublinear Update Time
Gramoz Goranci Peter Kiss?
Neel Patel
Martin P. Seybold 4 Eva Szilagyi David Zheng
Abstract
points in Rd and set B of n blue points in R4, the goal is to compute a bijection μ A - B such that Eae alla p(a))l2 is minimized. A key application of minimum-cost bi-chromatic matching in machine learning is estimating the 1-Wasserstein distance of two spatial distributions v and v' whose supports satisty Sy, S E gdd This distance, also known as Euclidean Earth Mover's Distance (EMD), Kantorovich-Rubieet metric, or Mallow's distance, is defined
We consider the Euclidean bi-chromatic match- ing problem in the dynamic setting, where the goal is to efficiently process point insertions and deletions while maintaining a high-quality solu- tion. Computing the minimum cost bi-chromatic matching is one of the core problems in geo- metric optimization that has found many appli- cations, most notably in estimating Wasserstein distance between two distributions. In this work, we present the first fully dynamic algorithm for Euclidean bi-chromatic matching with sub-linear update time. For any fixed & > O, our algorithm achieves -aproxiaiiooi and handles up- dates in On time. Our experiments show that our algorithm enables effective monitoring of the distributional drift in the Wasserstein distance on real and synthetic data sets, while outperforming the runtime of baseline approximations by orders of magnitudes.
inf Js
W (1,,))) PET(U.)) ss xS,,
||a (la-blla b1l2 dP(a,
were [(v, v') denotes the set of all probability measures on So x Ser with marginal distributions v and v/.
Estimating the 1-Wasserstein distance measure, or the Eu- clidean bi-chromatic matching problem, has been exten- sively studied across a wide range of disciplines, including machine learning (Liu et al. 2018 2017 2018 Cao et al.,
2019 Balaji et al. 2020 Tolstikhin et al. 2018 computer vision (Rbbner et al. 2000), statisties (Panaretos and Zemel, 2019 and economics (Galichon 2018). This broad inter- est has led different problem varaiants, e.g., the point reweighted setting, and the development of algorithms with high-accuracy approximation guarantees or those that only estimate the solution size without explicitly outputting the underlying matching (for a more detailed discussion, please see Appendix[A) AA).
1. Introduction
We consider the Euclidean bi-chromatic matching problem, one of the most fundamental optimization problems in low- dimensional Euclidean geometry. Given a set A of n red
"Equal contribution 'Faculty of Computer Science, University of Vienna, Austria. 'Faculty of Computer Science, University of Vienna, Austria. This research was funded in whole in part by the Austrian Science Fund (FWF) 10,55776/ESP0088024, University of Southern California, CA, US 'Faculty of Computer Science, University of Vienna, Austria. This project has received funding from the European Research Council (ERC) under the Eu- ropean Union's Horizon 2020 research and innovation programme (Grant agreement No. 101019564 "The Design of Modern Fully
In this work, consider the Euclidean bi-chromatic match- ing in the dynamic setting, where the input undergoes up- dates such as insertions or deletions of pairs of points. The goal is to process these updates as fast as possible while maintaining a matching that achieves a small, provable ap- proximation ratio compared to the optimal matching at any given time.
Dynamic Data Structures (MoDynStruer)" Austrian Science Fund (FWF) project "Fast Algorithms for a Re- active Network Layer (ReactNet)", P 33775-N, with additional funding from the netidee SCIENCE Stiftung, 2020-2024. 'Faculty of Computer Science UniVie Doctoral School Computer Science DoCS, University of Vienna, Austria. "University of Illinois at Sibana-Champali IL, US. Correspondence to: Eva Szilagyi <cva.szilagyi@@ ietieeyi@@niiii..em act.te
and from the
The dynamic setting is particularly relevant in the statisti- cal estimation of the W-Wasserstein distance (Gattani et al., 2023 Beugnot et al.) 2021 Liu al., 2018 Panaretos and Zemel 2019), when the closed-form expressions for // and v unknown. In such cases, one efficiently draw independent and identically distributed (i.i.d.) samples from both distributions. The resulting empirical distribution in and 1/,, assign a uniform probability mass of 1/n to each sample in A and B, respectively. It is well known that W1 (Vn, v/)) converges to W, (v, 2'), as n + xc. Rather
Proceedings of the 4nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
than computing a matching from scratch after new samples drawn, dynamic algorithm guarantees efficient mainte- nance of accurate estimates for min, # E - u(a) lla under such changes. Further real-world applications of dy- namic geometric matching include measuring the similarity between evolving data sets (Alvarez-Melh) and Fusi 2020), measuring the change in patient data (such as MRI images)
dynamic graphs.
Furthermore, dynamic graph matching algorithms focus on singular edge updates. As in the dynamic Euclidean bi- chromatic matching problem we have to work with updates affecting (?(n) point-wise distances, in order to obtain our results, we need techniques which are new to dynamic lit- erature and heavily exploit the guarantees of the geometric setting.
over a long period of time (Gramfort et al., 2015 Janati
et al., 2019). or applications of matching (or EMD) as
metric in time series analysis (Cheng et al.] 2021). Despite
To complement our primary contribution (Theorem we present experimental evaluations of the performance of our dynamic algorithm for estimating the 1- 1-Wasserstein distance in contrast to static state of the art approximate algorithms. Our experimental results reveal that in the dy- namic setting the total running time of our algorithm over a series of updates outperforms periodic static recomputation by orders of magnitudes. Furthermore, we find that the ap- proximation ratio obtained by our algorithm on real-world datasets is smaller than that of our theoretical guarantees and almost matches that of static algorithms. Our results confirm the expected approximation ratio drop-off and run- ning time decrease we would expect with the reduction of the e parameter.
the fundamental importance of this problem, dynamic geo- metric matching has only recently been considered in (Xu) and Ding 2024a) and the existing algorithms cannot break the linear barrier on the update time.
1.1. Main Contributions
We obtain the first sub-linear dynamic algorithm for Eu- clidean bi-chromatic bipartite matching. Specifically, our algorithm achieves (1) approximation in sub-linear in n update time. The exact guarantees of our result are summa- rized in the theorem below.
Theorem 1.1. For any D 1. there exists a fully dynamic algorithm that maintains an expected ((1/E)- approximate solution to the Euclidean bi-chromatic match- ing problem defined on the point-sets A, B C IR2,, |A| = [B), while pairs of points are inserted into and deleted from A, B in On , - worst-case update time. Here, n denotes the maximum size of A throughout the update sequence.
1.2. Technical Overview
In vein of related works, we discuss algorithms and analysis for d = 2 to simplify the presentation, though our theorems and implementation extend to constant d > 2.
Our dynamic algorithm builds on the static algorithm of (Agarwal and Varadarajan 2004). On very high level, the algorithm partitions the input points into a series of nested grid cells. The nested grid is constructed in a such way that each grid cell contains O(n²) sub-cells of the next level, bottom level cells contain at most (nn)) points and the data-structure consists of O(=-t) levels. Intuitively, the algorithm aims to match every point within the smallest grid cell possible, however, certain cells might not have the same number of red and blue points.
In the static setting, the Euclidean bi-chromatic matching problems admits known (1 + -approxximai algorithms that run in near-optimal running time (Agarwal et al. 2022b Raghvendra and Agarwal 2020). In comparison, the ((1/E) approximation guarantee in Theorem 1.1 may initially seem less competitive. However, as we show in Appendix Appendix[[[The (The- (F.I. simple lower bound construction reveals that, in the dynamic setting (even if restricted to only point inser- tions or deletions), there is no algorithm that simultaneously achieves (2- - 6) apprroxximtioonn and runs in sub-linear up- date time, for any 8 > O. This in turn implics that no non-trivial dynamic algorithm can maintain a solution to the problem with an arbitrarily tight approximation ratio.
A crucial difference between our implementation and that of (Agarwal and Varadarajan| (004) is that we describe process which constructs the matching in a bottom-up manner, starting with the smallest cells and progressing towards the cell containing all of the input, whereas the prior work iterates in top-to-botton order. This change allows for a shorter algorithm description; however, it requires a different but arguably simpler analysis, which we present in Appendix C
While the matching problem on graphs has received signifi- cant attention in the dynamic model (see (Azarmehr et al.] 2024) for a complete list of works), the only prior known algorithm for the geometric version of the problem was presented by (Xu and Ding 2024a) with update time (n)).
Notably, even disregarding edge weights, the maintenance of a perfect matching in the graph setting requires 2(n) update time under well accepted hardness conjectures (Hen- zinger et al.] 2015) Dahlgaard| (0166). On the other hand, there is a long line of work focusing on maintaining an ((1) apprrxxixe matching in near optimal update time in
As bottom cells contain only small number of input points, they may select a maximal color-balance subset of them, match its points optimally using a polynomial running time static algorithm, and forward the remaining points to their parent cell.
2

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
If a set of points can't be matched within a cell due color dis- dis-balance, can argue that an optimal cost matching must match the dis-balance with points outside of the cell, and hence with longer edges. This allows us to introduce some slack and represent these forwarded points implicitly by the center of their respective cell and their cardinality.
Appendix A
Paper Outline: We describe static algorithm and its dynamic implementation in sections [ and E| respectively. In Section we demonstrate the practicality and efficiency of the dynamic algorithm for estimating the 1- 1-Wasserstein distance.
Larger cells aiming to match these forwarded implicitly represented points can efficiently find an implicit matching of the underlying points by an optimal polynomial run- ning time geometric transportation algorithm. These im- plicit matchings can then efficiently be turned into an actual matching of the underlying points of roughly the same cost.
2. Preliminaries
The Euclidean matching problem Given a set of red points A and set of blue points B in the plane, where |A| = [ n, matching is a bijection of the form #: A B. The cost of matching p. is given by c(j) = EeEA ||a - u(a)]l2 The objective is to find a matching of minimum cost. We say that a pair of points e == (a, b), referred to as an edge, belongs to matching μ if b - p(a).
To dynamize the algorithm, we show that, due to an update, the matching structure might only change in cells containing the updated points. After an update for each affected cell in a bottom-up manner we carefully adjust the matchings calculated. Since in each cell we either have O(n)) points or a larger set of points represented with O(n=) sub-cell center points, this can be done with polynomial running time static sub-routines for each affected cell. As the nested grid contains O(s-1) levels, the update time works out to Onn=
We define matching μ to be a-approximate for 2 1 if c(u) a S c(u²) for any optimal solution μ²".
To simplify presentation, we assume that X = C x D for some = poly where D a power of 2. Throughout this paper we make the standard assumption that the spread of the input points, i.e. the ratio between the largest and the smallest inter-point distance in X, remains bounded by some large U poly(n).
Updating the matching of bottom cell is straightforward, as we may run an optimal polynomial time algorithm on their updated points. However, the sub-cell center-point of intermediate cells may represent ((n) input points implic- itly. Hence, if the implicit matching of the cell significantly changes due to an update, it may take ((n) time to update the actual output matching within the cell.
A generalization of the matching problem, supporting point multiplicities, is the so-called Euclidean transportation problem, defined as follows.
The Euclidean transportation problem Let A, B C RR? be any two sets of points such that each point a E A has a non-negative integer supply 8a 2 0 and each E B has a non-positive integer demand d, I 0 satisfying E.EA 8a +
To overcome this difficulty, we first carefully analyze the update process and show that the points matched in each cell may only change by 2 due to an update. Afterwards, we present graph shortest s path sub-routine based procedure which calculates updated version of the geo- metric transport solution stored in the cell. This solution only differs in at most Onz²) edges compared to its unup- dated state. In turn, this allows us to efficiently update the matching of the actual input points implicitly represented by the geometric transportation sub-solutions.
EsE de O. An assignment, y : A B - Zzu, is a non-negative set of weights on the edges such that 8a
ELEB (a, b), Va E A and -d, = Eaes (a, b), vb B.
The cost of an assignment y is c(7) EcEA r(a, bi]la - b112- An optimal assignment for the Euclidean transportation problem is an assignment of minimum cost.
Theorem 2.1. (Atkinson and Vaidya 1995) There ex- ists deterministic algorithm which, given an instance of the Euclidean transportation problem A, C. R2, s E ZA d E Z returns an optimal assignment 7
1.3. Related Work
Numerous results exist on special cases and generalizations of Euclidean bi-chromatic bipartite matching across various computational settings. In the static setting, an exact solu- tion can be found through minimum cost max-flow graph algorithms in Ó(n²) time (yan den Brand et al. 2020p. For restricted cases of the problem. (Sharathkumar 2013) has presented algorithm with 3/2) running time. In the approximate setting, (Agarwal et al. 2022a) presented the first deterministic (n)) time algorithm for arbitrarily tight approximations, Due to space restrictions, a more detailed overview and discussion on related work can be found in
B - Z:o in time (n log(n log M)), where M max a1100 [[d]l n n- IAU
The results of this paper rely a standard data structure for representing points in geometric setting (see (Agarwal and Varadarajan 2004 Agarwal et al. 2022al), the under- lying and restricted p-trees, two quadtree- like structures. Throughout the paper think of the parameter p as some n sized power of 2.
3

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
The underlying and the restricted p-tree We define cell of side-length L > 0 as square of the form (a, a L)) x [6, 6 + L) for some a, b E R. We will frequently subdivide cells of side-length L into subcells of side length L/p. We will call this a grid of side length T./p.
dren nodes in T. Starting from the root of T. traverse the children cells of each internal node row by row, then recursively traverse the subcells of each cell until a subcell contains only one point.
For integer p= 2 which is also power of 2, and point set X, we define underlying p-tree T on X to be a quadtree-like structure with larger branching factor (see e.g., (Agarwal and Varadarajan 2004), ) formally described as follows.
3. Basic Framework and the Static Algorithm
The starting point of our work is the static algorithm for computing bi-chromatic Euclidean matchings due to (Agar wal and Varadarajani 2004). Specifically, their top-down algorithm computes an -aprrxxappo matching for instance (A, B) in time (n for sufficiently small constants E >> 0. We begin by presenting variant of this algorithm for the static setting, referred to as STATIC- MATTC--MTCCINEG AACHINGG( A, B). The next section shows that our variant can be extended to a bottom-up algorithm for the fully dynamic setting achieving sub-linear update time.
(1) The root F consists of a bounding cell with side-length D, containing all points in X, and its subdivision into p² smaller subcells. The root F has p² children, de- noted by C,, corresponding to its subcells, each of side-length D/p. The root is at level of
modified
(2) Each internal node at level > 0 corresponds to one of the subcells of its parent node. Thus, an internal node v at level i consists of cell of side length D/p", and its subdivision into p² smaller subcells, each of size D/p x D/p
For the input sets A, B with bounded spread U we first apply a random shift to all the points by a fixed random vector r E [T]2 to get A' = {a + r : a E A} and B' = (b+r:be E B). From now on, we are only inter- ested in sets A' and B'. Slightly abusing notation, we will simply refer to these sets as A and B. Our first step is to construct a restricted p-tree T on the point set X = A U B. Using T in a bottom-up manner, we obtain matching be- tween the point set A and B as follows. Intuitively, as the number of unmatched red and blue points belonging to some node might not be balanced, we greedily match as many red-blue pairs of points as possible that node. Handling the unmatched points is deferred to the parent node, which matches many excess points from its children as possible, and so forth. Therefore, for any node of T, we keep track of the set of red points (.(A) and blue points Y- that we match at that node, and the set of unmatched monochromatic points that are in excess, denoted by Ey.
For node of the underlying p-tree T,, let X2 be the set of points from X belonging to the bounding cell of v. We say that a point x E X belongs to node v if & E X,, Now, notice that for some nodes v (on both low and high levels) the set X,, might be empty. To avoid storing unnecessary empty cells (which would result in slower data structure), define the restricted p-tree shortly p-tree, which is subtree of the underlying p-tree, where the root and leaves are defined as follows:
(1) The root r is a node of the underlying p-tree such that X, = X, but which does not have a child with X, X. In other words. is the node with X, = X located on the lowest possible level of the underlying p-tree.
An essential element of our algorithm is the choice of the matching algorithm used throughout the traversal of T. Ob- serve that the excess points forwarded to an internal node of T by its children could be of size ((n). To overcome this problem, our algorithm will rely on a data structure representing excess points forwarded by node of T by the center point of the cell of Using this data structure we are able to efficiently compute (and maintain) an implicit match- ing of the unmatched points at all a E V. which we can efficiently turn into an explicit matching of the underlying points when queried.
(2) A leaf of the p-tree is cell such that |Xo| S p², but whose parent u has [xul >> p².
As we insert or delete points in the dynamic setting, notice that the restricted p-tree might change over time: the root, as well as the leaves might change to higher or lower levels of the underlying p-tree. However, as we assume that the spread of the pointset X is bounded by U' at all times, we know that the height of the restricted p-tree is bounded by = log, U and is therefore constant if p is chosen to be sublinear in n.
It can be shown that a restricted p-tree w.r.t. set X, |X| n,
3.1. Sub-instances via transportation solver
can be maintained under insertions and deletions of points in Opp + log, U) amortized update time (see c.g. dde Berg et al. 2007
An implicit matching is an optimal assignment >: T associated with vertex of of an Euclidean transporta- tion problem defined by some inputs (Y,(4), Y,(B), d).
z-ordering The z-ordering of the point-set X w.n.t an underlying p-tree T is defined via specific ordering of chil-
To make sense of the above vague definition, we de-
4

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
scribe procedures IMPLICIT-MATCHING and EXPLICIT- MATCHING that will act as sub-routines in both static and dynamic algorithm for Euclidean matchings. The pur- pose of these functions is to obtain two representations of a maximal matching between specific points forwarded to an internal vertex from its children in a restricted p-tree.
EXPLiCIT-MATCHIN procedure The inputs of the pro- cedure are sets (A) and (B)) as well implicit match- ing [+,,,] between those points at an internal node Its output is matching between (AA) and Y, (B).
The procedure naturally translates the assignment >; [Y,] to perfect matching j4,[Y,, between the points (A)) and Y, (B)) follows: match >g [a,] (') many different pairs of points of the form (a, b). where a and b lies in the sub-cell whose centers are a' and 6' respectively. We refer to the latter matching as an explicit matching of node v w.r.t. its corresponding point set Y,. For an example involving both implicit and aggregated sub-instances, see Figure |
Given subset of points belonging to node in the p-tree, in the IMPLICIT-M WWW.ICIT-MATCHINGG we first define the notion of a ag- gregated sub-instance for that node. To aggregate the points of a given cell, we move all of them to the cell's 's center. We define an Euclidean transportation problem on this aggre- gated sub-instance, and the resulting optimal assignment will correspond an implicit matching of v. EXPLICIT- MATCHING converts this implicit matching into matching of the original subset of points, which we refer to as the explicit matching.
In the following lemmas, we summarize the above algo- rithms and their running times.
Lemma 3.1. Let be an internal node of p-tree T X == A U B with |A] |B| n. Let E {Eu : u.uEC be the pointers to some monochromatic sets of v's children forwarded to e, as well as their car- dinalities {|Eul : u E Cy.. There is an algorithm E, {(El : E C.)) that chooses balanced red and blue sets (A) and (B)),, computes an implicit (perfect) matching %.. [Y,)) between them, in time O(p² log p log Additionally, the function returns the (pointer to the) monochromatic set of unmatched points E, at node v.
IMPLICIT-MATCHING procedure The input of the pro- cedure is vertex u of p-tree T, the monochromatic sets of points {Eu : u E C.} (denote their union by E, and partition into blue Eg and red EA points) forwarded to v from its children and their cardinalities f|Eul u E C.). Without loss of generality, assume that E 2 k- = Exl.
Its output consists of an optimal assignment of transporta- tion problem (i.e., an implicit matching) defined by k points of Eg and all points of ΣA (i.e., the points we can match within the cell of v) and | - h blue points, which will be forwarded to the parent of e in T.
A crucial property of the Lemma 3.1 is that it has a running time only polylogarithmic with respect to n although E could consist of (2()) points. This can be naturally achieved through an implicit representation of the point sets ( (Fn, u E C.). We formalize this process in Appendix B
Choosing the points to match. Let (B) stand for the first k points of εn with respect to some ordering (say the z-ordering). Let E,- Eg (B) be the excess set, i.c. the set of points not matched at v, (A) = EA and Y.- = (A) U (B). Excess set E, will be forwarded to the parent of v.
Now, as described above, it is easy to convert an implicit matching to a matching of the input points.
Lemma 3.2. Let T be p-tree X AUA B and d its depth. The following statements hold:
Given a vertex v of p-tree T, and an implicit matching no [Y.] w.r.t. v and set Y,, the corresponding explicit matching (Y() can be obtained in time (((y))d).
Computing the implicit matching. We define an instance of Euclidean transportation problem over the center points of all p² sub-cells of the cell of c. Note that every such sub-cell E C, contains monochromatic subset of E. Define A and B to be center points of sub-cells containing only red and blue points of Y, respectively. We set &a' for a' A and -dy for b' E B to correspond to the number of red and blue points of Y, of the cells with centers a' and R, respectively.
2. Given tree T where each leaf stores a matching and each non-leaf node an implicit matching, the corre- sponding perfect matching on X can be reported in O(n) time.
3.2. Static Algorithm
Next, we present the formal description of the bottom-up algorithm. To this end, we define the subroutine MATCU(*) as follows. If node is leaf of the restricted p-tree:
The optimal assignment 7: [Y,] obtained by solving the Euclidean transport problem defined by (As, B", s, d) using the algorithm of Theorem 2.1 will be returned as the implicit matching of vertex V. Crucially, observe that setting up the Eucledian transport instance doesn't 't require the knowledge of the exact location of vertices in (A), (B) just their cardinalities in the sub-cells of v.
will
1. Determine sets (A),), (B(B) and E.: Let (A) = X, A be the red points, (B) = Xx, n B the blue points, and 7/- = [Y(( - |Y the difference in cardinality, at node v.
5

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
x
X
x
x
x
X
X
x -3
Figure 1. Examples of instance (left), implicit matching of the aggregated sub-instance (middle), and the resulting matching sub- instance (right). In the left picture, red points are represented by circles, while blue points are represented by crosses. The middle picture shows the implicit matching on aggregated sub-instance, with numbers next to vertices representing their demands and supplies, while numbers next to edges represent the weighting of the edges. The right picture shows the corresponding matching of the input points.
If ?o O, we set E,, V. If T7 > O, assign subset of Th red points chosen the last Th points in the z-ordering from (A)) to Ev,, and set (A) (A(A) \ E.p and Y, (B) = (B) (in the same way as presented in Section 3.T. Analogously, if no < 0, we assign a subset of T/. blue points from (((A) to E,,, and set (B) We refer to E, as the excess points at v.
are shifted randomly by the vector r. the probability of an edge being cut by a grid in a node of the p-tree is propor- tional to the length of the edge and inversely proportional to the spacing of the grid.
However, the algorithm of (Agarwal and Varadarajan 2004) works in top-to- to-tottooommmoooooooooooooooottt approach (building the matching starting from the root and progressing to the leaves) in con- trast to our algorithm, which takes a bottom-up approach. The bottom-up approach allows us to guarantee that for any node E the number of nodes not matched within the cell of simply corresponds to the difference in the number of blue and red nodes inside the cell of v. Intuitively, this im- plies that points are matched near-by neighbor whenever possible.
(B) Ey and (A()) (A).).
2. Compute matching: Obtain a matching (vv between red points (A) and blue points (B) using any efficient exact algorithm (e.g. the Hungarian algorithm Munnrres| 1957)), and compute its cost c(py))
3. Return the cost and the excess set (c(p)), (((p,),E.) E.).
Otherwise, node is internal:
This results in an arguably simpler analysis differing from that of (Agarwal and Varadarajan 2004), and crucially rely- ing the concept of implicit matchings. We provide self- contained analysis of the static algorithm in Appendix [
1. Compute implicit matching: For each child u of v, let Ey be the set of excess points returned Invoke algorithm IMPLICIT- E C.)) from Section 2 (see Lemma 3.1 Let Ye Yy] ] be the resulting implicit matching w.r.t. node and set Y, computed via the algorithm, and E, the set returned by the algorithm. Set the estimate of the cost for the matching in to c(Hu) = ((C)[(r)) + Σuec, c(pa).
by (u)C
MATCHING &, EuEal=
Theorem 3.3 (Static). For any Euclidean matching instance (A, with |A| [B] = n, and spread U S nn,, for some constant e > O, and for any constant E (0, ) algorithu STATIC-MATCHINGE MATCHING(A, B) computes an O(1/E)- approximate Euclidean matching with high probability in On 111) time.
2. Return the cost and excess set (c(pw), Ex).
Our algorithm MANTIC-MTCCHING(/ A, B) first invokes the MATCH(() subroutine, where r the root of T on A B. Note that this yields implicit matchings 70 [Y,] at internal nodes T. Using Lemma 3.2 Part (2), we can transform the implicit matchings > [Y.] of the aggregated versions of Y, into matchings of the original points μ' [Y,] of roughly the same cost (see Claim C.4] for a formal bound on the slack introduced by the aggregation). Finally, we return a maximal matching H4, on the entire point set X, along with its already computed cost c.
4. Dynamic Algorithm
In this paper, we present two dynamic data structures for maintaining perfect matchings under insertions and dele- tions of pairs of red and blue points. Both data structures naturally convert our static, bottom-up algorithm from the previous section into a dynamic one.
Our first algorithm EUCLMATCH1 implicitly maintains matching for each node of the tree. This allows to (i) estimate the cost of the entire matching after each update and (ii) report the current approximate matching upon a query in time proportional to its size. Our second algorithm EUCLMATTC2 explicitly maintains the output matching at
Our analysis of the algorithm builds on the following key idea from (Agarwal and Varadarajan 2004): when all points
6

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
all times. We defer the second algorithm Appendix due to page limitations.
nodes of T, that is nodes which either the implicit match- ing the excess set changes, lie along a path from a leaf to the root of T', Hence, the most costly operation of the algorithm, the updating of the implicit matchings, only has to be repeated at a small number of nodes.
Initialization For any given set of points A and B, we will next describe how to initialize a matching by essen- tially running the bottom-up algorithm SATTINGMATCIING on the instance (A, B). Formally, the INITIALIZE(A, (A B) initializes the restricted p-tree T w.r.t. X := A U B, and the following corresponding values: the set of unmatched points Y. (B), the excess set E,, each node of T, the implicit matching > [Y.] and the cost c(fo.) of the matching w.r.t. the subtree rooted at each internal node, and the matching He and its cost c(14p) at each leaf. This process
The conrectness of the algorithm follows from the algorithm always satisfying that after handling an update its output corresponds to the run of static algorithm on the current state of the input, up until the different choices the under- lying optimal bi-chromatic matching algorithm (Hungarian a algorithm) and optimal geometric transport algorithm (see Theorem 2.1 makes. Note that the choice of underlying ex- act algorithms is arbitrary long as they run in polynomial time.
takes O(n time by Theorem 3.3
Handling updates As the procedures for insertions and delctions are similar (see Appendix D) we focus on pre- senting the insertion procedure. Without loss of generality assume that we insert a red point a to A. The algorithm (N)ERT(a) (see Algorithm 2 proceeds in two steps: updat- ing the restricted p-tree T', and updating the matchings at each affected node, i.e. node u where a E X
To report the changes to the matching, the algo- rithm converts the implicit matchings at the affected nodes to an explicit matching of X via Lemma 3.2 To query the cost of the approximate solution the algorithm simply returns the approximate cost value c(u,) that is maintained at root of T.
Queries
First, updating T follows standard upproach, similar to the method described in e.g. (Har-Peled (01T)) Namely, one or both of the following cases can happen:
The following theorem summarizes the results of our sim- pler dynamic algorithm.
1. After insertion of a to A, a leaf of T might contain more than p² points. In this case, now becomes an internal node of T,, while new subtree has be built under u. We refer to these nodes as new subroots, while all other nodes belonging to this subtree are called marginal nodes. Finally, all nodes in T that are not marginal are refered to as non-marginal nodes.
Theorem 4.1. Let & E (0, and (A, two point sets in R² with |A| [B] = n, and spread U < n² for some con- stant > 0. There is dynamic data structure that supports insertions and deletions in O -E")) worst-case update match-
time and maintains an expected O(1 (11/))-aaaa /e/e)e-pprooiin ing between the dynamic set (A, B). Initialization takes O(n2+)) time. Reporting the matching, the changes to the matching when an insertion or deletion occurs (recourse), and the cost of the matching, take ((n), (n) and O(1) time, respectively.
2. If does not belong to the current root of the restricted p-tree, a internal node of the underlying p-tree that is on higher level than the previous root now becomes the new root T.
Second, we define the recursive procedure INSERT-UPDATE(C, a) which updates the matchings at each affected node of T in a bottom-up manner (see a)
Improved dynamic algorithm Ideally, the algorithm would maintain the solution matching explicitly. This does require however, that the explicit matchings implicitly repre- sented at all nodes of T are explicitly maintained at all times. This requires careful propagation of changes throughout the p-tree and the ability to update explicit matching effi- ciently when the underlying implicit matching changes. To dynamize the maintenance of explicit matchings, we rely on the fact that a single insertion might only change the excess set stored at each node by point.
Alorithm D. Specifically, we define INSERT-UPDATER for three types of affected nodes of T: non-marginal internal nodes, non-marginal leaves, and new subroots. We focus on presenting the case for non-margine marginal internal nodes, as the other two cases are simpler, and can be found in Appendix D
If a node e is a non-marginal internal node, the procedure INSERT-UPDATE(O, a) first updates set E., for the affected child of v, via INSERT-UPDATE(U.OO a). Then, we up- date the excess set F., accordingly. Finally, we invoke the IMPLICIT- MATCHING(.) procedure to obtain an implicit matching for the remaining points Y, at v.
This implies that the difference between the updated and pre-updated matching lies along a short alternating path of old and new matching edges. We find this augmenting path using a graph shortest- path based sub-routine, which, due to page limitations, we defer to Appendix E
Crucially, observe that throughout this process the affected
7

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
5. Experimenntl Evaluation
Speecup: Time State: (yyam Approratiion Goiss
1/003
With this work, we provide an opensource, single-threaded C++ implementation of our algorithm in Appendix E To demonstrate the practicality and effectiveness our pro- posed algorithms, we conducted experiments on 2.2 GHz Ubuntu 22. 04. 4 system with AMD Opteron 6174, us- ing synthetic and real world data sets. For reproducibil- ity and comparability, we report on measurements for the same 2D distributions and data sources used in the recent work (Gattani et al. 2023). For the synthetic data, we drew the point coordinates from the uniform distribution on in- tegers between and 500 to obtain the Uni form datasets. For the Gauss i an datasets, we drew the point coordinates from the normal distribution (μ = 0.5, α = 0.25) prior to scaling by a factor of 500 and rounding to integers. For the real data, extracted all trips from the Yellow-Cap Taxi datase December 2009, and, following (Gattani al. 2023), extracted those having at least 3 minutes duration, at most 110 mph speed, and latitude/longitde in the range of [ 74.5, - 73.5) x [40, 41], yielding 14 324 017 trips with one 'Pickup" and 'Dropoff' location each.
1000
1000
2000
3000
1000
5000
80000
1000
8000
Not. samples
Spocdus:: State Dyynaiic App coxmmations Piot.p Droport
100
10
5.1. Experimental Results
2000
Feot suurphee
Speedup: Dynamic Static Approximation Algorithms To measure the speedup of dynamic over baseline approx- imations, we measured the update time and time to compute a approximate matching from scratch, for various levels of approximation quality. Figure shows the speedup factors (p 8) for inserting S 000 samples from Uni rm vs Guassian (top) and P Pickup vs Dropoff (bottom).
Figure 2. Speedup of dynamic algorithm over static for approxima- tions with p 8.
Distance of Un torm and Uniior
150
330
H27
200
50
Ansertion-Worloadee Estimating Wasserstein distance using Exact vs Approximate Matchings To demon- strate the usability of approximations of minimum cost bi- chromatic matchings to estimate the Wasserstein distance, we estimated the empirical Wasserstein distance using exact and approximate matchings on identical and different dis- tributions for various levels of approximations and sample sizes. Figure 3 shows the empirical 1- Wasserstein distance of two Un Eo rm distributions (top) using exact minimum cost matchings and approximations with p = 2, p 8, and p = 32, for up to 200 samples. The bottom figure shows the same estimates for Uniform vs Gaussian
1.5
125
250
512
1004
Foot sampkes
300
70
Enserton-Wooklladee Convergence and Update Time for Very Large Sample Sizes To demonstrate the scalability of our approach, we measured the convergence and update time of our dynamic approximations for sample sizes up to 1 000 000, rendering the exact solvers impractical. That is, running the code of (Xu and Ding 2024b) on Uni form
40
310
10
20
20
sumpphes
https: / //// httpg//thub..cm///eeeeeeeee www.o/////////o//ooooocooo / mattc
Figure 3. Estimating the 1- 1-Wasserstein distance using exact and approximate matchings.
See https:: htwp:///wwwww.///. http://cwwwwaa/eec//////////
tle-trip- record- datcordddata.pag...
8

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
Gauss ian slowed down substantially with increasing number of samples. Inserting 200 samples took more than 12GB RAM and longer than 3 hours in total, which is when we had to terminate the process as it was running out of memory. (Note that this already > 5s per update on a very small instance.) In contrast, Figure 3 (top) shows the convergence of empirical I-Wasserstein distance, us- ing p 2, 4, 8 for approximations, for the Unife mirform On form benchmark (bottom three curves) and the Fnrm vs Gaussian benchmark (top three curves), The bottom figure shows the update time, in milliseconds, of our dynamic algorithm on both benchmark data sets, having sample sizes up to 000 000.
180
140
(4
Emprical WassersonDinnn
120
100
100
1000
30000
100000
Ax30P
Nof. sa sampis
Upoae::: (76)
Fully-Dynand Workloads: Drift of Pickup vs Dropoff distributions over time To demonstrate the effectiveness of our dynamic algorithms to monitor the drift of time- dependent, spatial distributions, we sorted the trips in the Taxi dataset by pickup time and use a sliding window, hav- ing a width of 10 samples, to obtain update sequences that contain insertions and deletions. (tg)) shows the empirical W-Wasserstein distance between the Pickup and Dropoff locations, using approximations with p 4 and 16. The bottom figure shows the update time of our fully dynamic algorithm on the Taxi benchmark.
01
0.01
Nof. samples
Figure Convergence and update time on large datasets.
Liatriuut anei Lindi ticeup Drspont
490001
5.2. Discussion of experimental results
The experimental results show that the dynamic algorithm is orders of magnitudes faster than computing static approxi- mations for minimum cost bi-chromatic matching (Figure In practice, the approximate matching costs are within very small factor (i.c. < 2) of the minimum bi-chromatic matching cost, which enables very accurate and effectiv estimation the 1-Wasserstein distance (Figure and top part of Figure 4. As suggested by our Theorem tradeoff between approximation quality and update time (bottom part of Figure 4). Updates of the dynamic algorithm took typically between one and ten millisecond and showed a clear separation, between inputs with zero and non-zero distance, at sample sizes around 10 000 samples. The time dependent spatial distributions of Pickup and Dropoff lo- cations in the Taxi dataset show a drift in the empirical 1-Wasserstein distance (Figure top), and the dynamic al- gorithm allows for effective monitoring of the change large window sizes on real data sets (Figure 5 bottom).
30000
Empreil Nioooret 30d6
30000
20000
2,000
we observe
2000100
400000
BUOOOD
Posilio of Slding Window (LOk samp bs WACk))
Updane Texi
30
Updes Time
831
Impact statement
0.001
230000
600020
#02000
Pf Slading Window (1Ok corppss wAX))
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences our work, none which we feel must be specifically highlighted here.
Figure 5. Drift in the Tickup-rropoff distributions of the Taxi datasct (top) and update time (bottom).

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
References
unit-capacity graphs. In 2020 IEEE 61 Annual Sympo- sium on Foundations of Computer Science (FOCS), IEEE, 93-104.
Pankaj K. Agarwal, Hsien-Chill Chang, Sharath Raghven- dra, and Allen Xiao. 2022a. Deterministic, near- near-linear e-approximation algorithm for geometrie bipartite match- ing. In STOC '22: 54th Annual ACM SIGACT Sympo- sium on Theory of Computing, Rome, Italy, June 20- 24, 2022, Stefano Leonardi and Anupam Gupta (Eds.). ACM, 1052-1065. https://doi org/10. 11457 3519935. 3519977
Amir Azarmehr, Soheil Behnezhad, and Mohammad Roghani. 2024. Fully Dynamic Matching:
Approximation in Polylog Update Time. In Proceed- ings of the 2024 ACM-SIAM Symposium on Discrete Algorithms. SODA 2024, Alexandria, VA, USA, Jan- uary 7-10, 2024, David P. Woodruff (Ed.), SIAM, dttpps://don ottp/..rgg11... 1137/1
Pankaj K. Agarwal, Sharath Raghvendra, Pouyan Shirza- dian, and Rachita Sowle. 2022b. An Improved E-Approximation Algorithm for Geometric Bipartite Matching. In 18th Scandinavian Symposium and Work- shops on Algorithm Theory, SWAT 2022, June 27- 29, 2022, Tórshavn, Faroe Islands (LIPlcs, Vol. 227). Schloss Dagstuhl- Lcibniz-Zentrum für Informatik, 6:1-6:20. /tttps://do.. otg//0..31100011E..... SWAT. 2022.66
3040 3061.
978 778161197797 109
Yogesh Balaji, Rama Chellappa, and Soheil Feizi 2020. Robust Optimal Transport with Applications in Genera- tive Modeling and Domain Adaptation. In Advances in Neural Information Processing Systems 33: Annual Con- ference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual.
Richard Bellman. 1958. On a routing problem. Quarterly of applied mathematics 16, (1958), 87-90.
Pankaj K. Agarwal and Kasturi R. Varadarajan. 2004. A near-linear constant-factor approximation for euclidean bipartite matching?. In Proceedings of the 20th ACM Symposium on Computational Geometry, Brooklyn, New York, USA. June 8-11, 2004, Jack Snoeyink and Jean- Daniel Boissonnat (Eds.). ACM, 247-252. https:// doi ori.org110..m 1 997797777777777 997856
Gaspard Beugnot, Aude Genevay, Kristjan H. Gree- newald, and Justin Solomon. 2021. Improving ap- proximate optimal transport distances using quanti- zation. In Proceedings of the Thirty-Seventh Confer- ence on Uncertainty in Artificial Intelligence, UAI 2021, Virtual Event, 27-30 July 2021 (Proceedings of Machine Learning Research, Vol. 161). AUAI Press, 290-300. https://proceedinng mlr /ress/ he1/beugnot22.. html
Jason Altschuler, Jonathan Niles- Weed, and Philippe Rigol- let. 2017. Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration. Advances in neural information processing systems 30 (2017).
David Alvarez-Melis and Nicolò Fusi. 2020. Geometric Dataset Distances via Optimal Transport In Advances in Neural Information Processing Systems 33: Anmual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, Hugo Larochelle, Marc' Aurelio Ranzato, Raia Hadsell, Maria- Florina Balcan, and Hsuan-Tien Lin (Eds.).
Sergio Cabello, Panos Giannopoulos, Christian Knauer, and Günter Rote. 2005. Matching Point Sets with Respect to the Earth Mover's Distance. In Algorithmss ESA 2005, 13th Annual European Symposium, Palma de Mallorea. Spain, October 3-6, 2005, Proceedings (Lecture Notes in Computer Science Vol. 3669). Springer, 520 5311. :ttp:///00i.0r000.00000111000011
Alexandr Andoni, Aleksandar Nikolov, Krzysztof Onak, and Grigory Yaroslavtsev. 2014. Parallel algorithms for geometric graph problems. In Symposium on The- ory of Computing, STOC 2014, New York, NY, USA, May 31- June 03, 2014, David B. Shmoys (Ed.). ACM, 574-583.
Sergio Cabello, Panos Giannopoulos, Christian Knauer, and Günter Rote. 2008, Matching point sets with respect to the Earth Mover's Distance. Comput. Geom. 39, 2 (2008), 118-133 https://doo http://oo...o/1..111111//// COMGEO. 200.0..2000.. 0011
http:/00.....//0..0...... 1145/
Jiezhang Cao, Langyuan Mo, Yifan Zhang, Kui Jia, Chun- hua Shen, and Mingkui Tan. 2019. Multi- marginal Wasserstein GAN. In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, De- cember 8-14, 2019, Vancouver. BC, Canada. 1774-1784.
David S Atkinson and Pravin M. Vaidya. 1995. Using geometry solve the transportation problem in the plane. Algorithmica 13, 5 (1995), 442-461. htttps://doi org/10. 1 007 /F01 00089900
Kyriakos Axiotis, Aleksander Madry, and Adrian Vladu. 2020. Circulation control for faster minimum cost flow in
Kevin C. Cheng, Shuchin Aeron, Michael C. Hughes, and Eric L.. Miller. 2021. Dynamical Wasserstein Barycenters
10

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
for Time-series Modeling. In Advances in Neural Infor- mation Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurlPS 2021, December 6-14, 2021, virtual, Marc' Aurelio Ran- zato, Alina Beygelzimer, Yann N. Dauphin, Percy Liang, and Jennifer Wortman Vaughan (Eds.). 27991-28003.
Monika Henzinger, Sebastian Krinninger, Danupon Nanongkai, and Thatchaphol Saranurak. 2015. Unify- ing and Strengthening Hardness for Dynamic Problems via the Online Matrix-Vector Multiplicatiion Conjecture. In Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC 2015, Port- land, OR, USA, June 14-17, 2015, Rocco A. Servedio and Ronitt Rubinfeld (Eds.). ACM, 21-30. https
Søren Dahlgaard. 2016. On the Hardness of Partially Dy- namic Graph Problems and Connections to Diameter. In 43rd International Colloquium on Automata, Lan- guages, and Programming. ICALP 2016, July 11 11-15, 2016, Rome, Italy (LIPles, Vol. 55), loannis Chatzigian- nakis, Michael Mitzenmacher, Yuval Rabani, and Davide Sangiorgi (Eds.). Schloss Dagstuhl- - Leibniz-Zentruum für Informatik, 18:1144:14 42.0/LIPP PCS. TCALP 2016.44 48
d/doi. 1111111452755311 2746609
Piotr Indyk. 2007. A near linear time constant factor ap- proximation for Euclidean bichromatic matching (cost), In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms. SODA 2007, New Orleans, Louisiana, USA, January 7-9, 2007. SIAM, 39- 42. 1283383.1283387
Mark de Berg, Herman Haverkort, Shripad Thite, and Laura Toma. 2007. 1/O-efficient map overlay and point location in low-density subdivisions. In International Symposium on Algorithms and Computation. Springer, 500-511.
Hicham Janati, Thomas Bazeille, Bertrand Thirion, Marco Cuturi, and Alexandre Gramfort. 2019. Group Level MEG/EEG Source Imaging via Optimal Transport: Mini- mum Wasserstein Estimates. In Information Processing in Medical Imaging 26th International Conference, IPMI 2019, Hong Kong, China, June 2-7, 2019, Proceedings (Lecture Notes in Computer Science, Vol. 1 1442),, Albert C. S. Chung. James C. Gee, Paul A. Yushkevich, and Siqi https: /ttp::://doi.ooo org/ 10. 10.1009777--11 3- 20351-203311122 1_
Harold N Gabow and Robert E Tarjan. 1989. Faster scaling algorithms for network problems. SIAM J. Comput. 18, 5 (1989), 1013-1036.
Alfred Galichon. 2018. Optimal transport methods in eco- nomics. Princeton University Press.
Bao (Eds.). Springer, 743-754.
Akshaykumar Gattani, Sharath Raghvendra, and Pouyan Shirzadian. 2023. A Robust Exact Algorithm for the Euclidean Bipartite Matching Problem. In Advances in Neural Information Processing Systems 36: Annual Con- ference on Neural Information Processing Systems 2023, NewrlPS 2023, New Orleans, LA, USA, December 10- - 16, 2023.
Harold W Kuhn. 1955. The Hungarian method for the assignment problem. Naval research logistics quarterly 2, 1-2 (1955), 83-97.
Huidong Liu, Xianfeng Gu, and Dimitris Samaras. 2018. A Two-Step Computation of the Exact GAN Wasser- stein Distance. In Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stock- holmsmässan, Stockholm, Sweden, July 10-15, 2018 (Pro- ceedings of Machine Learning Research, Vol. 80). PMLR, http://prooceedingg... wwl.a.gsmmlrrr...mm..mmm.....mmm www.prsssmmprpr..s//mmm
Alexandre Gramfort, Gabriel Peyré, and Marco Cuturi. 2015, Fast Optimal Transport Averaging of Neuroimag- ing Data. In Information Processing in Medical Imag- ing- 24th International Conference, IPMI 2015, Sab- hal Mor Ostaig, Isle of Skye, UK, June 28- July 3, 2015, Proceedings (Lecture Notes in Computer Science, Vol. 9123), Sébastien Ourselin, Daniel C. Alexander, Carl-Fredrik Westin, and Manuel Jorge Cardoso (Eds.). Springer, 261-272. https: http:////ooooog/..../o 1007/ 978-3-319-19992-2 20
v 0 1//liull0 8d. html
Yongchao Liu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. 2017. Primal-dual hybrid gradient method for distributionally robust optimization problems. Oper Res.
Lett. 45, 6 (2017), 625 https: /tte:://ooo..g////
1016/5. 2017 10.001
James Munkres. 1957. Algorithms for the assignment and transportation problems. Journal of the society for indus- trial and applied mathematics 5, (1957), 32-38.
Joachim Gudmundsson, Oliver Klein, Christian Knauer, and Michiel Smid. 2007. Small Manhattan networks and al- gorithmic applications for the earth mover's distance. In Proceedings of the 23rd European Workshop on Compu- tational Geometry. Citeseer, 174-177.
Victor M. Panaretos and Yoav Zemel. 2019. Statisti- cal Aspects of Wasserstein Distances. Annual Review of Statistics and Its Application 6, Volume 6, 2019
Sariel Har-Peled. 2011. Geometric approximation algo- rithms. Number 173. American Mathematical Soc.
(2019), 405-431.
annurev-state samurev-statissticc-0070771

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
Sharath Raghvendra and Pankaj K. Agarwal. 2020. Near- linear Time e-Approximation Algorithm for Geometric Bipartite Matching, J. ACM 67, 3 (2020), 18:1-18:19. http:////bi..gg1....111...
Xiaoyang Xu and Hu Ding. 2024b. A Novel Skip Or- thogonal List for Dynamic Optimal Transport Prob- lem. In Thirty-Eighth AAAI Conference Artificial Intelligence, AAAI 2024, Thirty-Sixth Conference on Innovative Applications of Artificial Intelligence, IAAI 2024, Fourteenth Symposium on Educational Advances
Yossi Rubner, Carlo Tomasi, and Leonidas J. Guibas. 2000. The Earth Mover's Distance Metric
Artificial Intelligence. EAAI 2014, February 20-27, 2024, Vancouver, Canada, Michael Wooldridge, Jen-
for Image Retrieval. Inl. J. Comput. Vis. 40, 2 (2000), 99-121. pettps://dok. 1otg/10.1000000 1026543900054
nifer G. Dy, and Sriraam Natarajan (Eds.). AAAI Press, nttp:://co1.orr/10..600
BAAT V38T18 18. 30073
R. Sharathkumar. 2013. A sub-quadratic algorithm for bi- partite matching of planar points with bounded integer coordinates. In Symposium on Computational Geome- try 2013, SoCG '13, Rio de Janeiro, Brazil, June 17-20, 2013. ACM, 9-16. attps://dbi.co 2462356. 248024
1145/
R. Sharathkumar and Pankaj K. Agarwal. 2012a. Al- gorithms for the transportation problem in geomet- ric settings. In Proceedings of the Twenty-Third An- mual ACM-SIAM Symposium on Discrete Algorithms, SODA 2012, Kyoto, Japan, January 17-19, 2012. http:///oooorr//1..
SIAM, 306-317. 1. 1.9786619970099..
R. Sharathkumar and Pankaj K. Agarwal. 2012b. A near- linear time E-approximation algorithm for geometric bi- partite matching. In Proceedings of the 44th Sympo- sium Theory of Computing Conference, STOC 2012, New York, NY, USA, May 19 22, 2012. ACM, 385- 394.
221 4014014
Ilya O. Tolstikhin, Olivier Bousquet, Sylvain Gelly, and Bernhard Schölkopf. 2018. Wasserstein Auto-Encoders. In 6th International Conference on Learning Represen- tations, ICLR 2018, Vancouver, BC, Canada, April 30 May 3, 2018, Conference Track Proceedings. Open- Review.net. htttp:://peerree
id H 7 0b
Jan van den Brand, Yin Tat Lee, Danupon Nanongkai, Richard Peng, Thatchaphol Saranurak, Aaron Sidford, Zhao Song, and Di Wang. 2020. Bipartite Matching in Nearly-linear Time on Moderately Dense Graphs. In 61 IEEE Annual Symposium on Foundations of Computer Science, FOCS 2020, Durham, NC, USA, November 16- 19, 2020. IEEE, 1109/FOC546 700 2020 00090
https: /ttp:://///i..r/11100
Xiaoyang Xu and Hu Ding. 2024a. A Novel Skip Orthog- onal List for Dynamic Optimal Transport Problem. In Proceedings of the AAAI Conference on Artificial Intelli- gence, Vol. 38. 20838 20846.
12

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
A. Related Works
allows computing an optimal Euclidean bi-chromatic matching in (nn") time. {Atkinson and Vaidya 1995) presented a Primal-Dual based scaling algorithm for points from R that, supporting points with multiplicity at most W, runs in On log(n) log(W)) time. Near-quadratic algorithms for computing an optimal matching (of points with multiplicity) can be derived from recent advances for minimum cost max-flows in graphs (van den Brand et al. 2020).
Static Algorithms: The classical Hungarian algorithm (Kuhn 1955)
Sub-quadratic algorithms are known for the following two important cases. In case d = 2 and all coordinates are integers <U, Sharathkumar 2013) presented an 0(n3/2+ log(nl))), for any fixed 5 E (0, 1]. In case both AA B are i..dd samples (from two distributions), (Gattani al. p023) recently presented an O(n2= za log log(U) (n)), where Φ(n) = poly (log m)) for d 2, expected time algorithm, and showed its practicality for estimating the 1- 1-Wasserstein distance of distributions on IR? for r v 50 000 samples in approximately 700 seconds.
Static Approximations: (Agarwal and Varadarajan 2004) presented a top-down algorithm that computes, with with high probability,, an Ollo !) >approximatioon in O +C) >-time,, for any fixed E E (0. 1]. (Sharathkumar and Agarwal, 2012b). (Raghvendra and Agarwal 2020) obtained (n poly(log =)) time algorithm that, using augmenting paths in randomly shifted quad-tree, computes with high probability a (1 + -apppoxxime matching. The poly n, nn,-- )-accoo
was further improved in (Agarwal et al.) 2022b). Recently, (Agarwal et al. 2022a) presented the first deterministic (1 (1+ -aprroxipioiinianinnn algorithm with near-linear n log )o(a) time.
Variants: Using spanners and fast minimum cost max-flow solvers, is possible to compute a (1+ (apprrooimiiait of
the cost-value in (12/22 or in ((12/3++(11) time, also supporting point multiplicity (see (Gabow and Tarjan 1989) and
(Axiotis et al. 2020)). See also e.g.. (Cabello et al. 2005 2008) and (Gudmundsson et al. 2007). (Indyk] 2007) showed that importance sampling can be used for obtaining. with high probability, an 0(1) (apppoxxxinitoo of the cost-value in (n pooy(Iogn)) n²) time, which however also does not allow report matching. (Altschuler et al.. 2017) presented an algorithm that also reports matching, but with an additive error e in time Ó(n²/ /ea). For the parallel and streaming model, (Andoni et al. 2014) presented an n 2,,,,11 (1) time algorithm for R² that allows to compute an appprxxxxxiiaaooonn for the cost-value of an optimal matching. (Sharathkumar and Agarwal 2012 also presented a dynamic approximation algorithm
that maintains a partial-matching
The more general problem of dynamic optimal transport has been considered in recent work (Xu and Ding 2024b). Specifically, the goal is to design data structure that accelerates executing one iteration of the Network Simplex Algorithm on an explicit complete bi-partite graph. Their work is concerned with maintaining an exact solution, within numeric computing accuracy, after the insertion/deletion of a vertex to the min-cost flow problem. Since their algorithm operates on a complete, bipartite graph, it requires I(n²) space to store the edge weights and M() update time to insert a new vertex. Moreover, in the worst case, the update time can take up to a quadratic time, As such, the algorithm and implementation cannot handle very large instances in practice.
B. Formal description of the implicit and explicit matching sub-routines
IMPLICIT MATCHING Recall that the input of the procedure is vertex e of p-tree T and the monochromatic set of points & := {Ea u } forwarded to u from its children and their cardinalities [|Eul E C.).
Choosing the points to match at v. Formally, first divide E into collections of red and blue monochromatic sets E := {Eu. Cv, F., is red} and & == {F., u E C. F,, is blue), respectively.
Then, determine the total number of red and blue points in E and Eg, respectively, by summing up the provided set cardinalities in each collection. Suppose that there are kj red points in E, k2 blue points in Eg, and that without loss of generality hn S kg. In this case, we aim to match all red points from sets E To pick k points from blue sets in &g, we loop over C, in a fixed order (e.g. z-order w.r.t. the underlying p-tree) and pick all elements E., &g until one child provides sufficiently many points to satisfy stopping (i.e. picking all |EE| points would lead to a set of blue points of size at least k).
From this child, the selection implicitly chooses the remaining set F., of k' S kj points with the smallest z-order. Define set Eg- := E., \ E of points that will not be matched e, which can be represented by an integer and a pointer to child u pointer to the (k' + -h element in set Eg. Finally, denote the set of chosen red and blue points by Ya(A) and Y. (B), respectively, and let Y, (A) U Y. (B). Note that (implicitly) choosing these subsets at node e takes O(p?) time,
13

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
though reporting the points in (A)., or in Y, takes O(k- d) time, where d O(log U) is the depth of the p-tree. Obtaining the aggregated sub-instance. For any cell u E Ce define zu to be the center point of u. For æ E Ye, let u(x) be the subcell of C, that x is contained in. We define a mapping m(x) from x to zu(x). to be the function that maps a point x to the center of its (bounding) subcell.
For any point set S, let m(S) == {m(a) | E S). We call mY.) aggregated sub-instance w.r.t to node v and the corresponding point set Y,.
Note that the function m(Y)) can be implicitly represented by the p² sub-grid cells of the cell of and their corresponding points in Y,-
Computing the implicit matching: We solve the following Euclidean transport instance on the points of For every point a' E mY,Y (A)), we set the supply of a' to be the number of points in (.(A) that map to a', i.e.. Se' := 1{ E (A)) mm()) = =. Similarly, for every point D' mY)((B)), we set the demand of b' to be the negative number of points in (B) that map to b', i.e., dy == lb E (B) m(b) b.). Let ne [Y.] be an optimal assignment (or implicit matching) obtained by solving the Euclidean transportation problem the sets m (A)) and mn(Y, using the algorithm of Theorem 2.1
(A)(A)).m(Y(())))
EXPLICIT MATCHING The inputs to the explicit matching procedure are vertex of the restricted p-tree v, the point sets Y, (partitioned by their corresponding cells in E C2 and an optimal assignment 7, YY.] to the Euclidean transport instance on (A))
s, d) described above. Its output is a matching of points Y..
The procedure iterates through the edges (a', E mA))(A)) mY,, (B)) with non-zero assignment in r[[Y.] in an arbitrary manner and adds n. [a] 8') arbitrary pairs of points a, b, where E Y, (A)) and b E (B),), and m(a) a' and m(b)= 6', as edges to the output matching from the cells with centers d' V.
C. Proof of Theorem
Theorem 3.3 (Static). For any Euclidean matching instance (A, B), with |A| = =IB] = and spread U S n", for some constant e > O, and for any constant & E (0, } algorithm STATIC-MATCHING(A. B) computes an O(1/)))approximate Euclidean matching with high probability in O(n 1+5 !e) time.
For this proof, we will assume that both n and are powers of 2. While a similar approach extends to the general case, we present the special case here for the sake of clarity. We show that p = n² yields the runtime and approximation bounds.
Proof of Theorem 3.3
Run-time analysis:
By definition each leaf of p-tree T contains at most p² many points. Now, since, we can compute the exact optimal matching between the points Y, for leaf node in O(Y, v 1s) time using the Hungarian algorithm, the total runtime of computing exact matchings of all leaves is bounded by (p n)..
Next, we discuss the running time of step 1 in IMPLICIT-MATCHINGG1 Y.). By Lemma 3.1 for each non-leaf node v, the MMPLICIT-MATCHNN in the execution of MATCH(U) can be implemented in O(p² + p² log" n) O -oog n) time. In addition, we need to solve AMPLICI-MATCHING for at most nodes in the p-tree T. The total run-time of solving the IMPLICCT-MATCHING at each node is bounded by O(n - p log" n).
For any node v. the total time to report the points in excess set F.g is bounded by O(FF,))- Thus, Lemma 3.2 yields the perfect matching on set X takes (n)) time. Hence, the runtime of STATIC-MATCHING(AAA B) is at most O (n (n+p" n + n. p² log² + n log since the number of levels in T is bounded by log U Setting p = n and logp'
using that U S n for a fixed e > O, we obtain the stated run-time bound.
Note that while the upper bound we present of the approximation ratio holds in expectation, in the static setting we may repeat the algorithm Ollog n. times and output the smallest cost solution to find an aaprrxiaatoo solution with high probability. In the dynamic setting, this simple trick guaranteeing a high-probability approximation ratio does not work.
Correctness analysis: Let n² be matching on set X with optimal cost c(l²".. Say that edge is crossing at level if it crosses the grid of level of T. Observe that an edge may be crossing at multiple levels. Let Cross(M, stand for the number of crossing edges of matching M at level i.
Note that if none of the edges of 1ı" is cut by grid at level e, this implies that all nodes at level have a balanced number of
14

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
blue and red points. Also note that no edge of μz is crossing at level 0, we assume input points are contained in the cell D D for some D n (both in the static and dynamic settings).
Let Pa be the matching computed with our algorithm w.r.t. T. Our goal is to prove that
Bic(H()) = -(=='_e(f") c(("))
where the expectation is taken with respect to the random shift introduced at initialization. Let T2 stand for the vertices of T at level f. Let X = D/pl stand for the length of the sides of the grid level Slightly overloading notation let X stand for the length of side of the cell of
Matching Ma has the following convenient property: every point of Y, is matched with a point of Y, (inside the cell of v) under Pa due to the construction of our implicit-explici matching procedure. Unfortunately, 1a might not have this convenient property, making the comparison between them difficult. We now describe a procedure which translates j²" into a matching z' with this property. Later, we compare the costs of y', ,²" and Ha.
Generating μl from j We describe iterative method which generates yl from μr" Initially, we set jé' = p=" Assume that the tree T has & levels. Iterating& E from l to O for every u E Tr, we keep repeating the following as long as possible:
If there exists edges of p", denoted (r1, b1), (r2, b2), such that rı, b2 E Y, for some v E Te, T1, by are of opposite color and bi, 12 & Y, then replace (r1,b)), (r2,
with (r1, b2), (r2, b2) in μم.
Let pt for Ee E 0l stand for the state of ye after not being able to repeat this process at level E + 1. Specifically, 1x- = Mee
Claim C.1. The above process terminates in n iterations and for any E it must hold that for all v E Tse the points of Y, are matched exclusively to each other by pt-
Proof. Observe that before the execution of any of the two steps none of the two edges run between points of Y,, for any w E T. However, after the execution we have (r1, b2) e pl and Tı, b2 E Yu for some v E T. Hence, the process must terminate in n steps.
Also, note that any Y, contains the same number of red and blue points. If, without loss of generality, any red point of some Y, is matched outside of Y, by p', then there must also exist a blue point of Y, matched outside of Y, and the process can take a step. Hence, once the process can't make any further steps for some Y, where rı, b2 E Y., then pz is matching the points of Y, to each other.
Claim C.2. Cross e') S Cross(ne, for all l' <e
Proof. Assume that the process is progressing through vertices of level E. We will argue that any step it may take will not increase the number of crossing edges at lower levels, implying the claim. If the process takes step for some v E Ti then it adds edges (r1, be), (r2, b1) to pL. (f1. b2) is in the same cell at level é hence it is not crossing at lover levels. (r2, b2) might be crossing at some level E' < f. However, in this case as (r1, b2) are in the same cell of level el it must hold that in this case either (ri, b1) or (r2, b2) was crossing at level before their removal from μ².
We will first relate the costs of matchings μ²" and p!. Claim C.3. clu') S c(u") + Ei 8- d Cross(u", E).
Proof. Consider the state of jı' after having finished with vertices of level E (+ 1, that is μ4z- Call a point bad at this point if it belongs to some Y, : u T3 and is matched to some point outside of Y, by v/:- In this case there are two options: it either matched outside of its cell at level f, or it is matched to point forwarded by node v to its ancestor but still inside C.. Note that it can't be matched to any point of Y,, for a node on higher level then v due to Claim C.1
The set of points forwarded by u to its parent is monochromatic. Hence, for any bad point of the later type we will be able to associate a bad point of the earlier type. This implies, that at least half of the bad points, denote their number by Be, are incident on edges crossing at level in plz that is 2ig has at least Br/2 crossing edges at level f. Observe that,
15

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
due Claim C.2 this implies that J also has least Be/2 crossing edges level t. Once the process gets level & in every step, it will reduce the number of bad edges by Furthermore, we claim that in each step may increase the cost of p by at most 4 Ag. By every step some (r1, b2), (r2: b1) are added to p' and some (1,, (r2, b2) are removed. We have that c(r1, b2) I 2 Re r1 and b2 are contained in the same cell at level L. By triangle inequality, we have c(r2,h) b1) S a(bı, r1) + c(ri, b2) + (2) and therefore c(r2, b1) + c(r,, (b)) - -c(b1, rı) - a(b2, r2) I 2. 2 b2). Hence the total increase is indeed at most 4 dg. That is, the process will take at most b3g steps, in each increasing μ by at most de, where B, S Cross((",,)) f2. Summing over all levels proves the claim.
We now relate the costs matchings y' and pa. Let M stand for the edges of matching M restricted to points Y. By construction of the algorithm and due to Claim we know that both of μ' and Pa match the points Y, to each other for any T. First, observe that if is leaf then our algorithm calculated an optimal minimum cost matching of Y,, that c(p/(Y)) ((Pa(Y))) for all leafs V.
Claim C.4. ((aalYa)) S ((.) + Ao/P- [Y.) for all non-leaf T.
Proof. Recall that the algorithm computes Ha [Y.] by aggregating all points of Y, to the respective center of their sub-cell in C, and computing an optimal solution of the resulting Euclidean geometric transport instance, then moving the points back to their original location. That is, assuming the points were already in the center of their respective cells, H4a [Y,] is a minimum cost matching of Y.. Moving a point to the center of it's respective sub-cell increases the cost of an optimal solution by most \//P (the length of the sides of the sub-cells of v), similarly the cost of moving them back their original position also upper bounded by An/p.
Claim C.5. EveTalo non-lenf Yol $ Cross(p", Cross(p",,+100 1)+ fort l fo.... -1].
Proof. Consider a particular vertex E Te. The points matched at v (the points Y.) are all inherited by by its children nodes due an imbalance of red and blue points in the respective cells of e-s children nodes. This implies that j matched at least |Yy| points in C, with points in separate sub-cell of layer e + 1. Hence, |Y,| points of C. are incident crossing edges at level + 1 in p² Summing over all vertices of (Y,/u Te} proves the claim.
We now have to tie the number of crossing edges to the cost of μ².
Claim C.6. [[rosss(p " ,0)//,))a>eeee = O(c(p")) 00/A>==(<(((")))6)) shift introduced at initialization.
all e kk where the expectation taken with respect to the random
Proof. Consider an edge e E μ of length L. With probability Q(L/Ac) it crosses the grid at level due to the random shift. Hence, it contributes T. to both sides of the inequality.
We are now ready to upper bound the approximation ratio of the algorithm in expectation with respect to the random shift at initialization.
16

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
B[c(pa)]
Σ
E[c(Ho[Di]]] +
[Y.)]
B leaf
vET v non-teaf
s B[c(H1)] +
Σ
2 2- 2-E| >xo/
(2)
vET l0 is non-teuf
k 8 xe E|Cross((u" ,0]+
< dyy") +
E
2 2.E||Yol] Xo/P
(3)
t=
vET is non leaf
R-1
dx²) +
Σ
8 de E(Cosss("' +
Σ
4 xe/p B|Cross(44 , 1)]
(4)
t= .k)
t=
O(d(x)).
(5)
O(dp))/E)
The first inequality follows from the fact that Ma matches vertices of Y, exclusively with each other. Equation holds due
Claim C.4 Equation 3 follows from Claim C,3 while equation follows from Claim C.5
follows from Claim C.6 Finally, the last inequality holds as the number of layers the p-tree is upper bounded by O(log U/ log P) O(1/s).
Equation 5
D. Basic Algorithm
In this section we first present the pseudo-code of our basic algorithm which can answer queries about the solution efficiently (see Theorem 4.T. Afterwards we argue about it's running time and approximation ratio.
D.1. Running Time Analysis
We will first argue about the update time of the algorithm. We will analyze the update time of an insertion the process of handling of deletions is similar. An update is handled in a top-to-bottom manner: the algorithm first calls the update procedure at the root and then always recuses towards the child which contains the newly inserted point. Hence, throughout an update the algorithm will visit nodes of 7 along a path form the root to some leaf. Apart from the visited leaf node the only non-constant time operation executed at any node is the implicit-matching procedure which runs a geometric Euclidean transport sub-routine on an aggregated sub-instance.
These aggregated sub-instances consist of points corresponding to the centers of the sub-cells of the respective node, that is p² points. The weights assigned to each center are at most n. Hence, by Theorem 2. they run in time O(p log(p log n.) ((e)). That is, updating the non-leaf nodes along the path takes ((e) /E) worst case update time, as will always have at most O(1/))-leves.
The leaf node might end up having more then p² input points in its cell after the insertion (otherwise we simply call the Hungarian algorithm on it). This might result in the leaf splitting into a new sub-tree. This sub-tree can similarly only have ((1/E)-level) and in total at most O(p²/e /e) nodes. On each of its non-leaf nodes similarly the only non-constant time operation will be the implicit matching procedure adding (ne)))// /æ) worst-case update time.
There might be~ y² leaves of the new sub-tree, but the number of points inside the cells of the new leaves is at most O(p"). In each such leaf the only non-constant time calculation the algorithm completes is a call to the Hungarian algorithm to obtain an optimal solution inside the leaf, hence this contributes at most O(p" = -(E)) to the update time.
maintain the p-tree T itself in O /E) O(n /E) worst-case update time. Hence, the total worst-case update time of the algorithm is (ne) /s).
Additionally,
D.2. Algorithm Correctness
The analysis of the approximation ratio of the static algorithm presented in Section C relies on the following properties of the output:
17

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
Algorithm (PDATEE,
if e is non-marginal internal node of T then u f child E Cp such that a E X Update Eu and d(uu) via INSERT-UPDATE((UA))
Update [Y.] and E.p via IMPLICIT-MATCHINN((C,C c(py) -Enec, c(pu)) c(][X)]
end if
if is non-marginal leaf in T then (A)) + (A(A) {a}
if r 2 0 then
E, - E, (a]
1) - 7 +1 (A) (A) E.
else
b any point from E.
E, t E, \ {b}
1) + 11 -
(B) + (B) \ E,
end
He
end if
if is new subroot of I' then
if 1 2 0 then
E. - E, fa}
17 - 7 +1
else
b any point from E,
F. F., {b)
17 17
end if
Update [r and c(ux) via MATCH(e) c(pv) + Enec. c((u)) + c(C((V))
end
Algorithm 2 INSERT(a) (a)
Update restricted p-tree T INSERT-UPDATE(r..) a)
18

---

Fully Dynamic Euclidean Bi-Chromasn Matching in Sublincar Update Time
Algorithm 3 DELETE- UPDATEEV, if e is an internal node of then 4 --hild и E T, such that & X Update Eu and c(Ha) via INSERT- USSET-UPDATT(UU,)) Update %o and Eg via (. c(pn) -uuc, c(ux) + (()((V)))
end if
if e is leaf in T then (A) + (((A) \ {a}
a E F, then
E. t E {a}
else
1 > 0 then
b < any point from E.
E, - E. {b}
1) n (B + (B) Ey
else
E E, u {a}
1) 7 +1
(A) & (A) E,
end if
phg < HHNNGARIAN-ALCORRTTMMMAA. (.)
end if
end
Algorithm 4 DELETE(a)
Update restricted p-tree T
MBLETE-UPDATE((R.. a)
(1) For any node E I points of Y, (which all lie in the cell of v) are matched exclusively to each other.
(2) For any leaf node E T the points matched in the cell of e, Y,, are matched via an optimal Euclidean bi-chromatic matching algorithm (Hungarian algorithm).
(3) For any internal node e E T the points matched in the cell of v, Y, are matched via the implicit-explice matching procedures relying on an underlying optimal Euclidean geometric transport sub-routine (Theorem 2.Tp.
(4) For any node v E the algorithm matches all but a monochromatic set of points out of the sets points forwarded to e by its children inside the cell of (or in the case of leaves among all the points the corresponding cell C.).
We will argue that once update is processed by the dynamic algorithm, all of these properties of the output are restored. The only vertices of 7 where any of the mentioned properties might be violated are the ones whose cell contains the updated point or ones which have been created due to the update to the p-tree.
Leaves, which have contained the updated point are created due the update, all select arbitrary maximal color balanced set of points in their respective cell. The algorithm Afterwards, leaves pass on monochromatic set of points their parent internal nodes.
the Hungarian algorithm on these sets restoring item (2).
Every internal node runs the implicit-explicit matching procedure on the points passed to by its children bottom-up order. This implies they will satisfy item (3) once the update has been processed. The implicit-explici matching procedure matches all but monochromatic set of points in cell, which set then forwarded to the parent node. This implies that by the time the root has finished its internal caleulations, both items (4) and (1) are satisfied.
Note that we may guarantee that the approximation ratio holds with high probability if we simply maintain O n) solutions in parallel and always query the one with smallest cost. However, in this case in order maintain the approximation
19

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
ratio with high probability over super-polynomially long update sequences the data-structure has to be periodically rebuilt reducing the update time guarantee to be amortized.
E. Advanced Algorithm
In this section, we describe an improvement of the algorithm described in Section 4 The main difference is that in the advanced algorithm, we maintain an explicit matching at each node. As shown previously, we may recompute the explicit matching at every node efficiently after an update through black-box Euclidean geometric transport algorithm. However, without opening the black-box, we can't 't guarantee that the implicit matching computed will not differ significantly from its previous version. This means that in order to maintain a matching of the actual input points we would have to re-run the explicit matching procedure, but that takes time proportional to the number of points matched at the specific node which could be (n)).
To overcome this difficulty, we show that we can update the implicit and explicit matchings stored at each node v E T in time (p²") when two points are added to the respective point set Y.. We will call this procedure CEEMENT-MATCIING We then show that the size of the excess set of each node of T changes by exactly after each update. We finally present an algorithm which maintains the explicit matchings stored at each node efficiently.
E.1. Updating Implicit/Exllict Matchings Under Insertions
Let be a matching of some red and blue point sets (A, B) and x; y be some arbitrary points of opposite color not in A, B. Define an augmenting path with respect to j to be a path between x, y consisting of edges between x. y, A, B, alternating between edges not in j and in μ²" Define the cost of such a path to be the sum of the costs of its edges not in j minus the cost of its edges in j..
For such an augmenting path II, define the matching p' consisting of all of the edges of μ not in Il and the edges of I1 not in to the matching j augmented by II. Observe that c(1²) = c(μ1 + c(II) and c(f²) is a perfect matching of U {x}, B {v). Claim E.1. Let A, B be a set of input points for the bi-chromatic Euclidean matching problem and μ² be an optimal matching of points A, B. Let x, y be some arbitrary points in the plane of opposite colors and let II the smallest cost augmenting path with respect to p starting and finishing at x and y, respectively. Then the matching obtained by augmenting μ² by 11* is a minimum cost perfect matching of AU {x}, B {y}
Proof. Define (n²) to be the matching obtained by augmenting μ* with IT.. Assume that it's not of minimum cost and let "ig stand for a min cost matching of the extended point set A U {x), BU {v). Look at the union of μ² and μ It consists of disjoint edges, cycles and a single augmenting path II with respect to μ² between x and y. Observe, that the cost of matching all points not part of II' is the same under both matchings due to their optimality in their respective point sets, hence c(pi)) = c(p²) + e(II') However, c((x²))) c(p²) + c(II') this would imply that c(II*) > e(TT') which is contradiction.
Now extend the problem bi-chromatic Euclidean matching to point sets where certain points might have the same coordinates. It can be thought as extending the problem to inputs where input points have multiplicities.
Claim E.2. Let A, B a set of red and blue points where some pairs of points might share locations and μº be a min cost bi-chromatic Euclidean matching between them. Assume that points A, B are located on k disjoint points of the plane, and points of A U B at the same location are monochromatic. For arbitrary red and blue points I, y, there exists a min cost augmenting path between x and y with respect to μ² which visits each of the k disjoint locations at most once.
Proof. For sake of contradiction assume that there is such min cost augment path and let x' be minimum length (in terms of number of edges) min cost augmenting path between r and with respect to μ². We will show that we can 'shortcut' T".
Assume that II visits one of the k disjoint locations twice. Then must be the case that the II contains sub-path C' which is cycle starting and ending with consisting of n² and non-ju" edges. There are three possible cases.
If the total cost of non-jr" edges of C' is the same as that of the μ edges, then we may remove this cycle meaning that IT" is not of minimum length.
20

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
If the total cost of non-jr edges of C' is larger then that of the μ edges, then if this cycle is removed from the path, we would get smaller cost augmenting path between and y-
If the total cost of non μ¹ edges of C' is smaller then that of the μº" edges, then the matching (n" |C) U (C matching received by augmenting μ² by C)) is of smaller cost then j and is perfect matching of A, B.
(1*) (the
Theorem E.3. (Bellman 1958) Belman-Ford: There exists an algorithm which given a graph G = (V, E) and distance function c: E + R where G doesn't 't contain a negative length cycle with respect to distance measure c, returns shortest path between two vertices of V in time O(V) |E))
Augment-Matching Procedure We are now ready to describe the procedure which augments the implicit matchings under vertex updates, while making sure only small number edges are affected by the update.
Recall that an implicit matching >" ["_), of v E T corresponds to a matching of the aggregated versions of points Y,.. Assume that points (x, y) are inserted into Yo. Define the following simple graph G (V, and corresponding edge distance function c: define V to consist of p² nodes, Ve representing the center pints of the sub-cells of v, and two additional nodes x',y'.
Without loss of generality assume and are red and blue respectively. Add an edge between x' and all nodes of ve E Ve to E where corresponds to sub-cell of containing only blue points in Y, (recall that the points of Y, in every sub-cell of v are monochromatic). Set c((x', vc)) to correspond to the distance between the centers of sub-cells of x' and vc.
Similarly, add an edge between y' and all nodes of vc E Ve to where v corresponds to sub-cell of v containing only red points in Y,, and set c((x', vc)) to correspond the distance between the centers of sub-cells of x' and vc-
In addition add an edge to E between any two sub-cell center nodes ve, we if their sub-cells contain points of the opposite color. Set c((re, Wc)) to be the distance between the respective sub-cell centers if there isn't 't an edge between the centers in > [Y=]. and otherwise set it to be the negative of this distance.
By construction, path between z', y in G' corresponds to augmenting path of 7 [Yo] between a, y. By Claim E.2 the shortest x', y in G' corresponds to minimum cost augmenting path between with respect to [Y]].. Observe that a negative length cycle of G could be used to reduce the cost of "[Y)) which was assumed to be optimal, hence no such cycle should exist.
Finally, due to Claim E.1 we may augment > [Xo] to be min cost perfect matching of the aggregated versions of Y, U UZr w) with such shortest ,/ path.
Constructing graph G may take O(IV] + [F])= = O(p2) time. We may find the shortest path in O(p²) time by Theorem E.3 Any shortest path in G can be of length O(p²) as G doesn't contain negative length cycles. Hence, we can update the implicit matching of Y, in O(p²) time. Once the the updated implicit matching is computed we can transform it to an explicit matching of Y, U {ar, y) in time O(p²) since the difference between the old and updated state of y [Yr. (and hence the corresponding explicit matching) is only O(p²) edges. We state this formally by the following lemma:
Lemma E.4. Let T be a p-tree and node in T, 7" [Y,]] the implicit matching corresponding node u, and x. y & Yo [Yn], two points of opposite colors. AUGMENT-MATCHING [Y.), x.,y) v) returns an optimal implicit mutching of the aggregated
versions of Y, U (z.v). in O(p²) time.
Note that while Lemma F.4 can be easily be extended to point pair deletions, there is no need for this. Assume we want to update an implicit matching "[Y]]] when points (x. y) are deleted from Y,. Let the neighbors of x, y under > [Y,] be x', y.. We may remove edges (x, x²) and (y, y²) from 7 [Y,) and re-insert point pair (x',y/) ,/)) using AUGMENT-MATCHING to simulate this deletion.
E.2. Dynamic Algorithm
We now describe INITIALIZE, the update procedures INSERT and DELETE, as well as QUERY-RECOURSE The procedure QUERY-COST is the same as in the basic algorithm. We will only describe how to do point insertions, as the procedure for point deletions is almost identical. Similarly as before, given p-tree T, our goal is to maintain sets E,, (B)). difference T) matching 7 [Y.), but also matchings p! (Y.) and Mv, together with its corresponding cost c(14,) for all the
21

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
nodes in T which is the key difference from the basic algorithm presented in the main body of the paper. We first begin by describing the initialization procedure.
The procedure INITIALIZE is almost the same as in the basic version of the algorithm. However, in the advanced algorithm, each node of T, we also convert the implicit matching to its corresponding explieit matching. This further allows us maintain the complete matching of subtree corresponding rooted each of its nodes, defined as Ma, here C, denotes the children of node in p-tree T,
Initialization
Me [Y.] U.wec,
Handling Insertions Without loss of generality, we assume that we are inserting red point to A. As in the basic algorithm, INSERT proceeds in two steps: updating the restricted p-tree and updating the matchings at each affected node. We define the recursive procedure INSERT-UPDATE(1, a) for three types of nodes in T: new subroots, non-marginal leaf nodes, and non marginal internal nodes.
If a leaf e is the new subroot of T, the INSERT-UPDATE(;, a) procedure is the same as in the dynamic algorithm presented in Algorithm[] I
For non morgmnala leaf nodes v, the procedure INSERT-UPDATER a) checks if there exists a pair of points of opposite color at v, and if so, augments the matching at v.
For non-marginal internal nodes v, we have two cases; the excess set E., of its affected child either increases by one red point a'. or decreases by one blue point b'. In the first case, we proceed similarly as in the leaves. Namely, if there exists a pair of points of opposite color at v; and if so, augments the matching at e. In the second case, we first have to check if the blue point b' was a part of a matching at v. so, we first have to break that matching. After this, we proceed similarly to the first case.
Finally, we prove the following theorem:
Theorem 1.1. For any 0 < & S 1, there exists a fully dynamic algorithm that maintains an expected ) aapproximate solution to the Euclidean bi-chromatic matching problem defined on the point-sets A, B c R2, |A| = BB), while pairs of points are inserted into and deleted from A, B in O( 111 worst-case update time. Here, n denotes the maximun size of A throughout the update sequence.
As the proof of this theorem is very similar to the proof of Theorem 4.1 we only note the key differences between these two algorithms.
Proof of Theorem T.1
Running time: Compared to the basic algorithm, the running time of INITIALIZE only increases by the time required to convert implicit matchings at nodes of T to their corresponding matchings of the input points, and therefore also store the complete matchings of the subtrees. By Lemma 3.2 this takes additional O(n) time, making the total run time of the procedure still +sn
Now we discuss updates. As mentioned before. the difference between the basic and the advanced update algorithm is that now, instead of recomputing the entire implicit matching at each affected node (using CEGEMENT-MATCHNNGG we use AUGMENT- MTTCHING procedure to augment both the implicit and the aggregated matching, if needed. The computation of AUGMENT- BUGMENT-MATCHING MATCHING takes O(p²) time by Lemma E.4 so we again obtain that (nn)) is running time required for the update procedure of this algorithm.
Regarding queries, in the advanced algorithm, we can easily report the perfect matching by simply returning the complete matching of subtree at root r of T, that is py. Further, we obtain more efficient algorithm for reporting the recourse of the matching. Namely, note that the aggregated matching at each affected node changes only along the augmenting path found in AUGMENT-MATCHING which is On /E) for p = fo.
Since the length of this path is at most p², the total recourse is at most O(p² log U log )
Correctness: In terms of correctness of the approximation, the advanced and basic versions of the algorithms behave exactly the same. That is, the advanced algorithm also satisfies the following. For any Y. for E T points of Y, are matched exclusively to each other. Further, for leaf v E T the matching of Y, is computed via an optimal algorithm, for any loa-leaf E I the matching of Y2 is computed by the implicit-explicit matching procedures. Finally, for any cell of the grid, the number of points not matched within that cell exactly corresponds to the the color dis-balance of the points present in that cell (see the similar argument concerning the basic algorithm see Section D.2). As argued in the analysis presented in Section C the algorithm correctness only relies on these properties of the output matching.
22

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
Algorithm 5 DASERT-UPDATE-((U.d
1: if is non-marginal internal node of 7 then
2: + child C2 such that E Xx.
3: E + Eu
4: Update E, via INSERT-UPDATE((U, u)
5:
if Eu = E U {a'} for some E A then
6:
TJ 2 then
E. t E. Ud'
8:
else
9:
+ any point from Eg
10:
Update y" [Y.] and j([o. via ATGMENT-MATCHING(Y" "
11:
end
12:
end
13:
if E. = E. \. } for some E B then
if n. 2 then
a' point from A such that (d'b6) E x(()) Delete (a',b') from [Y,] >"(Y,] b') < x-x [X.](d''' E. + U{a}}
16:
17:
18:
19:
else
20:
if E, then
Ey F. (b)
22:
else
23:
d' from A such that (a',8) [Y)] Delete (a', b') from ['Yo]]
24:
25:
Y (" 8) < > [Y,](a',
26:
8" + any point from E,
27:
Update > [Y,] and p'[Y]] via AUCMENN-MAXCHUNN((1X.....
F, Ew {b")
end if
30:
end if
31:
end if
32:
if E.,
E \. (b} for some E B then 1/ was matched in 7 [Y)] with some a' E A then
33:
34:
1/" < any point from F.g
35:
Update >" [Y,] and y u'[Y.] via AUGMENT--ATCHINNNTT [K,], a'), d',
end
37: end if
38: end if
39: if v is non-marginal leaf in T then
40:
(!(A) + (-) [
41: Th + (Y) (B)
42:
if 2 then
F.g < F, U {a}
else
45:
b any point from E, Ph ALCMENT-MATCHING((N.. E. < E, U (b}}
46:
47:
48:
end if
49:
end if
23

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
50:
if is new subroot of T then
51:
if 17
2 0 then
52:
E. < E, U {a}
53:
4 1-1
54:
else
55:
fー any point from E,
E, - E {b} (| 1-
57:
58:
end if
59:
Update [Y,) and c(pp) via MATCH(")
60:
c(pb) + Euec, c(pu) + c(1)lx)))
61:
end if
Algorithm (a)
1: Update restricted p-tree? T
2: INSERTUUPDATE(Y
F. Hardness Bound for Partially Dynamic Bi-Chromatic Euclidean Matching
In this section we present a simple construction which shows that even under point insertions no efficient dynamic algorithm may maintain arbitrarily tight approximation the bi-chromatic Euclidean matching problem. This is in contrast with the static setting, where multiple (1 approximmtee near-optimal running time algorithms are known (Agarwal et al.
2022b Raghvendra and Agarwal 2020).
Theorem F.1. For any 0, any (2 - 8) approxxmtttee algorithm that maintains a solution to the Euclidean bi-chromatic matching problem under either insertions or deletions of pairs of points must have (n)) amortized update time.
Proof. Consider an arbitrary dynamic algorithm which maintains a solution to the bi-chromatic Fuclidean matching problem in o(n) amortized update time. We will describe an instance with initially ((n) points and an insertion sequence consisting of further (n)) points such that when fed to the algorithm must output an at least (2 -a-pprroxeeate solution at some point during the insertion sequence for any arbitrarily small constant d 0. To extend the argument to algorithms fir to only handle point deletions may just reverse the insertion sequence.
The initial input consists of n red and blue points on a line. Specifically, it consists of red points A = (-n. n 2, n + 4..... n - 2). and blue points = (-n + 1, - 3,..... n 1). that simply a series of points with alternating colors with unit distance between them. Then, the next n operations consist of inserting pairs of points at positions -n and 1+ i, for i (I,..., n1,, such that at each time the ordering of points on the line is alternating between red and blue points. In other words, with each point pair insertion maintaining that the sequence of points are alternating in color. Observe, that at any point in time we may match the points greedily from left to right obtaining a matching of cost n + i after the i-th insertion.
simply extend the sequence by point each side while
Consider the i-th point pair insertion throughout the point pair insertions where insert a pair of points coordinates -x and x for x n n+ for some E [n). Let 14, and nt stand for the edges of the matching output by the algorithm before and after this specific insertion was handled respectively.
In the (multiset) 14; U n every point present before the i-th insertion incident on two edges where as the newly inserted points are only incident on one edge of μ.- Therefore, there exists a path II, μi nr," connecting the newly inserted pair of points. As the algorithms amortized update time is o(n) over the series of updates the algorithm has o(n2) time to maintain its output, hence there must exist some i [n] such that [II, o(n).
As the new pair of points are inserted at coordinates = and x - 1 it must hold that c(Il,) 2 2- x - 1 2- (n + i) 1. Furthermore, both n and μ consist of least n + - edges of length at least 1. Hence,
24

---

Fully Dynamic Euclidean Bi-Chromatic Matching in Sublincar Update Time
((1)++((
(Il.) + c(p. \II.) + c(p \ \1.)) 2 (n + - -1+ + lui + 1 - [11
4 (n + 1) o(n)
Therefore, one of μ or Si must have cost at least 2 (n+ i) -o(x). As the cost of an optimal solution after the i-th insertion i of these matchings must be least (2 2-00 -appoxaato proving the theorem.
25


Reviews:
[{"review_id": "IOPxWf3gdP", "reviewer": "Reviewer_MBvJ", "summary": "The paper studies the euclidean minimum cost bipartite matching problem: given n blue points and n red points in the 2d euclidean plane, we wish to compute a minimum cost bipartite matching between them, where the cost is measured in terms of euclidean distance. The novel component of the paper is to introduce dynamic updates to the point sets while requiring that an approximately optimal solution is always maintained. The trivial baseline is to recompute the matching after every update which can require linear time every update for a $1+\\epsilon$ approximation. The paper presents a tradeoff which leads to much faster update times. The paper obtains an update time of $n^{\\epsilon}$ but maintains a worse $O(1/\\epsilon)$ approximation. \n\nThe main idea seems to be inspired by the $1+\\epsilon$ approximation algorithm in the static case. We first divide the input space into many nested regions using a quad-tree like data structure, but the fan out of the tree must be carefully controlled. When points get updated, we try to match from bottom up, but it is not clear if at some point $\\Omega(n)$ updates need to be made. However, a careful update step introduced by the authors only needs to do work proportional to the fan out of the quad-tree like data structure.\n\nEmpirical evidence is given showing that their algorithm can obtain multiple orders of speed improvement over the naive recalculating baseline while maintaining an accurate solution, in both synthetic and real world data.\n\nOverall, while the setting is a bit limited, I think it is a very solid contribution to an important problem.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their time and for providing valuable comments. \n\n- “Do the authors have an estimate for the constant in the $O(1/ \\varepsilon)$ approximation. I understand that asymptotically it doesn’t matter but it would be nice to understand what the factor is for say $\\sqrt{n}$ update time.”\n\n\nWe haven’t done an exact analysis of the constant in the $O(1/ \\varepsilon)$ approximation guarantee. However, if one would want to achieve $O(\\sqrt{n})$ update time, we expect the approximation ratio to be $< 20$.\n\n\n- “Does the analysis work for weighted matchings e.g. something like optimal transport? Of course one can replicate the points by their weight but it leads to some blowup in the parameters.” \n\nThis is a great question. We believe it works and our algorithm should extend to this setting, but details need to be worked out.\n\n- “The case of $\\varepsilon \\approx 1/ \\log⁡(n)$ seems especially interesting since one gets $O(\\log ⁡n)$ update time. Is there any hopes of improving the running time in this regime?”\n\nThis is a natural question to consider, but we’re not aware of any techniques of pushing down the runtime to something sub-logarithmic.\n\n- “Can the authors quickly survey what is known about the dynamic problem in high dimensional settings? i.e. when $d$ is not constant? What is the dependency on $d$ of the current algorithms? Does an exponential factor in $d$ show up in the approximation or in the update time?” \n\nTo the best of our knowledge, the dynamic version of our problem in higher dimensions hasn’t been studied yet. However, we can easily extend our algorithm to work in any dimension $d$. The best possible trade-offs need to be worked out, but we believe that our approximation ratio has a multiplicative factor of $O(\\sqrt{d})$ and the update time should be of the form $O(n^{\\varepsilon d})$.", "claims_and_evidence": "Yes the proofs seem convincing and correct.", "methods_and_evaluation_criteria": "I am not an expert but the experiments seem sound and show many magnitude speed gains over the naive solution which recomputes the matching every time for both real world and synthetic datasets.", "theoretical_claims": "I checked to the best of my ability.", "experimental_designs_or_analyses": "I did not check very carefully.", "relation_to_broader_scientific_literature": "No", "essential_references_not_discussed": "None that I know of.", "other_strengths_and_weaknesses": "I appreciate the authors attempt to make cler a complicated algoirthm in figure 1", "other_comments_or_suggestions": "None.", "questions_for_authors": "- Do the authors have an estimate for the constant in the $O(1/\\epsilon)$ approximation. I understand that asymptotically it doesnt matter but it would be nice to understand what the factor is for say $\\sqrt{n}$ update time. \n- Does the analysis work for weighted matchings e.g. something like optimal transport? Of course one can replicate the points by their weight but it leads to some blowup in the parameters.\n- The case of $\\epsilon \\approx 1/\\log(n)$ seems especially interesting since one gets $O(\\log n)$ update time. Is there any hopes of improving the running time in this regime? \n- Can the authors quickly survey what is known about the dynamic problem in high dimensional settings? i.e. when d is not constant? What is the dependency on $d$ of the current algorithms? Does an exponential factor in $d$ show up in the approximation or in the update time?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ru7CeXCOMu", "reviewer": "Reviewer_zchS", "summary": "This paper gives an algorithm for dynamic bi-chromatic matching in Euclidean space with $O(\\frac{1}{\\varepsilon})$-approximation ratio and sublinear update time $O(\\frac{n^\\varepsilon}{\\varepsilon})$ with theoretical guarantee and this algorithm is the frist sublinear update time algorithm for geometric dynamic bi-chromatic matching. In addition, this paper gives a proof of approximation lower bound for dynamic bi-chromatic matching in Euclidean, which shows that even if only insertion and deletion operations are considered, there does not exist a dynamic algorithm that implements both the $2-\\delta$ approximation and the sublinear update time.\n\nThe main technique is based on the static algorithm by [Agarwal-Varadarajan 2004], the idea is to construct a nested grid cells structure(p-tree) to get an $O(\\frac{1}{\\varepsilon})$-approximation matching and maintain the implicit matching on the tree with bottom-up manner. The theoretical results are validated with experiments on real-world datasets.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their time and for providing valuable comments.\n\n- \"The contribution is not clear enough in relation to [Xu-Ding 2024].\"\n\nThe Xu-Ding paper studies a slightly different problem: dynamic maintenance of optimal transport. Specifically, the goal is to design a data structure that accelerates executing 'one iteration' of the Network Simplex Algorithm on an explicit complete bi-partite graph in $O(n)$ time, where $n$ is the current number of points. Their work is concerned with maintaining an exact solution, within numeric computing accuracy, after the insertion/deletion of a vertex to the min-cost flow problem. Since their algorithm operates on a complete, bipartite graph, it requires $\\Omega(n^2)$ space to store the edge weights and $\\Omega(n)$ update time to insert a new vertex. Moreover, in the worst case, the update time can take up to $O(n^2)$ time. As such, the algorithm and the implementation cannot handle very large instances in practice. \n\nIn contrast, our paper studies the problem of quickly maintaining an approximate solution to the Euclidean bi-chromatic bipartite matching (i.e., optimal transport with uniform demands), in time that is *sub-linear* in $n$. For example, taking $\\varepsilon = 1/2$ gives an update time of $O(\\sqrt{n})$ for updating an $O(1)$-approximate matching. Our proposed tree structure uses $O(n)$ space, regardless of $\\varepsilon$, which makes the approach practical for very large inputs. Moreover, our lower-bound in Thm. F.1 shows that the approximation ratio must be at least $2$ for any dynamic algorithm to achieve a sublinear update time.  \n\nIn the related works section, our focus has been to compare our contributions with the known approximation algorithms that run in near-linear time and space, as those are also applicable to larger instances. That being said, as per the reviewer’s suggestion, we will include the above comparison to the Xu-Ding work in the next version of the paper.\n\n**Experimental comparison with [Xu-Ding 2024]**\n\nWe ran the code of [Xu-Ding 2024] on our benchmark dataset (Unif-Vs-Gaussian). We have observed that the insertion time slowed down substantially after inserting $1$k point-pairs (source and sink). In fact, already inserting $5$k point-pairs took more than $12$GB RAM and longer than $3$ hours in total, which is when we had to terminate the process as it was running out of memory. Note that this is already $>1$s per update on a small instance. \n\nIn contrast, our experimental results show that the update time of our proposed algorithm remains $<1$ms for very large instances (up to $n=10^6$ in Figure $2$; right), with the quality of the solution being within a factor of two from the optimum (see Figure $5$; right).", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "Yes", "experimental_designs_or_analyses": "Yes.\n\nIt is mentioned that [Xu-Ding 2024] studied the same problem. Thus, it would be desirable for this paper to have an experimental comparison against Xu-Ding 2004), specifically, the update time comparison.", "relation_to_broader_scientific_literature": "This paper states some connection between this question and the 1-Wassertein distance estimation (e.g. we can apply this matching to estimate 1-Wassertein distance measure). Besides, the paper presents experimental results on real-world datasets and shows how the algorithm performs when we use it to estimate the 1-Wassertein distance.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "Strengths:\n* The algorithm is simple and intuitive, with a solid theoretical guarantee.\n* The problem of bi-chromatic matching is well-studied and practical, and this paper gives some nice approaches to analysis that may be extended to some related work.\n\n\nWeaknesses: \n* The contribution is not clear enough; e.g., the approximation ratio of [Xu-Ding 2024] is not mentioned and used to compare with the theoretical results in this paper.", "other_comments_or_suggestions": "No", "questions_for_authors": "No", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {"ethics_expertise_needed": ["Other expertise"]}}, {"review_id": "8sLUi6QnnP", "reviewer": "Reviewer_PFgM", "summary": "This papers presents a dynamic data structure that maintains a bipartite matching and supports insertions ans deletions of points. Given two point sets in $2$-dimensional space and a parameter $\\varepsilon>0$, the data structure computes an $O(1/\\varepsilon)$-approximate matching and handles updates in $O(n^{-\\varepsilon})$ time. \n\nThe data structure stores a $p$-tree, which is a hierarchical partitioning similar to quadtrees but with a higher fan-out of $p^2$. For each cell of the hierarchical partitioning, the data structure maintains an \"implicit matching\", which can be then converted to an explicit matching in linear time. The data structure gathers the excess demand/supply of the points in the sub-tree of that node at the center of cell, which are then matched to the excess masses at the centers of the sibling cells. Inserting/deleting a point into/from each set requires updating the matching computes for the cells along the path from the root to the leaf node containing the inserted point. By picking $p=n^{-O(\\varepsilon)}$, the height of the tree will be $O(1/\\varepsilon)$ and updating the matching of each cell would take $O(n^{-\\varepsilon})$ time.\n\n## update after rebuttal\nI keep my positive score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their time and for providing valuable comments. For the minor comments and typos not addressed below, we will incorporate them in the next version of the paper.\n\n- “Missing the impact statement.”\n\nWe will include the following impact statement, and we apologize for our oversight: \n\nThis paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.\n\n- “In your experimental results, you used $p=2,8,32$, but your paper mentions the value $p$ as $n^{−\\varepsilon}$, which is less than $1$. I assume that the diameter of the space plays a role in the choice of $p$ in your experiments, but I could not get a sense of what is in fact the value of $\\varepsilon$.”\n\nTo clarify, we choose the value of $p$ as $n^\\varepsilon$ in the paper and obtain an approximation ratio of $O(1/\\varepsilon)$. Here are some examples to give a sense of the scale of $\\varepsilon$: \n\n- our largest experiment had $n=1 000 000$ and $p=2$ corresponding to $\\varepsilon \\approx 1/20$, \n\n- our uniform synthetic data set where points were integers between $0$ and $500$, the choice of $p=32$ corresponds to $\\varepsilon \\approx 0.55$. \n\nEmpirically, all choices of $p$ achieve a significantly better approximation ratio in practice than our worst-case theoretical guarantee (see Figure 5 in Appendix G).", "claims_and_evidence": "The paper introduces a fully dynamic algorithm for Euclidean bi-chromatic matching that achieves an $O(1/\\varepsilon)$ approximation with sublinear update time. The authors prove that no algorithm can achieve a $(2-\\delta)$-approximation while maintaining sublinear updates, establishing a fundamental trade-off. Their method leverages a hierarchical grid-based partitioning scheme to efficiently handle insertions and deletions in $O(n^{-\\varepsilon})$ time. Through experiments, they show that their algorithm significantly outperforms static recomputation methods while maintaining high accuracy. Real-world applications, such as tracking spatial distribution changes in taxi pickup/dropoff data, demonstrate its effectiveness in monitoring Wasserstein distance and evolving datasets.", "methods_and_evaluation_criteria": "The authors propose a hierarchical grid-based partitioning approach that processes updates in a bottom-up manner, ensuring efficient maintenance of an approximate bi-chromatic matching with provable guarantees. Their evaluation, which includes theoretical analysis, synthetic benchmarks, and real-world datasets (such as taxi pickup/dropoff locations), is well-designed to demonstrate both the algorithm’s efficiency and its practical applicability in tracking dynamic spatial distributions.", "theoretical_claims": "I checked the correctness of the algorithms and lemmas.", "experimental_designs_or_analyses": "The experimental results makes sense.", "relation_to_broader_scientific_literature": "The problem of computing bipartite matchings in a dynamic setting has not been extensively explored in the literature, as there are few known algorithms for it. However, I believe this problem could have potential applications in machine learning methods.", "essential_references_not_discussed": "All relevant references that I am aware of are included.", "other_strengths_and_weaknesses": "The paper is well-written and easy to follow.", "other_comments_or_suggestions": "* Line 95 LC: The sentence needs reformatting\n* Line 114 LC: There are two O() notations\n* Line 119 LC: input puts points -> input points\n* Line 168 RC: spread bounded $U$ -> bounded spread $U$\n* Line 411 RC: use a use a -> use a\n* Missing the impact statement", "questions_for_authors": "* In your experimental results, you used $p=2, 8, 32$, but your paper mentions the value $p$ as $n^{-\\varepsilon}$, which is less than $1$. I assume that the diameter of the space plays a role in the choice of $p$ in your experiments, but I could not get a sense of what is in fact the value of $\\varepsilon$.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3H9YLH8H5w", "reviewer": "Reviewer_45u7", "summary": "This paper studied geometric bi-chromatic matching (aka the optimal transportation problem) for discrete distributions in the dynamic setting, where points could undergo insertions and deletions. Here, we are given $n$ points inserted and deleted dynamically, and we want to always maintain an approximation to the optimal matching between the red and blue points. \n\nThe main result of the paper is a dynamic algorithm that maintains an implicit representation of a $(1/\\varepsilon)$-approximation to the optimal bi-chromatic matching in $O(n^{1+\\varepsilon})$ pre-processing time and $O(n^{\\varepsilon}/\\varepsilon)$ update time. The paper further complements the algorithm with a lower bound of $\\Omega(n)$ update time for $<2$ approximation, indicating that we could not expect $1+\\varepsilon$-approximation as in the offline setting.\n\nThe paper further conducted experiments, showing that the proposed algorithm has significant advantage for the efficiency over the static algorithm.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their time and for providing valuable comments.\n\n- “The texts in Figures 2 and 3 are generally very small. My suggestion would be to have the same figures with better resolutions in the appendix and add a pointer.”\n\nAll figures will be replicated in the appendix, allowing for a more readable presentation.\n\n- “Do you have any opinion on how tight your result is? Is it possible to, say, get any approximation in $\\text{polylog } n$ update time?”\n\nFor a choice $\\varepsilon = 1/ \\log n$, our algorithm reports an $O(\\log n)$-approximation solution in $O(\\text{polylog } n)$ update time. We believe that there are instances where our current algorithms attain the claimed trade-off between approximation ratio and update time. It’s an interesting open question whether such a trade-off can be improved, and will likely require new ideas.", "claims_and_evidence": "Yes. The algorithm is quite natural, and the detailed descriptions and the proofs are included in the appendix.", "methods_and_evaluation_criteria": "Yes. Since this paper is the first to consider dynamic bi-chromatic matching, comparing it with the offline algorithm seems reasonable.", "theoretical_claims": "The main technical idea of the paper is to take advantage of the $p$-tree structure to maintain an implicit representation of the matching. Here, if we pick $p=n^{O(\\varepsilon)}$, the $p$-tree will of of depth only $O(1/\\varepsilon)$. The algorithm first tries to resolve matching inside each leaf node; if there is any surplus of a color, the algorithm aggregates the ‘surplus’ and matches the mass to sibling nodes in the $p$-tree. The process could be done recursively only on $O(1/\\varepsilon)$ nodes per point update (since we only need to follow a path from leaf to root), which results in a sublinear update time algorithm.\n\nThe design of the algorithm is quite natural and straightforward. Due to the simplicity for understanding, I did not check the details in the proof. I’m confident the algorithm is correct.", "experimental_designs_or_analyses": "Yes, nothing looks bad to me.", "relation_to_broader_scientific_literature": "Bi-chromatic matching (also known as optimal transportation) has a wide range of applications in machine learning. Solving the problem in the dynamic setting is an important problem in my opinion.", "essential_references_not_discussed": "No essential reference is missing, as far as I know.", "other_strengths_and_weaknesses": "My overall opinion of this paper is positive. Optimal transportation is an important problem in machine learning, and the dynamic setting is an important scenario increasingly common in modern applications. I’m surprised the problem was not studied before. From a technical perspective, the problem also correctly identified the gap between geometric matching vs. graph matching and obtained sublinear update time dynamic algorithm for the former. The paper is also well-written, and I could follow most of the arguments without checking the details.\n\nThe techniques are relatively straightforward, and this might become an issue for technical contributions. However, I personally like simple and cute algorithms and do not have a problem with them.", "other_comments_or_suggestions": "The texts in Figures 2 and 3 are generally very small. I had to zoom in 400% to read. I understand this might be a problem with the template. My suggestion would be to have the same figures with better resolutions in the appendix and add a pointer.", "questions_for_authors": "Do you have any opinion on how tight your result is? Is it possible to, say, get any approximation in $\\text{polylog}(n)$ update time?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper presents the first dynamic algorithm for Euclidean bipartite (bi-chromatic) matching with provable approximation guarantees and sublinear update time. The algorithm maintains an approximate minimum-cost matching between two dynamic point sets in low-dimensional Euclidean space (e.g., ℝ²), under insertions and deletions. The core contribution is a hierarchical space partitioning approach that enables matching to be maintained efficiently via a bottom-up process.

All reviewers are positive about this paper.  They appreciated the result, empirical validation and many applications of the problem in machine learning and beyond.
Decision: Accept (oral)