A Combinatorial Algorithm for the Semi-Discrete Optimal Transport Problem
Authors: Pankaj K Agarwal
Sharath Raghvendra
Pouyan Shirzadian
Keegan Yao
Venue: NeurIPS 2024 poster (2024)
Keywords: Semi-Discrete Optimal Transport
Approximation Algorithm
Cost-Scaling Framework
Primary area: optimization

Abstract: Optimal Transport (OT, also known as the Wasserstein distance) is a popular metric for comparing probability distributions and has been successfully used in many machine-learning applications.
In the semi-discrete $2$-Wasserstein problem, we wish to compute the cheapest way to transport all the mass from a continuous distribution $\mu$ to a discrete distribution $\nu$ in $\mathbb{R}^d$ for $d\ge 1$, where the cost of transporting unit mass between points $a$ and $b$ is $d(a,b)=||a-b||^2$. When both distributions are discrete, a simple combinatorial framework has been used to find the exact solution (see e.g. [Orlin, STOC 1988]). 
In this paper, we propose a combinatorial framework for the semi-discrete OT, which can be viewed as an extension of the combinatorial framework for the discrete OT but requires several new ideas. We present a new algorithm that given $\mu$ and $\nu$ in $\mathbb{R}^2$ and a parameter $\varepsilon>0$, computes an $\varepsilon$-additive approximate semi-discrete transport plan in $O(n^{4}\log n\log \frac{1}{\varepsilon})$ time (in the worst case), where $n$ is the support-size of the discrete distribution $\nu$ and we assume that the mass of $\mu$ inside a triangle can be computed in $O(1)$ time. Our algorithm is significantly faster than the known algorithms, and unlike many numerical algorithms, it does not make any assumptions on the smoothness of $\mu$.
As an application of our algorithm, we describe a data structure to store a large discrete distribution $\mu$ (with support size $N$) using $O(N)$ space so that, given a query discrete distribution $\nu$ (with support size $k$), an $\varepsilon$-additive approximate transport plan can be computed in $O(k^{3}\sqrt{N}\log \frac{1}{\varepsilon})$ time in $2$ dimensions. Our algorithm and data structure extend to higher dimensions as well as to $p$-Wasserstein problem for any $p \ge 1$.

A Combinatorial Algorithm for the Semi-Diserete Optimal Transport Problem
Pankaj K. Agarwal', Sharath Raghvendra', Pouyan Shirzadian', and Keegan Yao
Duke University, North Carolina State University, *Virginia Tech
Abstract
Optimal Transport (OT, also known as the Wasserstein distance) is popular metric for comparing probability distributions and has been successfully used many machine-learning applications. In the semi-discrete 2-Wasserstein problem, we wish to compute the cheapest way to transport all the mass from continuous distribution μ to discrete distribution in Rd for d 1, where the cost of transporting unit mass between points a and is d(a, b) = |a - bll?. When both distributions discrete, simple combinatorial framework has been used find the exact solution (see e.g. [Orlin, STOC 1988)). In this paper, we propose combinatorial framework for the semi-discrete OT, which can be viewed an extension the combinatorial framework for the discrete OT but requires several new ideas. We present a new algorithm that given μ and in R2 and a parameter 8 >> 0, computes an e-additive approximate semi-discrete transport plan in Onn log n log !) time (in the worst case), where the support-size of the discrete distribution and we assume that the mass of inside a triangle can be computed in O(1) time. Our algorithm is significantly faster than the known algorithms, and unlike many numerical algorithms, it does not make any assumptions on the smoothness of μ. As an application of our algorithm, we describe a data structure to store large discrete distribution (with support size N) using (N) space so that, given a query discrete distribution v (with support size k), an E-additive approximate transport plan can be computed in Olk log 1) time in 2 dimensions. Our algorithm and data structure extend to higher dimensions as well as to p- Wasserstein problem for any p 2
Introduction
Optimal Transport (OT) is powerful metric for comparing probability distributions and is used in many machine-learning applications. The semi-discrete optimal transport problem asks for the cheapest transport plan transport mass from (a possibly continuous) distribution j that stored compactly (say, using a deep neural network) to a discrete distribution v. In recent years, the semi- discrete OT has been used data mining (31], image processing 28
29
33
40 computational
biology [53] variational inference (7), blue noise generation 21
51 optics
46
solving PDEs [30]
34], and generative models 9
More formally. a semi-discrete transport plan T between a continuous probability distribution compact support A Rd and discrete distribution with support set B of n points
defined
IRd is distribution over x B whose marginals are dominated by y and v, i..., T: Ax B - Rg0 is transport plan between y and y where Eoe n T(a, b) & p(a) for all A and JA T(a, b) da < vbb) for all b E B. A transport plan is complete if r transports all mass of f to v. For any fixed 2 1,
"The authors are listed in alphabetical order.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

the semi-discrett p- p-Wasserstein problem, the goal is to compute a complete transport plan r between I and y minimizing the cost <(()) := EsEB ||a (a, b) da. An optimal plan is referred to as a W5 plar I For any parameter s > 0, a transport plan T is called an s-close Op-OT plan e(r) - <(TT() EE, where r" is OT plan.
An optimal solution for the semi-discrete OT problem can be compactly represented as the weighted Voronoi diagram (also called the Laguerre diagram) with respect to a weight assignment y(b) on every point E [2] The choice of weights guarantees that mass at b equal to the mass of Ja inside the Voronoi cell (also called the Laguerre cell) of b and the optimal transport plan simply transports the mass at to the mass inside its Voronoi cell. For 2 dimensions, the weighted Voronoi diagram under the squared Euclidean distance (also known as the power diagram) can be constructed in O(n logn) time 24 57 For higher dimensions d > 2, the construction time would be O ((1[()/2/22) M.
The compact and connected nature of the semi-discrete optimal transport plans makes them attractive for many ML applications; for instance, they can help achieve stability in training GANs and avoid issues such as discontinuities and mode collapse [10 19 555, they improve the mapping between continuous latent spaces and discrete data in Variational Autoencoders (VAEs) [ and also have been used in diffusion- based generative models [42]. However, there are no known exact algorithms for computing semi-discrete optimal transport. Additionally, computing an g-close transport plan is known to be # H--d with respect to d and log 1/e, i.e., an algorithm with an execution time that is polynomial in d and log 1/s seems unlikely [58]. Due to the intractability of the semi-discrete optimal transport in high dimensions, researchers have considered taking n samples from the model, which is continuous distribution Jr, and computing an optimal discrete transport plan between the empirical distribution jh, defined on the n samples and v 433. It has been shown that this empirical p- Wasserstein distance converges to the true semi-discrett p- Wasserstein distance a rate of n 1/12p 22 35 50]; note that the rate of convergence does not depend on the dimension. However, the optimal transport plan from samples is not necessarily a good approximation of the semi-discrete
transport plan and may cause biased gradients [14]]
For the discrete OT problem, there are several near-optimal scalable exact and approximation algorithms I 4 S 20 25 26 36 39 [2 [54], some of which extend to very high dimensions. For the semi-discrete OT in low dimensions, despite extensive work, scalable algorithms to find optimal semi- discrete transport plans remain elusive. There are algorithms to compute an E-close semi- dimi-discrete entropic regularization 66 16 32. and multiscale approaches 44 44 The execution time of all these algorithms is exponential in both and log l/c. Furthermore, their convergence relies on a smoothness parameter of j. For instance, a notable numerical algorithm by Oliker and Prussner 48 assumes that for point d E B, a small change in y vb) will change the mass of μ inside the Voronoi cell of b by a proportionately small amount 45 Remark 22]. Under this assumption, their algorithm executes poly (n, 1/e) iterations, where each iteration requires the computation of a weighted Voronoi diagram which takes n (d)) time. Their algorithm slows down when js is non-smooth and does not even converge when μ is a discrete distribution. Furthermore, these methods approximate the transport cost but the transport plan that they compute may not be an approximation of the optimal weighted Voronoi diagram.
transport plan using numerical solvers [2]
L5
37
Recently, Agarwal et al. [5] described a cost-scaling paradigm to compute an E-close semi-discrete transport plan. Their algorithm executes log scales, where A is the diameter of A U B. Within each scale, they create an instance of the discrete OT problem of size O(n²) in TR2 (and ((d) in IR4). Using any existing strongly polynomial primal-dual discrete OT solver, such as the algorithm by Orlin [49], their algorithm computes an exact discrete OT plan for each instance and updates the weights for B. The OT plan computed in the final scale is an E-close transport plan. They show that the weight assigned to any point b B is away from the optimal weight assignment. The overall execution time of this algorithm is (n log A in 2R22 and n ((d) log A in d-dimensions. Furthermore, in the limit, the transport plan of their algorithm converges to the optimal weighted Voronoi diagram. Note that their algorithm does not make any assumptions on μ.
Our Contributions. The following theorem states our main result.
Theorem 1.1. Let μ be continuous probability distribution defined on compact set A C IR, v discrete probability distribution with a support B C IR2 of size n, and E > a parameter. Suppose
'The p- Wasserstein cost of is defined as e(r)A/P, and the p- Wasserstein distance between μ and is the Wasserstein cost of Op plan.
2

---

there exists an oracle that, given a triangle returns the mass of μ inside yp in & time. Then, an E-close Wz-OT plan between I and v can be computed in On nllog n) log A time in the worst-case, where Å is the diameter of A U B.
Similar to the algorithm by Agarwal et al. [5]. our algorithm is also based on a cost-scaling approach and executes O(log 2) scales. However, the algorithm within each scale is different. Agarwal et al. create a discrete OT instance with n² vertices and n² edges and use an exact discrete OT algorithm to solve this instance. Instead, we extend the combinatorial primal-dual framework of discrete OT to the continuous space and present an algorithm to find the desired cemi-discrete transport plan in O(n* log n time in the worst case. There are several challenges in extending the combinatorial framework to semi-discrete settings, and overcoming these challenges one of the main technical contributions of the paper.
In more detail, in each scale, our algorithm maintains (partial) transport plan, iteratively computes a set of augmenting paths, and augments the transport plan along such paths until all of the mass is transported. In order to assist in finding augmenting paths, we define a residual graph of size ((n)) in the continuous space. Algorithms for discrete OT maintain dual weights for all points in AU B. Unlike in the discrete setting where the vertex set is fixed, the vertex set of our residual graph includes "continuous regions", which evolve over time, and the vertex set of the residual graph changes. Therefore, we are able to maintain weights only for points in B and not for the regions in A. Our primal-dual framework (especially the definition of admissibility in Section 2), as well as our algorithm (in particular Sections 3.2 and 3.4) contain a number of novel ideas, carefully designed to address various challenges that arise due to the dynamically changing continuous regions of the residual graph. Like all existing algorithms, our algorithm also requires access to an oracle that, given a query triangle, returns the mass of μ. inside the triangle. We note that Dijkstra's 's shortest path algorithm has been extended to continuous space (47), but we are not aware of any previous work that extends combinatorial discrete OT framework to continuous space.
Our algorithm extends to any dimension d 2 2 and any p 2 1 in straightforward way. For d >>>> 2, the algorithm in Theorem 1.1 can be shown to have an execution time of OO ( ( log n) log A Note that the runtime of DTiker-Prusseer' 's algorithm has factor 1/E while ours has only log l/e. Unlike their algorithm, ours does not make any assumptions on the smoothness of μ. Furthermore, similar Agarwal al. [5]. our transport plan approximates an optimal weighted Voronoi diagram within a small additive factor.
One consequence of our algorithm is that e-close optimal transport queries efficiently. More precisely, consider a large distribution with support of size N. We can preprocess this distribution into data structure that can return the total weight of points inside query triangle in O(N1-1/d) time [59]. By using Theorem I.1 we can report an c-additive approximate transport plan to any query distribution & with support of size ls in 1111/d N1-1// log ) time. Note that this query time is sub-linear in N. There has been some work on computing 1- I-Wasserstein distance approximately in sub-linear time and for answering nearest neighbor queries under 1- W-Wasserstein distance [3] but these algorithms have larger error and do not extend to 2- Wasserstein distance.
it to design data structure that answers
O(kd+1
Theorem 1.2. Let μ be discrete probability distribution with a support A C R& of size N. The distribution p can be preprocessed, in ON log N) time, into an O(N) size data structure so that for a discrete probability distribution with support C Rdd size k. and parameter e-close W2-OT plan between μ and can be computed in O (Nl 11/d + le klog k) log è time, where A the diameter of A B.
0, an
If the points in the supports of fu and v have integer coordinates and the masses on them are rational numbers, we can adapt our data structure in Theorem 1.2 to compute an exact Wp-OT plan between and v in ((kd+I(N2-1)))) + k log k) log A time for any fixed even value of p.
We also note that in 2-dimensions, for k < N2/4, we can use the data structure to compute an exact discrete OT plan for any distribution v in O(K N1/2 log A) time that is faster than any existing OT algorithm that takes (NK))++(1) time 3 18
For simplicity in presentation, we restrict presentation to d = and p 2.
3

---

h
Figure 1: (left) The 6-expanded Voronoi cell V of b (green shaded area), (middle) the partitioning Xg, and (right) the green region shows the mass of μ that is transported b, the red points show the representative points of regions, and the purple segments show the compressed transport plan i.
Combinatorial Framework
In this section, we extend the combinatorial framework used by discrete OT algorithms to the semi- discrete settings. Let be continuous probability distribution defined over a compact support A C 2 and y be discrete distribution with support set B of points in R
Weighted Voronoi diagram. For any pair of points a, b E IR? let d(a, b) = ||a - blI?. Given a weight function y B + Rzo the points in B, the weighted distance dy (a, b) from point b B to any point E R2 is defined as (a, 6) d(a, b) y(b). For point b B, the Voronoi cell of b the locus of all points with b as their weighted nearest neighbor; more formally,
Vorg (b) = {a E TR2 dy b) & (a, b'), Vb' E B)-
The weighted Voronoi diagram (B) of the points with weights y(-) is the decomposition of IR² induced by Voronoi cells. There exists weight function y(-) for B such that yu( ((Vor,(b) = =(b) for every point b E B, and an optimal semi-discrete transport plan transports the mass of each point E to the mass of j in (or,(b) (b) O A 12
5-expanded Voronoi cell. Consider a weight function y(-) for the points in B. For any point b e B and parameter O, consider the following weight function vg.
[v()) + 5, (1(1')
b b
(b')
V b b.
The 6-expanded Voronoi cell of b. denoted by Vg". is simply the Voronoi cell of in the weighted Voronoi diagram VD,g yy (B) of the point set B with weights (-). See Figure (lef).
o-faasibility. Any (possibly partial) transport plan between y and along with weight function v()) for the points B is 5-feasible
(F1) for any pair (a, b) E Ax B with r(a, b) > the point a lies inside the 26-expanded Voronoi cell V,28.
For any 8-feasible transport plan т, v(-), if is complete transport plan between j and v, then T, y(-)) is called a 6-optimal transport plan. Recall that any optimal transport plan transports the mass at b to its weighted Voronoi cell. In 5-optimal transport plan, however, the mass of each point b B is transported inside a 26-expansion of the Voronoi cell b. This introduces an additive increase of most 20 in the cost of the transport plan.
Lemma 2.1. Any 8-optimal transport plan T, v(()) between p and is 28-close.
See Appendix A for proof. Next, given a d-feasible transport plan T, v()) we define a residual graph and an augmenting path. We also introduce the process of augmenting т along an augmenting path, which allows to increase the mass transported by т.

---

Residual graph. Given 6-feasible (possibly partial) transport plan T, v(-), we construct residual graph G3 by first partitioning the support A of μ into regions to form the vertex set of G3 := ( y. 6) and then defining set of directed edges.
Vertex set. For each point b E B, consider the Voronoi cell of b and its 8- and 26-expansions V and V,25 Let X' denote the arrangement 2 of these 3n cells across all n points of B. See Figure 1 (middle). For each region y in this arrangement, pick an arbitrary representative point T inside yp and assign it mass of pr. == u(p)), where μu = Jgg pf(a) da denotes the mass of μ inside yp Let As denote the set of representative points of all regions in Xg.. The vertex set of Gg is the point set As UB along with source vertex and a sink vertex t. We refer to any point b E B whose mass is not fully transported by т as a free point and define its excess as the amount of mass of that not transported by T, i.e., cx(b) = v(b)- - SacA ( b) da. Similarly, any point rip E AA is free if T does not fully transport the mass into the region sp, and its excess is defined as ex(rx) = Pr, - EvEB T(p, V).
Edge set. For each pair (rp.b)= E As x B, if T(p,, b) > O, we add a backward edge directed from Tip to in the residual graph. Furthermore, if r E V2 we add a forward edge directed from b to rip in Gs. Additionally, we add forward edge from the source to every free point E B and backward edge directed from every free vertex rip to t. This completes the description of the residual graph.
Lemma 2.2. For any > 0 and a d-feasible transport plan f, 2.)( the residual graph G8 has O(n²) vertices and (n²) edges.
While describing our algorithm, is useful to have the definition of weighted distance for all the backward edges, including those incident on t. Therefore, extend the definition of weighted distance to any edge (re, f) as follows. Let by denote the weighted nearest neighbor of rig in B, i.e. bp = minse B dy (rip b). Define dy (rip: t)= == (r bp) 6.
Compressing a semi-discrete transport plan. Given a semi-discrete transport plan 7 that transports mass from B to A. we construct transport plan i from B to As as follows. For each pair (r b) E As x B, let f(r,, b) == T(p, b), i.e., we assign the entire mass transported from b to y to the pair (T.p.) b) (Figure I (riihht). We refer to the transport plan f as the compressed transport plan. Lemma 2.3. For any ô-feasible semi-discrete transport plan T, v(-), the compressed transport plan 7 along with weights v(-) is also 6-feasible.
Conversely, consider a transport plan f. v-)) from B to Ag- One can compute a semi-discrete transport plan that, given pair br E As x B, arbitrarily transports a mass of fr.pe b) from b to y. Lemma 2.4. Any ō-feasible transport plan i, ( B to As can be converted into a 8-feasible semi-discrete transport plan T from B to A.
We say that any compressed transport plan f is a forest if the edges transporting a positive mass in f do not create an undirected cycle.
Augmentation. Given the residual graph G5 for a 8-feasible transport plan f, v-). an alternating path (resp. alternating cycle) is a directed path (resp. directed cycle) in Gs. Note that, in any directed path (resp. cycle) in the residual graph, the edges alternate between forward and backward edges. An augmenting path P- = (s, (1, (22,..... Pk, L bx+1) is directed path from the source vertex s to the sink vertex t in the residual graph. By construction. the vertex b and the vertex TA are free vertices in the residual graph. One can augment along P as follows. Define the bottleneck capacity of the augmenting path P as bc(P) == min min(ex((1), ex(rr),,miniinn..
min,e f(r,, b441). To augment i along P, set f(r,, b.) + f(ri, b,) + bc(P) for each forward edge (b,, ra) E P and (ris, bb+1) + f(r;, b1+1) bc(P) for each backward edge (r:, bi+1) E P.
Lemma 2.5. The transport plan obtained after augmenting a 8-feasible transport plan f, v(-) along any augmenting path P in the residual graph 8 8-feasible.
Consider the following straightforward way to augment a semi- cemi-discrete transport plan. Given a 5-feasible remi-discrete transport plan t, v()), we can compute an augmenting path in the residual graph and augment the compressed transport plan i along P. From Lemmas 2.3 and 2.5 f remains 8- feasible and from Lemma 2.4 the updated 8-feasible transport plan 7 can be converted to 8-feasible semi-discrete transport plan as desired. To obtain complete transport plan, one can iteratively apply this procedure until there are no augmenting paths in the residual graph; however, this may result in an unbounded number of iterations. To obtain efficient algorithm, we iteratively compute set of
5

---

22
b3
5
b1
Figure 2: (left) Three admissible triples (b1, r1, be), (61, r2, b2), and (61, r3, b2), where solid (resp. dashed) lines show backward (resp. forward) edges, and (right) an admissible augmenting path.
special augmenting paths called the "admissible" augmenting paths and augment the transport plan along these paths. By doing so,
reduce the number of iterations to O(n).
Admissibility. Suppose i, x(-) is a d-feasible transport plan between μ and v. For any pair points
61,62 b2 E B and any region & X such that f(rpp 62) 0, the triple (b1, b2) is admissible if dy(r:pbbi bi) < (ry, be I See Figure 2 (left). Note that admissible triple (bi, b2) forms by forward edge followed by a backward edge in the residual graph satisfying dry(r, b1) < dy (Tee by). Intuitively, for any admissible triple (h1, 12), the mass of inside rip is transported from by but bj nearer to Tip than b2 (with respect to the weighted distances).
We extend the definition of admissibility to augmenting paths and alternating cycles as follows. Any augmenting path (resp. alternating cycle) P (b,.r1, be...., Tk. bk+1) is admissible if all triples (bi, rs. 1, E [1, k), are admissible. In Figure riighh,) by transports mass to and by is its weighted nearest neighbor. By augmenting along this admissible path, we increase f(11,b) and reduce f(r,, b2), thereby transporting more mass to rı from its weighted nearest neighbor.
3 Algorithm
In this section, we present our soat-ncaling algorithm that uses the combinatorial framework from Section [ to compute an s-close semi-discrete OT plan. Classical discrete OT algorithms assign weights to points in
and use them identify a large set of augmenting paths. Their efficiency critically relies on the acyclicity of the "search" graph. In contrast, during the execution of our algorithm. a change in the weights of creates new weighted Voronoi diagram, which in turn changes As. the discrete representation of A, and thus the vertex set of the residual graph. Since As may change significantly in each iteration during the exccution our algorithm, we cannot maintain weights for them. This creates significant challenges the algorithm searches for augmenting paths (See Section 3.2 for details). Furthermore, the updated residual graph may have cycles. We introduce additional steps in algorithm eliminate these cycles (See Section 3.4 for details). First, present an overview of our algorithm and then present the details.
3.1 Overview
The algorithm runs for Olog ) scales, where A is the diameter of A U B. In each scale 8, our algorithm maintains transport plan is and weight ybb for every point b B. Initially, set d = 2 and define y(b) = for all E B. each scale 5, execute the following steps.
1. Initialization: Set Tı be an empty transport plan. Compute the residual graph G3 and the compressed transport plan is with respect to Ta. v()).
2. Iterations: While f8 is not complete transport plan:
Discrette OT algorithms use the weights assigned to A U B to define admissible edges. Since As evolves during the execution of our algorithm. we cannot maintain weights for them, forcing to define admissibility sequence of two consecutive edges (bz, r.)) and (r b2) rather than per edge.

---

(i) Compute set of admissible augmenting paths in the residual graph G3 and augment fg along these paths using the SEARCHANDAUGMENT procedure described in Section [.22 At the end of this step, there are no admissible augmenting paths in the residual graph.
(11) Adjust the weights of all points of B that are reachable from the source by admissible paths by 8 and recompute the set As, the residual graph Gs. and the compressed transport plan fs using the ANCREASEWEIGTT
procedure described in Section 3.3
(iii) Update the compressed transport plan and the residual graph using the ACYCLIFY procedure described in Section [.4, so that the transport plan Fg is forest and the residual graph does not have any admissible cycles.
3. Scale Update: Set ā 6/2.
After the execution a scale S e/2, our algorithm terminates by returning a complete semi-discrete transport plan Ta obtained from the compressed complete transport plan fa (using Lemma 2.4)
Invariants. As shown in Section 4 below, our algorithm iteratively updates the weights y(-) and the transport plan ig while maintaining the following invariants:
(I1) The transport plan 4s; v() is 5-feasible, and
(12) At the start of each iteration, the transport plan řg a forest and there are no admissible cycles in the residual graphs
Remark. The algorithm by Agarwal al. [5] creates discrete instance in each scale of the algorithm by computing the arrangement of the 8- 26-..., and (4n + 1)6-expanded Voronoi cells of each point E B. Instead of using such a fine partition to create discrete instance with O(n") edges, work directly with the continuous space, maintain a much smaller residual graph with ((n²) edges,
and use our semi-discrete combinatorial framework to find transport plan.
3.2 SEARCHANDAUGMENT Procedure
The SEARCHANDAUGMENT procedure executes a partial DFS-style search to identify a set of admissible augmenting paths and augments the transport plan along these paths. The SE-RCHAN- DAUGMENT procedure is somewhat similar in style to the blocking flow procedure in Dinic's max-flow algorithm (23 or the partial-DFS procedure in Gabow-Tarjan's algorithm [27), both of which rely on the property that there are no admissible cycles in the residual graph. Unlike these algorithms, in our case, if we augment along an arbitrary admissible augmenting path, we may create an admissible cycle. See Figure 3
We overcome this challenge by carefully calibrating the search algorithm in two ways. First, begin our search from the free regions instead of the free points of B. Thus, we reverse the direction of all the edges of the residual graph and begin our search from the sink l. Second, we explore all forward edges incident on a region in the increasing order of their weighted distance. This order of processing edges ensures that no admissible cycles are created and that there are no more admissible augmenting paths in the residual graph after the SEARCHANDAUGMENN procedure terminates. We provide the details below.
Let Go be the graph formed by reversing the direction of all the edges of G8. We conduct our search starting from the sink in the graph G. Initially, mark all points of B and all backward edges as unvisited, define U := B as the set of unvisited points of B, and Q = = by) as the search path that the procedure grows. Execute the following steps until the search path Q becomes empty.
1. If Q = (t -t=b0 bo, Ti, bi.b......... ...l.......... bi) for some 2 0.
(a) If there is an edge from by to s To (i.c.. by is free point), then P (s. by F)) is an augmenting path in Go. Augment fa along P and set
Q (t bu).
(b) Assume there is not edge from by to 8. If there an unvisited edge (b,, r) in Go then add Ti+1 to Q. Otherwise, mark b, as visited and remove b, from U and Q.
7

---

x
b1
40
P2
40
P22
b2
b4
b3
b4
b3
3
f3
Figure 3: (left) An arbitrary admissible augmenting path (s, b1, ri b2) (blue path), and (right) an admissible cycle (bi, Typa ba, rips ba, ris)) (red cycle) formed after augmentation.
2. Q = bo. F1, bi...., ......,,, bi, Ti+1) for some i 2 O, let b := arg minyeu dy (ri+1,6)) denote the unvisited point with the minimum weighted distance to ri+
(a) If (b, ri41, b,) is admissible, i.e. (ri+1, b) < dy (ri-1, (1t-1,,b),, b,), then add b as b;+1 to Q.
(b) Otherwise, remove 1i+1 from Q and mark the edge (b,, re+1) visited.
The algorithm terminates when the search path Q becomes empty, ie.. the procedure marked as visited and removed it from Q. The following lemma shows the properties of the SEARCHANDAUG- MENT procedure.
Suppose invariants (II) and (12) hold at the start of the SPARCHANDAUGMENT procedure. Then, during the execution of the SEARCHANDAUGMENT procedure,
Lemma 3.1.
(S1) the transport plan 98. v(-) remains 8 ô-feasible,
(S2) any point E (resp. backward edge (r. b)) marked as visited will not form an admissible augmenting path during the execution of the procedure, and
(S3) there admissible cycles in the residual graph.
3.3 ExCReASEWEIGHES Procedure
The MCCEASSWEIIGHT procedure adjusts the weights of B leading to the formation of admissible augmenting paths. The CERREASEWEIGHTS procedure performs DFS on the residual graph Go to compute a set IC of all vertices of B that are reachable from the source vertex by admissible paths. It increases the weights of all points in K. by 5 (expand their Voronoi cells). Recall that the set As. and thus Go as well as the compressed transport plan îs depend on the weights of B. The procedure then recomputes As. Gs. and is from 7g. See Appendix B.2 for complete details of the SCRRASSEEIGHNN procedure. The following lemma states the properties of the INCREASEWS SEREASSWEIIHHS procedure.
Lemma 3.2. Suppose invariant (II) holds at the start of the INCREASEWEIGHTT procedure. Then, during the execution of the INCREASEW EIGHTS procedure,
(WI) the transport plan 78. v(.) remain 5-feasüble,
(W2) the weight of each free point b B increases by 5, and
(W3) the weight of each point e with fiee regions inside Vi remains unchanged.
3.4 ACYCLIFY Procedure
The change in the weights of B in the INCREASEWEIGHE procedure requires us to recompute the residual graph and the compressed transport plan. This recomputation may potentially create cycle in the transport plan or an admissible cycle in Gj. The ACYCLIFY procedure eliminates such cycles and ensures that the invariant (12) holds at the start of the next iteration. Converting the procedure
transport plan into forest is critical for the efficiency of the while eliminating admissible cycles is essential for the correctness of the SEARCHANDAUGMENT procedure. The procedure executes the following steps. First, use the dynamic tree structure by Sleator and Tarjan (56] to make the transport plan fg forest. Then, execute a partial DFS search from
"To perform this step efficiently, in the construction of the residual graph, for each E As, we store the list of all neighbors of sorted in increasing order of their weighted distance.

---

each unvisited point b E B similar to the one described in the SEARCHANDAUGGEENT procedure to detect admissible cycles. Upon finding an admissible cycle, cancel the cycle right away. remove the vertices of the cycle from the search path, and continue the search. When all vertices are visited, no admissible cycles are remaining the residual graph. Since canceling admissible cycles could have introduced new cycles in the transport plan, repeat the first step to update is and make it forest. See Appendix [.3 for details. The following lemma shows the properties of the ACYCLIFY procedure.
Lemma 3.3. Suppose invariant (II) holds at the start of the ACYCLIFY procedure. Then, during the execution of the ACYCLIFY procedure,
(AI) the transport plan fs: v(-) remains 5-feasible, and
(A2) the transport plan 98 is a forest and there are no admissible cycles in the residual graph.
This completes the description of our algorithm.
4 Algorithm Analysis
In this section, we first prove the correctness of our algorithm and then analyze its runtime.
Proof of invariants (II) and (I2). For any scale 8, the initial transport plan is is empty. Therefore, ig along with the weights y() is 5-feasible. By properties (S1), (W1), and (A1), the transport plan fs. v(-) remains 8 f-faasible in each iteration of our algorithm, and therefore, invariant (II) holds. The invariant (12) is a direct consequence of property (A2) in Lemma 3.3
Proof of Correctness. From Invariant (II), in each scale 6, our algorithm maintains a d-feasible transport plan y-) during its execution. The while loop in Step 2 breaks when is is a complete transport plan. Therefore, igg along with weights y(-) is 8-optimal. From Lemma 2.4 one can convert is into a 8-optimal semi-discrete transport plan Ty- Given that our algorithm terminates when 5 - s/2, from Lemma 2.1 the transport plan returned by our algorithm is E-close, as desired.
Efficiency of the algorithm. The SEARCHANDAUGMENT procedure runs a partial DFS on the residual graph with O(n²) edges. Upon finding an augmenting path P. the procedure augments the transport plan along P in O( time. In Lemma C.I we use invariant (I2) to show that the total length of all augmenting paths found by the procedure is O(n²). Hence, the running time of the procedure is Ó(n²). The INCREASEWEICE EIGHTS procedure stores a sorted list of neighbors for each region and executes a DFS procedure in the residual graph, which takes O( log n) time. After increasing the weight of a subset of points of B, the procedure recomputes the residual graph G5 and the transport plan is. which can be done in O(n ( + n log n)) time. See Lemma C.2 in the appendix for details. In the ACYCLIFY procedure, converting a transport plan into a forest can be done using dynamic tree data structure [56] in On² log n)) time. To eliminate the admissible cycles, the ACYCLIFY procedure relies on a partial DFS that runs in ((n)) time. Therefore, the total time by the ACYCLIFY procedure is ( log n) (see Lemma C.3 for details). Combining the running times of all three procedures, each iteration of step 2 of our algorithm takes Onn2 1- + n log n.) time. In the following lemma, we show that in each scale, the total number of iterations of step 2 of our algorithm is at most (n)).
REARCHANDAUGMENN
Lemma 4.1. For each scale ö, the total number of iterations of step 2 of our algorithm is O(n).
Proof Sketch. Let 726, (-)(-) denote the 25-feasible semi-discrete transport plan computed by our algorithm for scale 25, and let Ta, (-) denote a partial semi-discrete transport plan maintained during the execution of step of our algorithm. Let X2s (resp. X's) denote the partitioning of the set A with respect to weights yes (-) (resp- V)(-)). Let )/ be the arrangement of all 3n cells used to construct X28 with all 3n cells used to construct Xs. Let P22 (resp. fs) denote the compressed transport plan for 725 (resp. Tτ) using the partitioning V. It is well-known that one can transform fs to î25 by augmenting is along a set of augmenting paths P on ) x B and rearrange the transported mass along a set of cycles C on y' x B. Consider an augmenting path P = (r1, bi...., rk. by) in P. Since P is a simple path, it contains each point of at most once and therefore, it has length at most 2n 1. Additionally, for all i E [1, k], F26 (ra b,) > O and for each i E [2, k), (r.. b;-1) > O. For each edge (ra, bz), since (r., b,) > 0, by 26-feasibility of 726, (-2(-) (condition (F1)), the point by is a 45-weighted nearest neighbor of r; with respect to weights y2s( (-). Similarly, for each edge (b-,,, ri). since is( b--1 > O. by 8- d-feasibility of T8, (n(-) (condition (F1)), the point b is a 26-whidt nearest neighbor of r, with respect to weights yo (.). Since the length of P is at most 2n and the

---

weight of by does not change (by property (W3) in Lemma 3.3). we show that for the free point by in P. (bx)) V2s (b.] 6nô. Our algorithm increases the weight of by by 8 in each iteration (property (W2) in I Lemma 3.3 and therefore after 6n iterations, the point by cannot be free, i.e., after (n)
iterations, all points of B are fully transported in ÿs. We provide the full proof in Section C.4
Using Lemma 4.1 the total time spent in step 2 in each scale of our algorithm is On 4 + n log n)). Since there are (log ) scales, the overall runtime is ( + n log n) log A ), thereby proving Theorem
5 Applications to the Discrete Optimal Transport Problem
In this section, we extend our combinatorial semi-discrete OT algorithm to the discrete OT problem and design a data structure that preprocesses and stores large discrete distribution and efficiently computes an E-close OT plan between js and any query distribution support size of μ. More precisely, given discrete distribution μ with a (possibly large) support A of N points in IR?, we design a data structure that, given a query discrete distribution v with a support B of points, computes an e-close transport plan between ll and v in (k- + log k) log è time. Additionally, we show that if the support points have bounded integer coordinates and the masses rational numbers, our data structure can efficiently compute an exact discrete OT plan.
sub-linear time relative the
At a high level, our data structure interprets the large discrete distribution and uses simplex ange-searching data structure as an oracle to compute the mass of μ inside a query triangle. In this way, for any query distribution V, one can execute the steps of our combinatorial semi mini-discrete algorithm to compute an E-close transport plan between μ and v. More formally, our data structure preprocesses the distribution into a simplex range-searching data structure RS-DS, which takes O(N) space, can be built in O(N log N) time, and returns the mass of μ inside a query triangle in = O(VN) time [59]. Given query discrete distribution V, one can use our algorithm from Section in conjunction with the RS-DS to compute an e-close transport plan between μ and in O(k* (VN & log k) log è) time leading to Theorem 1.2
continuous distribution
Consider the special case where the points in A B have positive integer coordinates bounded by the mass of j (resp. v) on each point a E A (resp. b E B) is rational number of form Es (resp. 45) for positive integers and Ta (resp. It), and p is an even number. In this case, the p- Wasserstein cost of any transport plan between fı and v/ is an integer multiple of 7 and therefore, any 27 -close transport plan between μ and v would have a minimum cost. Thus, one can compute an exact discrete OT plan between μ and by setting
1
2T in our data structure, which would have query time of O Olk + le log k) log(XT)), leading to the following corollary.
Corollary 5.1. When the points in the supports of the distributions μ and have integer coordinates bounded by X and the mass on each point is a rational number of form 4, our data structure computes, for any even number p. an optimal solution for the p- Wassersteii problem between p and in O (VN + le log k) log(XT)) time.
6 Conclusion
In this paper, we designed a novel combinatorial framework for the semi-discrete optimal transport problem and used it compute E-close semi-discrete transport plan. We also used this framework to design a data structure that stores a discrete distribution Jı over a large support of size N and can compute e-close OT cost between μ and query discrete distribution v in time that sub-linear in N. We conclude with the following question: Can we use our combinatorial framework to compute an s-close semi-discrete transport plan between high dimensional distributions in poly(n, 1/e) time.
Acknowledgement
Work by P.A. and K.Y. was supported by NSF grants CCF-22-23870 and IIS-24-02823 and by BSF grant 2022131. Work by S.R. and P.S. was supported by NSF grant CCF-2223871. We thank the anonymous reviewers for their useful comments.
10

---

References
[1] Pankaj K. Agarwal and R. Sharathkumar. Approximation algorithms for bipartite matching with metric and geometric costs. In Proc. 46th Annual ACM Sympos. on Theory of Computing. STOC 14, pages 555-564, 2014.
[2] Pankaj Agarwal and Micha Sharir. Arrangements and their applications. In Handbook of Computational Geometry, pages 1191 Elsevier, 2000.
[3] Pankaj Agarwal, Kyle Fox, Debmalya Panigrahi, Kasturi R Varadarajan, and Allen Xiao. Faster algorithms for the geometric transportation problem. In 33rd Internat. Symposium Computational Geometry, 2017.
[4] Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Rachita Sowle. A higher precision algorithm for computing the 1- 1-Wasserstein distance. In Proc. 11th Internat. Conference on Learning Representations, 2023.
[5] Pankaj K Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Keegan Yao. Fast and accurate approximations of the optimal transport in semi-discrete and discrete settings. In Proc. 2024 Annual ACM-SIAM Sympos. on Discrete Algorithms, pages 4514-4529. SIAM, 2024.
[6] Jason M Altschuler, Jonathan Hiles-Weed and Austin Stromme. Asymptotics for semidiscrete entropic optimal transport. SIAM Journal of Mathematical Analysis, 54(2):118--7741 2022.
[7] Luca Ambrogioni, Umut Guclu, and Marcel van Gerven. Wasserstein variational gradient descent: From semi-discrete optimal transport to ensemble variational inference. arXiv preprint 181 2018.
[8] Luca Ambrogioni, Umut Güçlü, and Marcel van Gerven. k-gans: Ensemble of generative models with semi-discrete optimal transport. arXiv preprint arXiv:1907.04050 2019.
[9] Dongsheng An, Yang Guo, Na Lei, Zhongxuan Luo, Shing- Tung Yau, and Xianfeng Gu. AE-OT: new generative model based on extended semi-discrete optimal transport. In Proc. 32nd Internat. Conference Learning Representations, 2019.
[10] Dongsheng An, Yang Guo, Min Zhang, Xin Qi, Na Lei, and Xianfang Gu. AE-OT-GAN: Training gans from data specific latent distribution. In Proc. 16th European Conference on Computer Vision, pages 548-564. Springer, 2020.
[11] Franz Aurenhammer. Power diagrams: properties, algorithms and applications. SIAM journal on computing, 16(1):78-9))) -6,,9,,,11 1987.
[12] Franz Aurenhammer, Friedrich Hoffmann, and Boris Aronov. Minkowski-yy type theorems and least-squares clustering. Algorithmica, 20(1):61-76, 1998.
[13] Arturs Backurs, Yihe Dong, Piotr Indyk, Ilya Razenshteyn, and Tal Wagner. Scalable nearest neighbor search for optimal transport. In Internat. Conference Machine Learning, pages 497-506. PMLR, 2020.
[14] Mare G Bellemare, Ivo Danihelka, Will Dabney, Shakir Mohamed, Balaji Lakshminarayanan, Stephan Hoyer, and Rémi Munos. The Cramer distance as a solution to biased Wasserstein gradients. arXiv preprint arXiv: 10743, 2017.
[15] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the Kantorvccto mass transfer problem. Numerische Mathematik, 84(3):375-393, 2000.
[16] Bernard Bercu and Jérémie Bigot. Asymptotic distribution and convergence rates of stochastic algorithms for entropic optimal transportation between probability measures. Annals of Statistics, 49(2):968-98) 2021.
[17] Rick Chartrand, Brendt Wohlberg, Kevin Vixie, and Erik Bollt. A gradient descent solution to the MongeeKantoroovia problem. Applied Math. Sci., 3(22): 2009.
11

---

[18] Li Chen, Rasmus Kyng, Yang P Liu, Richard Peng, Maximilian Probst Gutenberg, and Sushant Sachdeva. Maximum flow and minimum- cost flow almost-linear time. Proc. 63rd IEEE Annual Sympos. Foundations of Computer Science, pages 612-623. IEFF, 2022.
[19] Yucheng Chen, Matus Telgarsky, Chao Zhang, Bolton Bailey, Daniel Hsu, and Jian Peng. A gradual, semi-discrete approach to generative network training via explicit Wasserstein minimization. In Internat. Conference on Machine Learning, pages 1071-1080. PMLR, 2019.
[20] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, pages 2292-2300, 2013.
[21] Fernando De Goes, Katherine Breeden, Victor Ostromoukhov, and Mathieu Desbrun. Blue noise through optimal transport. ACM Transactions of Graphics (1(6):1-11, 2012.
[22] Eustasio Del Barrio, Alberto González Sanz, and Jean-Michel Loubes. Central limit theorems for semi-discrete Wasserstein distances. Bernoulli, 30(1):554-580 2024.
[23] Efim A Dinic. Algorithm for solution of a problem of maximum flow in networks with power estimation. In Soviet Math. Doklady, volume 11, pages 1277-1280, 1970.
[24] Steven Fortune. Voronoi diagrams and Delaunay triangulations. In Handbook of Discrete and Computational Geometry, pages 225-265, 1995.
[25] Emily Fox and Jiashuai Lu. deterministic near-linear time approximation scheme for geomet- ric transportation. In Proc. 64th IEEE Annual Sympos. on Foundations of Computer Science, pages 1301-1315. IEEE, 2023.
[26] Kyle Fox and Jiashuai Lu. A near-linear time approximation scheme for geometric transportation with arbitrary supplies and spread. In Proc. 36th Annual Sympos. on Computational Geometry, pages 45:1-45:18, 2020.
[27] H. N. Gabow and R. E. Tarjan. Faster scaling algorithms for network problems. SIAM Journal of Computing, 1815)10013 October 1989. ISSN 0097-5397.
[28] Bruno Galerne, Arthur Leclaire, and Julien Rabin. Semi-discrete optimal transport in patch space for enriching Gaussian textures. In Proc. 3rd Internat. Conference on Geometric Science of Information, pages 100-108, 2017.
[29] Bruno Galerne, Arthur Leclaire, and Julien Rabin. texture synthesis model based on semi- discrete optimal transport in patch space. SIAM Journal of Imaging Sciences, 11(4):2456-2493, 2018.
[30] Thomas O Gallouët and Quentin Mérigot. A Lagrangian scheme la Brenier for the in- compressible Euler equations. Foundations of Computational Mathematics, 18(4):835-865, 2018.
[31] Jing Gao, Peng Li, Asif Ali Laghari, Gautam Srivastava, Thippa Reddy Gadekallu, Sidra Abbas, and Jianing Zhang. Incomplete multiview clustering via semidiscrete optimal transport for multimedia data mining in loT. ACM Transactions of Multimedia Computing, Communications and Applications, 20(6):1-20 2024.
[32] Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach. Stochastic optimization for large-scale optimal transport. Advances in Neural Information Processing Systems, 29, 2016.
[33] Antoine Houdard, Arthur Leclaire, Nicolas Papadakis, and Julien Rabin. A generative model for texture synthesis based on optimal transport between feature distributions. Journal of Mathematical Imaging and Vision, 65(1):4-25 2023.
[34] Martin Huesmann and Khro-Th Sturm. Optimal transport from Lebesgue to poisson. The Annals of Probability, 41(4))2426-2478 2013.
[35] Shayan Hundrieser, Thomas Staudt, and Axel Munk. Empirical optimal transport between different measures adapts lower complexity. In Annales de I'Insituu Henri Poincare (B) Probabilites et statistiques, volume 60, pages 824-846, 2024.
12

---

[36] Andrey Boris Khesin, Aleksandar Nikolov, and Dmitry Paramonov. Preconditioning for the geometric transportation problem. arXiv preprint NXiv:19022.8884,, 2019.
[37] Jun Kitagawa. An iterative scheme for solving the optimal transportation problem. Calculus of Variations and Purtial Difserential Equations, 51(1):243-263, 2014.
[38] Jun Kitagawa. Quentin Mérigot, and Boris Thibert. Convergence of Newton algorithm for semi-discrete optimal transport. Journal of European Math. Sociery, 2109):603-2651, 2019.
[39] Nathaniel Lahn, Deepika Mulchandani, and Sharath Raghvendra. A graph theoretic additive approximation of optimal transport. Advances Neural Information Processing Systems, 32, 2019.
[40] Arthur Leclaire and Julien Rabin. A stochastic multi-layer algorithm for semi-discrete optimal transport with applications to texture synthesis and style transfer. Journal of Mathematical Imaging and Vision, 63(2):882-308, 2021.
[41] Bruno Lévy and Erica I. Schwindt. Notions of optimal transport theory and how implement them on computer. Computers & Graphics, 72:135-145 2018.
[42] Zezeng Li, Shenghao Li. Zhanpeng Wang, Na Lei, Zhongxuan Luo, and David Xianfeng Gu. DPM-OT: a new diffusion probabilistic model based on optimal transport. In Proc. IEEE/CVF Internat. Conference on Computer Vision, pages 22624-22633, 2023.
[43] Huidong Liu, GU Xianfeng, and Dimitris Samaras. A two-step computation of the exact GAN Wasserstein distance. In Internat. Conference on Machine Learning, pages 3159-3168, 2018.
[44] Quentin Mérigot. A multiscale approach to optimal transport. Computer Graphics Forum, 30 (5):1533-1592, 2011.
[45] Quentin Merigot and Boris Thibert. Optimal transport: discretization and algorithms. In Handbook of Numerical Analysis, volume 22, pages 133-212. Elsevier, 2021.
[46] Jocelyn Meyron. Initialization procedures for discrete and semi-discrete optimal transport. Computer-Aided Design, 115:13-22, 2019.
[47] Joseph SB Mitchell, David M Mount, and Christos H Papadimitriou. The discrete geodesic problem. SIAM Journal of Computing, 16(4):647-668, 1987.
[48] Vladimir Oliker and Laird D Prussner. On the numerical solution of the equation
8 Uzóy
= and its discretizations, i. Numerische Mathematik, 54(3):271-293, 1989.
[49] James Orlin. A faster strongly polynomial minimum cost flow algorithm. In Proc. 20th Annual ACM Sympos. Theory of Computing, pages 377-387, 1988.
[50] Aram- Alexadde Pooladian, Vincent Divol, and Jonathan Niles-Weed Minimax estimation of discontinuous optimal transport maps: The semi-discrete case. In Internat. Conference on Machine Learning, pages 28128-28152 150. PMLR, 2023.
[51] Hongxing Qin, Yi Chen, Jinlong He, and Baoquan Chen. Wasserstein blue noise sampling. ACM Transactions of Graphics, 36(5):1-13, 2017.
[52] Sharath Raghvendra and Pankaj K Agarwal. A near-linear time e-approximation algorithm for geometric bipartite matching. Journal of the ACM, 6793):1-19, 2020.
[53] Nathan A Seifert, Kirill Prozument, and Michael J Davis. Computational optimal transport for molecular spectra: The semi-discrete case. The Journal of Chemical Physics, 156(13), 2022.
[54] R Sharathkumar and Pankaj K Agarwal. Algorithms for the transportation problem in geometric settings. In Proc. 23rd Annual ACM-SIAM Sympos. on Discrete Algorithms, pages 306-317, 2012.
13

---

[55] Ilia Shumailov, Zakhar Shumaylov, Yiren Zhao. Nicolas Papernot, Ross Anderson, and Yarin Gal. Ai models collapse when trained on recursively generated data. Nature, 631(2)
755-759, 2024.
[56] Daniel D Sleator and Robert Endre Tarjan. A data structure for dynamic trees. In Proc. 13th Annual ACM Sympos. on Theory of Computing, pages 114-122, 1981.
[57] Kokichi Sugihara. Laguerre Voronoi diagram on the sphere. Journal for Geometry and Graphics, 6(1):69 81, 2002.
[58] Bahar Taşkesen, Soroosh Bhaficezadeh-Abadee and Daniel Kuhn. Semi-discrete optimal transport: Hardness, regularization and numerical solution. Mathematical Programming, 199 (1):1033-1106, 2023.
[59] Emo Welzl. Partition trees for triangle counting and other range searching problems. In Proc.
4th Annual Sympos. Computational Geometry, pages 23-33, 1988.
14

---

A Missing Proofs of Section 2
In this section, we present the missing proofs and details of our combinatorial framework described in Section 2
Lemma 2.1. Any ô-optimal transport plan T, v(-) between p and v is 2ô-close.
Proof. Define the weighted cost of a transport plan 7' as [y(7") = ElEB SA (y(a, da. For any complete transport plan T',
(T( Σ dg(a,
Σ (d(a, - 9(6)) = <(7)) Σ ((6)((b).
(2)
For any point a E A, let ba denote the weighted nearest neighbor of a, i.e. ba == arg minse dy (a, (). By property (F1) in the definition of 5-feasibility,
(T) = Σ B
(a, (a, b) da S
Σ
(dy(a, bo) + ( b) da
DEB
Σ
(a) b) da 25.
(3)
BEB
Let 7* denote optimal transport plan from p Then,
c,( = Σ dy (a," da > Σ E dy(a, ((((,bb)T"(aa da.
(4)
B
BEB
Combining Equations and 44 and plugging and 7" in Equation
<(T) (T) v(b) Σ / da da 28 + Eu(b)u(b)
BEB
S (T) + +
Σ v(b))(b) = <(r") 28.
Therefore, <(7) S <(T*) + 26 and r is 26-close.
Residual graph. Next, we show that for any 5 > O and any d-feasible transport plan f, v(). the residual graph Ga has On??) vertices and O edges, leading to Lemma 2.2 To do so, we show below that the partitioning X's consists of Onn²) regions. We then conclude that the number of vertices of Gs is Onn²). Furthermore, since the residual graph is a bipartite graph between set B of size n and set As of size O(n²), the number of edges of Go would be at most O(1²).
Recall that the partitioning Xs is constructed the arrangement of all Voronoi cells, 5-expansions, and 2o-expansions of the Voronoi cells of all points in B. Let V denote the set of vertices of this arrangement. Since the arrangement is a planar graph, the number of faces (i.c., regions) in X is ( Therefore, to show that X's has Onn²) regions, we show that the number of vertices of this arrangement is On?²).
For each point b E B, let Vi (resp. V V:2 denote the Voronoi cell (resp. d-expanded Voronoi cell, 25-expanded Voronoi cell) of b, and let y (resp. vi, yo denote the weights B used to compute the cell. Note that each Voronoi cell V for each E B and E {0, 28} has at most n vertices. Hence, the total number of Voronoi vertices the arrangement X5 is O Next, we count the number of intersection points of these Voronoi cells. For each pair of points bi, b2 E B, consider the pair of cells V and V2 for some i.j E {0, 6, 26}. We show that Vi and V intersect each other in at most two points.
Define the weighted bisector of two points b and b/ with respect to weights y(.) as the locus of points that have the same weighted distance to and b' i..., all points r RR² such that (x b) dy (z. B"). Note that under the squared Euclidean distances, the weighted bisector of two points is a straight line. Let denote intersection point of Vž, and Vž and suppose v lies on the segment of Vg
15

---

22
v
b
Figure 4: For two cells Vi (green cell) and Vi (orange cell), any intersection point lies on their weighted bisector.
corresponding to the weighted bisector of b and a point by E B (See Figure 4). Note that if by = b2. then Vi and V/2 share a segment containing which means that
dy bi) - == =dgi, (0.b1) dug, (v, b2) = (1(v,
and also
dy b2) - dvi (v, b2) dzi (v, = (v
Therefore, in this case, j number Voronoi vertices of the arrangement. Hence, we assume by & b2. In this case,
j== = O, and the two endpoints of the shared segment is counted toward the
dg. (v. b1) dsi, (v. ba) min dii ().
(5)
BEB
By the construction of the weight function vi, (-) in Equation I for any point b E B (b),b2). vi, (b) = whg (b) = y(b;; thus, using Equation (5) d 1/2 (v. ba) = dgi, (v. by) min Be
dgi (u, = min Be b2)
d (v. b).
(6)
Using Equation 66. the point also lies in the segment of V2 that corresponds to the weighted bisector of b2 and by with respect to weights vbg. Thus,
dgg (v, bi) doi, (v, ba) = d (v, ba) d (v,
In other words, the point satisfies
d(v,bi) - (v(b1) + = = (2((22) +
ie., the point lies on the weighted bisector of b and by (the blue dashed line in Figure 4. Since the weighted bisector a straight line, it intersects the convex polygon Vř, in at most points; hence, the two Voronoi cells V2, and V intersect each other in most points. Since the O pairs of such Voronoi cells, each intersecting each other in at most 2 points, the total number of intersection points is at most O(n²), claimed.
Lemma 2.2. For any > 0 and ô-feasible transport plan f, v.))), the residual graph Gs has ( vertices and (n?²) edges.
Next, we extend our bounds for the number of vertices and edges of the residual graph in Lemma 2.2 to any dimension 2 and show that for any pair d-dimensional distributions f and and any 5 5-feasible transport plan f. (-), the residual graph G& has O(n²) vertices and Oz edges. We prove bounds by showing that the number of vertices of the arrangement X's is O(n)) Using this
16

---

bound, we get that the number of regions in the arrangement is upper-bounded by O(n"); therefore, |Asl= = (n²)) and the residual graph has ((n)) vertices. Furthermore, since 98 is a bipartite graph between set B with n points and set Ag with (nn)) points, the number of edges would be bounded by O(nd+1).
For each point b B and each i {0, ō, 26), the Voronoi cell V has O(n (d22) ) vertices; since 3n such cells are used in the construction of the arrangement, the total number of Voronoi vertices in the arrangement is (n [4/2]+1). ). Next, we bound the number of intersection points of these Voronoi cells. For each d-tuple (V)) Voy,..,, Vid of Voronoi cells, for d distinct points bi.-.., bu and d values i1,...,4 ia E {0, 6, 28). let e denote a point in the intersection of these cells. Similar to our analysis for 2-dimensional distributions, we show that e lies on the weighted bisector of bi.... ba. where the weight of each point bj is y(b;) + ij. Since this weighted bisector is a straight line, and since a straight line intersects a convex polytope in at most 2 points, each d-tuple of Voronoi cells intersect in at most 2 points, where the number of such d- d-tuples are O(n Therefore, we conclude that the number of intersection points is Onn)). as desired.
.-...
Let be an intersection point of the d Voronoi cells Voz V/2..... and Vi. Suppose v lies on the (d - -)-dimensional hyperplane of V that is shared between b1 and point ba+1 E B. In this case,
d Ubg (u, b1) = d i (v, bat1) min a (v, be Ub
We claim that for any l: E (1, .......,. d]. yor (ba+1) = y(bat1). Consider the following two cases:
If b841 = by for some j E 12,...,.................... d}. then
dg(0, - i1 = d (u.b1) = =d 96j
(v. by) dy(r, bj),
and also
dy(w:bs) - in = 0/
(1) = d i, b1) = dy(u;
Therefore, i = ig = 0 and for any A E {1....,d}.
d}. 1 (but1) = v(ba+1).
Otherwise, b141 $ (bi..... ......,bl ba) and by the construction of yhr (.), for any h E (1,...,d].
vit (ba+1) = y(ba+1).
Therefore, from Equation (7)
dy, ba+1) = d (v. ba+1) ) = min d (u, b) = -imming dy b), be US1
(8)
9b1
BEB
and for any k E {l....,dd d}. since v lies on a (d - diddinnnsional hyperplane of Vida it has to lie on the hyperplane of Vidg that is shared between by and bat1; therefore, d(u, ba+1) y(ba+1) d(u, bb) - (b(bx) + it). Thus,
= From Equation ( the point v lies on the weighted bisector of b,,..., ba with weights y(bi)
d(v,bi) (v(b1) ) + +i) d(v,b2)) (1(b2) + 12) ... =...= d(v, ba) (y(ba) + ia)-
(9)
..,..., .bbba) + dd and the d-tuple (V% Vig. ...., Vild) intersect each other in at most points.
Lemma A.1. For any dimension d 2 2, a parameter 8 O. and ô-feasible transport plan 4, y(-),
the residual graph Gs has O(n²)) vertices and (nd+1) edges.
Compressing the transport plan. Recall that for any semi-discrete transport plan 7 from B to A, of weights v(-), and parameter > O. the compressed transport plan è is discrete transport plan from B to Ag. Note that the definition of d-feasibility (and more precisely, condition (F1)) is applicable for the discrete OT as well; the discrete transport plan along with weights y(-) for B is 6-feasible if f transports the mass of each point b E B to points of As within its 26-expanded Voronoi cell. Any complete transport plan i from B3 to As that is ô-feasible along with weights y(-) is called 6-optimal.
Here, the constant is O(n hides a factor of 3d
17

---

Lemma A.2. Any 8-optimal compressed transport plan i, y-)) B to As is 28-close.
Proof. For any transport plan r' from B to Ag define (,y = E(rD)) dylr
by(, (r, b). Then,
(T(T))=
E (nD)EAsBB
(r(,,,))(r,,b))
Σ
b) - ((()))))))(((bbb) =
Σ (6))((b) BEB
(rblEAsx
(10)
For any point E As, b, denote the weighted nearest neighbor of if., b, == arg m d, b). By the definition of 6-feasibility, for any pair (r, b) with f(r, b) > O. we have dr, S dg 4 +28. Therefore,
(x)
E
b)(r,))
Σ
(dg(1;br) + 2
(r.))AAX
(r.b)EAsxB
E
d,(r,
(,),,,
+ +
(11)
(r,b)EAxBB
Let 7" denote an optimal transport plan from to Ag. Then,
e (T(
E
b)
Σ
br)(,,))r"
(12)
(n.b)EAssP BB
(5x)
Combining Equations 1I
and 12 and plugging i and 7" in Equation 10
<(f) = (t) + (b))
Σ
+ 28
+ Σ vb6)r(6)
BEB
(x
BEB
S (Ty + 25 +
Σ
v(b)u(b) = <((7() + 26.
beB
Therefore, <(f) - <(T") + 25 and is 28-close.
The following observation, which is straightforward from the construction of the partitioning Xs. helps in proving Lemmas 2.3 and 2.4
Observation A.3. For any region yp E X5 and any point E B, the region ip either completely lies inside v/2 it is completely outside V226.
Lemma 2.3. For any ô-feasible semi-discrete transport plan r, v.)) the compressed transport plan along with weights y(-) also 8-feasible.
Proof. Consider any point b E B and any region y E Xs such that ((p, > 0, i... in the compressed graph, f(rp. b) > To prove this lemma, we show that rip V/24. Since 7(p. b) > there exists point a E y such that (a, b) > and by the 6-feasibility of 7, y-)), we have E V,2 therefore. using Observation A.3 the region yp has to completely lie inside V,25, and the representative point Tip- which is point inside y also lies inside V,2 Hence, ł, y(-) is d-feasible.
Lemma 2.4. Any ò-feasible transport plan f, v.-)) to As can be converted into a 8-feasible semi-discrete transport plan from B to A.
Proof. Consider any transport plan from B to A such that for any point b E B and any region E Xs. transports a mass of fr b) from to the mass of μ inside the region sp One such construction is to assign T(a,b) ((r..)) u(a) for each point sp. We next show that the transport x(Φ)
plan 6-feasible.
Consider any point E B and any point a E A such that T(a, b) > 0. We prove this lemma by showing that a E V,26. Suppose
Xs is the region containing the point a. Since T(a, b) > 0, we should have Fr b) > 0, and by the d-feasibility o y.), we have r E V226, By Observation A.3 the whole region yp has lies inside 1/26 and therefore, the point also lies inside V Hence, v-) is d-feasible.
18

---

Augmentation. In this section, we first show that upon augmenting transport plan f along an augmenting path P. the transport plan remains valid. We then show that upon augmentation, either backward edge gets removed from the transport plan or free point will become fully transported. We finally prove Lemma 2.5
For an augmenting path P = (s, Tk. from a free point b to a free point Tk. recall that
tk.t)
the bottleneck capacity of P iss defined as
be( P) == min(ex(bi), ex(ra), (x)(()))))x()))))))) min f(ri, b.41)}
vE(1,8-1]
where ex(b1) (resp, ex(rx)) denotes the excess mass b1 (resp. rι). To show that the transport plan after augmentation is a valid one, we show that no point in G8 transports more mass than the mass at and that each edge transports a non-negative amount of mass.
In the augmentation process, for any i {l...., k - 1}, we decrease fr,, b(41) by be( (o(P), where by definition, be bc(P) S f(r,, be41); hence, F(r,, b11) 2 after augmentation. Furthermore, for any forward edge (i), ( i kl, we increases Fr:;, b,) by be( (e)P) and f(r;, b.) remains non-negative.
For any i E (2...., k.. we increase F(6,, r4) (resp. decrease f(r;-1, b.)) by bc(P); hence, the total amount of mass transported from b, remains unchanged. Similarly, for each E {1, K 1}. we increase F(b;, rz) (resp. decrease frri, bi+1)) 1) by be( and the total amount of mass transported into rs remains unchanged. For the endpoint bj (resp. rg), we only increase the amount of mass transported from bj (resp. into r4:) by bc(P), where by definition, bc(P) I cx(bı) (resp. be(P) I ex(rx)). Therefore, the total mass transport from b1 (resp. into r3) after augmentation is at most v(b) (resp. Mrg), as desired.
Note that if backward edge br; P determines the bottleneck capacity of the augmenting path P. then Frr b) = after augmentation and the backward edge removed from the residual graph. Otherwise, if the endpoint bı (reps. r*) determines the bottleneck capacity of P, then bj (resp. rς) will be fully transported after augmentation.
Lemma 2.5. The transport plan obtained after augmenting a 8-feasible transport plan 7, v(-) along any augmenting path P in the residual graph is 5-feasible.
Proof. Let P (s,b), Fi,.... bb, Tk: 1) denote an augmenting path in the residual graph. When augmenting the transport plan along P, we increase the mass transportation forward edges (b,,r,) r.) for each E [1. ks] and decrease the mass transportation on the backward edges (r,, bx+1) for each i [1, 1]. Therefore, any pair (r, b) that transports mass after augmentation but was not transporting mass before augmentation has to be a forward edge of P. Since we only add forward edges from the point to the representative points in V,26 n.Ag the edge (r, b) satisfies the d-feasibility condition (F1), and f, y(.) remains feasible after augmentation.
B Missing Details of Section 3
In this section, we provide the missing details of the implementation of the INCREASEWEIGHTS and ACYCLIFY procedures and also prove the properties of the three procedures.
B.1 Missing Proofs of the SEARCHANDAUGMENT Procedure
Lemma 3.1. Suppose invariants (II) and (12) hold at the start of the SPARCHANDAUGMENT procedure. Then, during the execution of the SEARCHANDAUGMENT procedure,
(SI) the transport plan f6. v(-) remains ō-feasible,
(S2) any point b E B (resp. buckward edge (r, b)) marked as visited will not form admissible augmenting path during the execution of the procedure, and
(S3) there are no admissible cycles in the residual graph.
Proof. We prove the properties separately in the following-
Property (SI). By the construction of the search path, any augmenting path computed by the SEARCHANDAUGMENT procedure is admissible augmenting path. From Lemma 2.5 augmenting
19

---

b4
b4
b
b2
b5
b5
b6
b6
b1
b1
b3
D3
Figure 5: (left) An augmenting path found by the SEARCHANDAUGMEN procedure (blue path), and (right) two admissible triples (highlighted in pink) formed after augmentation.
is along admissible augmenting path does not violate the 6-feasibility condition; hence, the trans- port plan fs. v() is a d-feasible transport plan during the exccution of the SLARCHANDAUGMENT procedure and (S1) holds.
Next, we present an overview of a new property of the SEARGHANDAUGMENT procedure, which we formally state and prove in Lemma below. Using that, we first prove (S3) and then prove (S2).
In Figure 5 (left), suppose the blue edges show an admissible augmenting path found by the ANGAAMENAT procedure, and suppose the green (resp. red) points show the unvisited (resp. visited) points of B. For the region yp, suppose (b1.b2,... b2...... b6) denote the set of neighbors of To in Go. sorted in increasing order of their weighted distance to T.p. When the partial DFS procedure processes ri., the two weighted nearest neighbors of oo (i.c., 61 and b2) are already marked as visited, leading the procedure to add by to the search path. After augmentation (Figure 5 (right)), for the newly created backward edge (r.g,3)) the only admissible triples containing (r.p 03) are (bi,rg. and (bg, rpp by) (the triples highlighted in pink), where by is unvisited and both b2 and by are visited. More formally, as shown in Lemma [B.]] for an augmenting path P found by the procedure, assuming that (S3) holds before augmentation along P. for any newly created admissible triples (b, Tip. 1/) after augmenting the transport plan along P. the point b (resp. b') is marked as visited (resp. unvisited).
Property (S3). We begin by presenting an overview of our proof. Consider any augmenting path P found by the SEARCHANDAUGMEEN procedure. Assuming that (S3) holds before augmentation along P. all vertices that are reachable from a visited point by an admissible path in the procedure are also visited, since those points were also added to the search path, did not lead an admissible augmenting path, marked as visited and removed from the search path. Hence, all points having an admissible path to the visited points (i.e., all points that are reachable from the visited points in our backward DFS) are also visited. Therefore, there are no admissible paths from an unvisited point to a visited point in the residual graph. After augmenting along P, by Lemma B.T I for any newly formed admissible triple (b, To B/) the point b (resp. 6') is visited (resp. unvisited), and since there are no admissible paths from any unvisited point to any visited point, the newly formed admissible triple (b, r.p. 6') does not form cycle of admissible triples. Hence, (S3) holds after augmentation well. We provide the details of the proof below.
Let pl...., pk denote the sequence of augmenting paths computed by the denote the sequence of transport plans computed by the procedure,
procedure, and let +8, $ ....... #
i.e., 7g is the transport plan computed in the previous iteration, and for each [1, k). f3 is obtained by augmenting + along p Let Go denote the residual graph corresponding f3 for each i [0, k).
Let V (resp. U³) denote the set of visited (resp. unvisited) points in B when the procedure augments - along ps.
Initially, from invariant (12), there are no cycles of admissible triples in the residual graph, and (S3) holds for go. For any E k), assuming "....., Gi-1 satisfies (S3), we show that G also satisfies (S3). Suppose (b, r, b') is any admissible triple in Gi formed after augmenting -i-1 along pi, i.e., the triple (b, b') is admissible in G but not in gi-1. We show that the triple (b, r, 1/) does not participate any admissible cycles; hence, using property (S3) on gi-1, there are no admissible eycles in Go and property (S3) holds for G as well.
20

---

Since (S3) holds for .".....9--,, for any visited point b E Vi-1, the point b has been added to the search path Q by the SEARCHANDAUGMENT procedure, did not lead to the computation of an augmenting path, marked as visited and removed from the path. In this case, any point 8" that
is reachable from b by an admissible path in g' (and therefore, is reachable by our partial DFS procedure from b) would have also been added to the path, marked as visited and removed from Q. i... any point 6" E B that has an admissible path to the visited point b in the residual graph G is also visited. By this observation, there are no admissible paths from any unvisited point in B to any visited point in Gi From Lemma B.1 for any newly formed admissible triple (b, r, V'), we have b E v and '' E - Thus, all admissible triples formed after augmenting Fg p/-1 along P are from a visited point to an unvisited point, while there are no admissible paths from any unvisited point to any visited point; therefore, the newly formed admissible triples do not form any admissible cycles and (S3) holds for G as well.
Property (S2). We use the property (S3) to show that (S2) holds. For any point b E B that is marked as visited, as discussed above, (S3) holds, all vertices that are reachable from b by an admissible path in our backward DFS (i.e., all points having an admissible path to b in the residual graph) are also visited. Since any free point by E B is unvisited, by does not have an admissible path any visited point b E B and therefore, the visited points do not participate in an admissible augmenting path. Furthermore, the procedure marks a backward edge (r, b) as visited if, for each admissible triple (b', r, b). the point b' is visited. Since the point b' cannot be included in an admissible augmenting path, the visited backward edge (r, b) also does not form an admissible augmenting path.
Lemma B.1. During the execution of the AEDAGGENN procedure, suppose P is an admissible augmenting path found by the procedure, and let G (resp. G') denote the residual graph before (resp. after) augmenting the transport plan along P. Let (b, r, b)) denole an admissible triple in G' that is not admissible in G. Assuming that (S3) holds prior augmentation along P, the point b is marked as visited, and b/ is marked as unvisited.
Proof. Consider any triple (b, r, b') that is admissible in G' but not in G. Recall that by the definition
of the admissible triples, (r, b') is a backward edge in G' and d, (r, > (r, Since the SEARCHANDAUGMEEN procedure does not change the weights y(-), the only case where (b, r. b') is not admissible in 9 is when (r. 8') is not a backward edge in G, i.c., the pair (b', r) is in P as a forward edge, and augmenting fs along P results in transporting mass from b' to r. On the other hand, by step 2(a) of the SEARCHANDAUGMENT to the search path only if D' is the weighted nearest unvisited neighbor of r: in other words, since dy b) > (r, b') and the procedure added D' to the search path (instead of b), the point was marked as visited by the procedure. Therefore, for any newly formed admissible triple (b, r)), point b (resp. 8') is marked as visited (resp. unvisited).
procedure, a forward edge (6)r) will be added
B.2 Missing Details and Proofs of the INCREASE W WEIGHTS Procedure
After the execution of the SEARCHAND remain the residual graph, i.e., there are no admissible paths from the source vertex to the sink vertex L. The INCREASEW WERRETS procedure increases the weights of the subset of points in B that are reachable from & by admissible paths to expand their Voronoi cells and to create new admissible triples in the residual graph. For instance, in Figure (left), the path from 8 to L is not admissible (the triple (b2, T2. ba) is not admissible, as by has a lower weighted distance to r2 than b2). As shown in Figure 6 (right), the SCREASSWEIGHE
procedure, no admissible augmenting paths
procedure then increases the weight of the points 6 and be (which are reachable from s). leading to the formation of an admissible augmenting path (note that upon updating the weights, the regions corresponding 71 and r2 have slightly changed). The details of the INCREASE WEIGHTS procedure are described below.
For each point E As. let (rrr S B denote the set of points b E B with (a(r.b) > O, sorted in decreasing order of their weighted distance to r, i.e., (rr) = (bI,..., ..... bx) where (r, bn) > (r, b2) > ... .d (r: bz). Mark all points b E B and all forward edges (b.r) as unvisited and set K 0. let U = denote the set of unvisited points of B, and define Q := (s) as the search path that the algorithm grows. Execute the following steps until Q becomes empty.
21

---

2
b3
22
b1
Figure 6: (left) After the execution of the SEARCHAANAUGMEER procedure, there are admissible augmenting paths in G5. and (right) by increasing the weights of the points that are reachable from by augmenting paths (points bj and b2), new admissible triples are created (e.g. (62, r2, b2)), which might lead to the formation of admissible augmenting paths.
Q = (s), then if there exists an unvisited point b e U such that (s, b) E Go. then add to Q as b1. Otherwise, remove from Q.
2. If Q = (s,
bi) for some 2 1.
(a) there exists an unvisited forward edge (b,,r) in Gs. add r to Q as T++1-
(b) Otherwise, mark b; as visited, remove b; from U, add b; to K. and remove 6, from Q.
3. If Q (s, bi. .... by, Tx) Fi........
let b arg mingeurw(r)) dy 6') denote the unvisited point of b with the minimum weighted distance to among all points of B that transport mass to r.
for some 1,
(a) If (b,, Ti, b) is admissible, i.c., dy > (ra, b.), then add b as bi+1 to Q.
(b) Otherwise, remove r, from Q and mark (b,, r;)) visited.
After the DFS procedure terminates, for each point b E K, set y(b) + y(b) + 5. This completes the description of the DFS step. We next describe how recompute the residual graph and the compressed transport plan with respect to the updated weights.
Let v(-) (resp. y' (-) denote the weights of the points in B after (resp. before) the weight updates, let Xs (resp. X²) denote the partitioning of the set with respect to weights y(-) (resp. y((.), and let As (resp. A²)) denote the of representative points of the regions in X's (resp. 2))) Furthermore, let igg denote the transport plan maintained by the algorithm for partitioning Xz To compute the new transport plan řş for the point set As, the INCREASEW EIGHTSS procedure first computes the arrangement y of all 3n cells used to construct X% with all 3n cells used to construct Xs. i.e., ) is the arrangement of Voronoi cell, 8 5-expanded Voronoi cell, and 25 26-expanded Voronoi cell of each point b E B both before and after weight updates. See Figure 7 For each region C(p) S y denote the set of regions of y that lie inside jp For each region e E Y, pick an arbitrary representative point rg inside g. We denote the set of all representative points of the regions in y by Y.
X' U X's let
The EVCRAASEEEIGHT procedure first converts i/ to a transport plan ł over the finer partitioning Y x B by simply splitting each region y' in X (and its mass transportation) to fine regions of y inside o.. The procedure then uses f to construct transport plan ta over As x B by merging the regions of Y' (and their mass transportation) to regions of X'. The details are provided next.
For each point E B, each region vp' E Xy with fg(rrbb)) b) > O, and each d' E C(p'), set F(e', , b) m(P) #f b). This completes the description of the split step and the construction of A. Next, the procedure constructs is by setting, for each point b E B and each region yp E Xs. fT b) Epec(e) f(e.b) The transport plan ts is defined over As B.
22

---

Figure 7: (left) The partitioning X before updating weights, (middle) the partitioning X after updating the weights of the points of that are reachable from s by admissible paths (the blue cross points), and (right) the combined partitioning y.
Finally, for each region E Xg. our algorithm stores a list (rrp) of all points h E B with Thp E V2225 sorted in increasing order of their weights distance to Tip, i.e., (rrp)) ((r)== (bi...., bu), Thp E V22 for each i E [1, ki). and dy (rip: bb) S (r for each S < S k. This completes the description of the INCREASE WEIEEWEIGHTT procedure.
Lemma 3.2. Suppose invariant (II) holds at the start of the ANCREASEWEIGHT p Then, during the execution of the INCREASEW SEEEASSEWIGGTSS procedure,
(WI) the transport plan is. v()) remain 6-feasible,
(W2) the weight of each free point b E B increases by 6, and
(W3) the weight of each point E B with free regions inside Vo remains unchanged.
Proof. Let y(-) (resp. y'() denote the weights of the points after (resp. before) the execution the INCREASEW WEGIITS procedure, and let X (resp. X²) be the partitioning with respect to weights y(-) (resp. s (-). For any point a E A, let Sa (resp. (()) be the region in X (resp. X3) that contains a, and let ba (resp. by)) denote the weighted nearest neighbor of a in B with respect weights y(.) (resp. y').). For any pair of points (a, b) E A x with (a, b) > O, by the f-fassibility T/ ( we have dy (a, b) 25 < dy (a, b²). To prove property (W 1), show that dy(a, b) 26 S dy bu).
Recall that the INCREASE W FIGHTS procedure finds the subset K. C B of points that have admissible paths from the source vertex of the residual graph and increases the weights of all points in K. by 8. Consider the following cases:
i E IC is among the points whose weights are increased by the procedure, then
(a b) (a, b) - < dy (a, (/)) + < dy (a, ba) + 8 S dl (,) + 25,
where the second inequality holds from 5-feasibility of T/S, 3-). the third inequality holds since by the weighted nearest neighbor of with respect to y (-), and the last inequality holds since y(h.)) y (bbb) + Consequently, dy b) 25 S (a,
Otherwise, b && K. and (a, (a,
If by E K. then dy (r...; b) < dr ba) (since otherwise, the triple (ba, Tg.. b) would have been admissible triple and b, IC would have resulted in b E K)). In this case, since the weighted nearest neighhor of and r.p. are the same,
(T b.)) I dy (Te., b) < dy (T.p., b.) S dy (N..: 5,.) + 8,
where the last inequality holds since increasing the weight of V/2 by 5 made it the weighted nearest neighbor of a. Therefore, the region Sa and consequently, the point lie inside V and V Thus,
dy b) = dy (a, b) s (a ba) + dy bu) + 28.
Otherwise, by is also not in K., then (a ha) dy (a, ba) and,
da b) 28 (a, 28 d (a, b.) dl bo).
23

---

bi
ba
bb
Figure 8: An example of cycle in transport plan that is created while increasing weights: (left) transport plan that is forest, (middle) the new Voronoi diagram and partitioning (red dashed lines) after increasing the weight of by, b2, and by, and (right) cycle formed in the new residual graph.
As a result, the transport plan Ty along with the updated weights y(-) is 8-feasible and property (W1) holds.
Note that each free point by has an edge from the source vertex s; therefore, by E and the INCREASEWEIGHTT increases the weight of by by 8, proving (W2). Furthermore, for each point b B with free region yp Vgg since defined (r miny e dy (r 8), the triple (b, r.p.t) t) is admissible. Therefore, the SEARCHANDAUGOE GMENT procedure should have added re and b to the search path, not found admissible augmenting path, and marked b (resp. (r, t)) visited.
Therefore, there are no admissible paths from the source to b, or equivalently b & IC, and the weight of b remains unchanged, leading to (W3).
B.3 Missing Details of the ACYCLIFY Procedure
The goal of the ACYCLIFY procedure is to ensure that the invariant (I2) holds, i.c., that at the beginning of each iteration of our algorithm, the transport plan is forest and there are admissible cycles in the residual graph. The procedure runs three steps: (1) make îs a forest, as described in Section B.3.1 (2) cancel any admissible cycles from the residual graph, described in Section B33.22 and (3) acyclify the transport plan again, as described in Section B.3.1 Note that our procedure acyclifies the transport plan (to make it a forest) twice, in steps (1) and (3). Making the transport plan a forest in the first step is essential for the efficiency of the second step, and making it forest in the third step is essential for invariant (12), canceling admissible cycles might introduce cycles in the transport plan.
B.3.1 Acyclifying the Transport Plan
Similar to the Acyclify procedure introduced in [54] Section 3.3], we use a dynamic tree structure to make Fa forest as follows. Let È = (e1, e2..,e...... Ea) denote the set of all edges
(r, b) E As x B with 4s b) > O. For any S u, let E == (e1, e2. ,e)) Define Fo @ as an empty forest and (r, b) = 0 for all pairs (r, b) E Aa x B. Starting from k = 1, for any & < the algorithm computes forest F& and transport plan i defined over F₁ using Fx-1 and k-1 as follows. If
adding the edge ex to Fk-1 does not create cycle, then the algorithm simply sets F3 + Ffx) (ex) + fe(ex) and (e) + Th-1(E) ) for all edges E Fx Otherwise, adding e to Fk-1 results the creation of an even-length cycle C. Let e denote the minimum capacity of the edges in and let e" denote the edge with the minimum capacity. Consider an ordering of the edges of the cycle that starts with e", i.e.,
(e* = C/, Ezz). The algorithm increases (resp. reduces) the mass transported along the edge E'2 (resp. 1 by e for each E [1,j]. Finally, the algorithm sets FR + Fa-1 U [ex} \ ( }. This completes the description of step 1. Using the dynamic tree structure by Sleator and Tarjan [56), each operation takes Ollog n) amortized time, and since |8] = O(n²),
this process takes total of On log n) time.
24

---

D3
by
122
b/2
bb
Figure 9: An example of the formation of admissible cycles due to updating the weights and the residual graph: (left) a transport plan with no admissible cycles in the corresponding residual graph, (middle) the Voronoi diagram and partitioning (red lines) after increasing the weights of by and b2, and (right) an admissible cycle (bi, Fi, b2, r2, b3, ra) formed in the new residual graph.
B.3.2 Acyclifying the Admissible Triples
To remove all admissible cycles from the residual graph, we use partial DFS similar the one described in the direction of all the edges of Gs. The procedure first marks all points of and all backward edges as unvisited and defines U := B as the set of unvisited points. While there exists an unvisited point b E B, the procedure initializes partial DFS by setting Q (b b1) and searches as follows until Q becomes empty.
procedure. Let So be the graph formed by reversing the
1. If Q (bI. b;) for some 1 1.
(a) there are no unvisited backward edges (b,, r) in To. then mark b, as visited and remove by from and U.
(b) Otherwise, there exists unvisited backward edge (b,, rp). Add ni as T++1-
2. If Q = (b,, (1..,, bt....,,,,,
rı) for some i 2 1, let b == arg mingeu dy (r;, 1)) be the unvisited point with the minimum weighted distance rj.
(a) If (b, rv, b.) is admissible, i.c., dy (r., b) < (r,r,, bi),
If b already exists in the path Q as 6j. then C' = (by,, jj...., bis Ti, b;+1 = = b,) is an admissible cycle. Cancel the cycle C described below) and set Q (b1.71...... ..r.............. bz).
Otherwise, add as by+1 to Q.
(b) Otherwise, remove r; from Q and mark the backward edge (ri, b:) as visited.
Given an admissible cycle C' = (by, b,, bi+1 = by), the procedure cancels the cycle as follows. Let bo(C)) ;= min fa(re, 6,) denote the bottleneck capacity of the cycle C. For any forward edge (br+1, P2) (resp. backward edge (r,, bc)) in C, set i (rr, b++1) ia (ri, bet1) + be(C) (resp. (Fr, b,) < (r,, b,) be(C)). In this way, at least one of the backward edges of the cycle Cs removed from the residual graph and the cycle has vanished.
Lemma 3.3. Suppose invariant (II) holds at the start of the ACYCLIFY procedure. Then, during the execution of the ACYCLIFY procedure,
(Al) the transport plan 68, v(-) remains ò-feasible, and
(A2) the transport plan 9, is a forest and there are no admissible cycles the residual graph.
Proof. In the first step of the ACYCLIFY procedure, the algorithm makes the transportation network to be a DAG. Note that the resulting transportation network is a subset of the transportation network before the ACYCLIFY procedure, and therefore. for any transporting edge (r. As x B, the point b is weighttt nearest neighbor r, i.e., the transport plan obtained after the first step of the ACYCLIFY procedure is d-feasible.
Next, we show property (A1) and (A2) for the second step of the procedure. By construction, all triples the search path maintained by the procedure are admissible, and therefore, any cycle
25

---

computed by the second step of ACYCLIFY procedure is admissible. Since all forward edges on the computed cycles are from points E to regions in V,20 canceling cycle C' does not violate the d-feasibility condition (FI); hence, the transport plan Ts, v(.) is d-feasible transport plan during the execution of the second step of the ACYCLIFY procedure and (A1) holds.
To prove property (A2), we show that
(A3) any point E B (resp. backward edge (r, 6)) marked as visited will not form admissible cycle during the execution of the procedure, and
(A4) during the execution of the second step of the ACYCLIFY procedure, the subgraph of the residual graph induced by visited vertices and their neighboring regions does not have any cycles of admissible triples. Furthermore, there are no admissible paths from an unvisited point to a visited point.
Assuming property (A3) holds, any point E B that is marked as visited does not participate in an admissible cycle. Since the ACYCLIFY procedure stops when all points in B are visited, there are no admissible cycles in the residual graph. Furthermore, since the transport plan is maintained by dynamic tree structure the transport plan is forest. Therefore, prove property (A2), we show that (A3) holds.
Furthermore, the property (A3) is direct corollary of property (A4), explained next: For each visited point E B, all points 8' E B having admissible paths to are also visited (note that the ACYCLIFY procedure searches on the residual graph in the reverse direction of the edges), and if (A4) holds, there are no admissible cycles solely formed by visited vertices; hence, b does not form admissible cycles. Furthermore, the procedure marks backward edge (r, b) visited if, for each admissible triple (b', r, b). the point b' is visited. Since b' is not part of any admissible cycles (assuming (A4) holds), the triple (b', r, b) also cannot be a part of an admissible cycle. Therefore, to prove property (A3), we prove that (A4) holds in the following.
We use inductive argument prove (A4). Let C',..., denote the sequence of admissible cycles found by the procedure, and let ff...., f3 denote the sequence of transport plans, where #g is the transport plan maintained by the algorithm at the beginning of the second step of the procedure and fJ is obtained by canceling f3 along C². The property (A4) trivially holds at the beginning of the execution of the second step of the ACYCLIFY procedure.
C*
Note that the ACYCLIFY procedure marks point b E B as visited only if the search from b did not lead to finding an admissible cycle, i.e., for all pairs (b/r)) E B x As such that (b', r, b) is admissible, the point b' is marked as visited and the backward edge (r, b) is marked as visited. Hence, all points having an admissible path to the visited points (i.e., all points that are reachable from the visited points in our backward DFS) are also visited. Therefore, Assuming that property (A4) holds before marking as visited, property (A4) holds after marking as visited as well. Next, we show the same when we cancel cycle C
First, note that all points of B that are on the search path are unvisited. Therefore, for any admissible cycle found by the procedure, the points of B on the cycle are unvisited. Using an identical proof as in Lemma B.I one can show that after canceling a cycle, for any newly formed admissible triple (b, r, b¹), the point b is visited and the point // is unvisited. In other words, canceling a cycle C" only creates additional admissible paths from visited points to unvisited points. Assuming that (A4) holds before canceling C³, since no new admissible triples are created from visited point to another visited point, the subgraph induced by visited points and their neighboring regions remain free of admissible cycles after canceling C'é. Furthermore, since new admissible triples are created from an unvisited point to a visited point upon canceling C³, there will be no admissible paths from an unvisited point to visited point after cancellation; hence, (A4) remains satisfied.
C Missing Details of Section 4
In this section, we analyze the efficiency of the three procedures SEARCHANDAUGMENTT IN- CRIGSEWEIGHTS, and ACYCLIFY.
26

---

C.1 Efficiency of the SEARCHANDAUGMENT Procedure
The SEARCHANNAUGGENT procedure runs a partial DFS on the residual graph to find a set of admissible augmenting paths. The partial DFS procedure, upon backtracking from a point b E B (resp. 7 E As). marks the point b (resp. the backward edge (b',r) used to reach r) as visited and does not visit the point 6 (resp. the backward edge (b', r)) again in the same exccution. Upon finding an augmenting path P. the procedure augments the transport plan along P in O((PI) time. Let procedure. (P...... Ps) denote the set of all augmenting paths found by the SEARCHANDAUGMER In 2 dimensions (resp. d dimensions), since the residual graph has O(n²) edges (resp. O(nd+1)), &
the running time of the procedure would be O(n² + - (PP)) (resp. O( + + (P.))). In Lemma C.I we show that Ei=1 PP = O(n²) (resp. Σi=1 | = =((n&t1)). Hence, each execution procedure takes ((n)) (resp. O(nd+1)) time.
of the SEARCHANDAUGMER
Lemma C.1. The total length of augmenting paths computed during the execution of the SEARCHAN- DAUGMENT procedure is ((n)) in dimensions and ( in d dimensions.
Proof. Let ff denote the transport plan maintained by the algorithm at the beginning of execution of procedure. To prove this lemma, we categorize the augmenting paths found by the procedure based on the source of their bottleneck capacity, namely (1) set P. consisting of augmenting paths whose bottleneck capacity is determined based on the residual capacity of its endpoints, and (2) set Pe consisting of augmenting paths whose bottleneck capacity is determined based on mass transportation over its backward edges. We first show that |Pul = Ón?) and then show the same bound for Pe. Since each augmenting path has length of at most 2n, we then conclude that the total length of all augmenting paths is O(n²).
the
Let P be augmenting path in P,.. If the bottleneck capacity of P is determined by a free point b E B (resp. free region r E As). then the mass of b (resp. r) will be fully transported after augmentation; therefore, since |As U B O(n"), we have iPol On??). Next, let P be an augmenting path in Pe: in this case, the backward edge (r, b) determining the bottleneck capacity of P will be removed from the transport plan after augmentation. Note that by Lemma B. I for any newly formed admissible triples (V/ r, b). the point is already marked as visited. By property (S2) in Lemma 3.1 the point B/ cannot form admissible augmenting path during the same execution of the
procedure. Hence, the edge (r, b) determining the bottleneck capacity of P was backward edge of the initial transport plan 7 and augmentation along each path / E Pe removes one of the transporting edges of the transport plan fg. Using invariant (I2), |Pel O( as claimed. Hence, the total number of augmenting paths found by the procedure is O(n)), and since each augmenting path has a length of at most 2n, their total length O(n²).
Next, we extend our analysis to d dimensional space, for any d 2 2. Note that in d dimensions, the residual graph has O(n²) vertices and O( edges. Hence, PP = O(n"). Since the transport plan #8 is a forest over the point set As U B, the total number of edges transporting a positive mass in f8 would be Onn"); since augmenting the transport plan along each augmenting path in Pe climinates one of the edges transporting positive mass in fg. |Pel = On")). Finally, since each augmenting path has length of at most 2n, the total length of all augmenting paths found by the SEARCIAANDAUGMENT procedure would be O(nd+1).
C.2 Efficiency of the INCREASEWEIGHTT Procedure
In this section, we show that for 2-dimensional (resp. d-dimensional) distributions, the IN- CREASEWEIGITTS procedure runs in O(n² + n l)) (resp. On + n log n))) time. The INCREASEWEIGITT procedure runs a DFS that visits each edge of the residual graph at most once and has a total rumning time of O(n²) (resp. O(r+1))) )). Furthermore, in the arrangement used to construct partitioning ), each point b E B has at most 6 Voronoi cells (three cells that are used in the construction of X's and three that are used in the construction of X²). Using a similar discussion as Section A one can show that the total number of vertices in the arrangement used to construct y is O(n² (resp. O(n)),, and the number of regions in y at most O(n?) (resp. O(n")). The construction of the transport plan & can be done in O ( + n)) (resp. ( + +))) time since (1) the mass of all regions in y can be determined in O(n?d) (resp. O(n))) time (partitioning the regions in y into simplices remains arrangement with Onz²) ) (resp. O(r))) vertices and therefore,
27

---

have (n²)) (resp. )(n")) regions), and (2) the mass transported on each pair (e. 6) ) x B can be determined in O(1) time. Converting 5 to ig. as is done in the merge step, also takes O (resp. O(n*t)) time, given that the total complexity of is ( (resp. O(n*11). Finally, storing a sorted list of neighbors for each region As takes On log n) (resp. O(1=~1) log n) time in total. Hence, the execution of the SCCRAASWEIGGHT procedure takes Onn" + n log n)) (resp. ( + n log n)))) time.
Lemma C.2. Each execution of the INCREASEWEIGHT procedure takes On ($ + n log time in dimensions and On ( + logn)) time in d dimensions.
C.3 Efficiency of the ACYCLIFY Procedure
The first step of this procedure uses dynamic tree structure to acyclify the transport plan fa- Using the dynamic tree structure by Sleator and Tarjan [56], since the total number of edges of the graph is (n²) in 2 dimensions (resp. (n in d dimensions), the running time of this step would be (n² log n) (resp. O(n log n)). In the second step, the procedure runs partial DFS procedure on the residual graph and cancels the admissible cycles. The partial DFS procedure, upon backtracking from a point b B (resp. r E As), marks the point b (resp. the backward edge (b',r) used to reach r) as visited and does not visit the point (resp. the edge (b',r) again in the same execution. Furthermore, upon finding an admissible cycle C, it cancels the cycle in O(CC) time. Let (CI...... Ck) denote the set of all cycles found in the execution of the ACYCLIFY procedure. In Lemma C.4 we show that Σ
|Cd (n²²) in 2 dimensions (resp. Σx= |Ci| = O(n in d dimensions). Given that the size of the residual graph is at most (nn)) (resp. ((d the second step of the ACYCLIFY procedure takes a total of O Σ (Cal) (n²) (resp. (( time, leading to the following lemma.
Lemma C.3. Each execution of the ACYCLIFY procedure takes On logn) time in 2 dimensions and (n d+1 logn), time in d dimensions.
Lemma C.4. The total length of admissible cycles computed during the execution of the second step of the ACYCLIFY procedure is O(r²)) in 2 dimensions and ((nd+)) time in dimensions.
Proof. Let f denote the transport plan maintained by the algorithm at the beginning of execution of the ACYCLIFY procedure. To prove this lemma, we show that the ACYCLIFY procedure finds (n²²) admissible cycles, where each cycle has a length of at most 2n; hence, the total length of all cycles found by the procedure would be (n²).
Let C be an admissible cycle found by the procedure; in this case, the backward edge determining the bottleneck capacity of C' will be removed from the transport plan after cancellation. For any admissible triple (b, r, b') formed after canceling C, using an identical argument as Lemma B. one can show that the edge (r. B)) is backward edge that was on the cycle C as a forward edge, and the point b is marked as visited: hence, by Lemma 3.3] the point does not form an admissible cycle in the same execution of the ACYCLIFY procedure and therefore, the newly formed backward edge (r, b') cannot be included in any admissible cycles. Therefore, each cycle cancellation removes one of the backward edges of 7gg which is the transport plan obtained after the first step of the ACYCLIFY procedure, i.e.. the transportation network of iff is a forest and the number of its transporting edges is Onn²). Therefore, the total number of cycles found by the ACYCLIFY procedure is O(n²), and their total length is O(n²). as claimed.
We next show that the total length of admissible eycles in d dimensions is O(n in d dimensions, for any d. 2 2. Note that in dimensions, the residual graph has ((14+1) edges. Since the transport plan # is a forest over the point set As B, the total number of edges transporting a positive mass in +g would be O(n"); since canceling each admissible cycle eliminates one of the edges transporting positive mass in :f. |C| ( Since each admissible cycle has length of most 2n, the total length of all admissible cycles found by the INCREASEW EIGHTS procedure would be ((d+11). Li
C.4 Number of Iterations
Lemma 4.1. For each scale ö, the total number of iterations of step of our algorithm is ((n).
28

---

Proof. Let 728, (26(-) denote the 28-febee transport plan computed by our algorithm for scale 26, and let Ts- ya (-) denote a partial transport plan maintained during the execution of step 2 of our algorithm. Let X28 (resp. X2) denote the partitioning of the set A with respect to weights V-s(.) (resp. v=(-)) Let y be the arrangement of all 3n cells used to construct X2s with all 3n cells used to construct Xs. i.e., 3/ is the arrangement of Voronoi cell, 8-expanded Voronoi cell, and 26-expanded Voronoi cell of each point b B with respect to weights ys(+) along with the Voronoi cell, 29-expanded Voronoi cell, and 46-expanded Voronoi cell of each point b E B with respect to weights y-s(.). For each region @ E V pick an arbitrary representative point Te inside p. We denote the set of all representative points of the regions in J by Y. Let f25 (resp. 58) denote the compressed transport plan for 726 (resp. Tτ) using the partitioning V Note that the partitioning y is a refinement of both partitionings X₃ and X25. Define T' := f28 fg.
We construct a bipartite graph G' over Y x B, where for any pair (r, b) E Y x B, there exists an edge directed from to b if (r, b) < 0 and directed from b to if (r, b) >> 0. Consider any directed path = (r1, b.....: Tk. by) from free region r E to a free point b E (with respect to ř&).
The path / is augmenting path in the residual graph corresponding to řs, ys(.).
Similar to the standard graph algorithms, we define the net-cost of the path P as ((P) E d(rk, by) E d(rk, 1. Let by = b,, be the weighted nearest neighbor of T1. Then, can rewrite the net-cost of P as
((P) a(r1, 0n) +
[d(re, b.) - -drs, be-1)
= dys (r1,
[dge (ri, b.) dye (r., b b(-1)] + vb(b)).
Due to d-feasibility of the transport plan fs. (-). for all edges (ri, 1) i E [1, k)., the point b;-1 is a 26-weighted nearest neighbor of r;; hence, dys (Ta, bi) - dhys (ri, b;-b);1)) 2 22.. Since the length of P is at most 2n- 1,
((P) = dys (r1, bp) +
dg (ri, (rrrbb)-che - dya (rss b)] + +16(bb) 2 dys (r1, bo) + vb((b)) 2nd. (13)
Similarly, we can rewrite the net-cost of P using weights y/25 as follows.
of (P) = dlyss (r1, bo) + Σ ddys (ri, b.) days (ris 1)] + 126(bx).
According to the 26-feasibility of 428, (2)(.), for all edges (rs, b:), E [1, k), the point by is a 48- weighted nearest neighbor of rx (with respect to weights yas (-)); hence, dyyo (ré, by) cyys (ri, b;-1) S 48. Since the length of P is at most 2n 1,
((P) dyzz (TI, +
ddgs (Ta, bb) dyys (ris bi-1) + y2s(bb) I dyys (r1, bo) + (bs(ba) + 4nd.
(14)
3.3 during the execution of step 2 of our algorithm, the weights of the points in B with free regions in their d-expansion remain unchanged; therefore, since bn contains free regions inside its d-expanded Voronoi cell, y25(bn) ybs(b)) and dzys (ri,bo) dyy (ri, bo). Combining with Equations (13) and (I4p,
By property (W3) in Lemma
(b(b)) 2nd s (PP) dys (r1,bo) = o(P) - dyza bo) S (bx + 4nd.
Equivalently,
(bb(b)) - y2s(bb)) S Gnd.
By property (W2) in Lemma 3.3
algorithm increases the weight of the point by by 8 in each iteration while its mass is not fully transported by tx; therefore, the point by cannot remain free after 6n iterations, i.e., after (n)) iterations, there cannot be any remaining free points in B and the step of our algorithm terminates after (ni)) iterations.
29

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: This paper presents a new combinatorial framework for the semi-discrete OT problem. It used design algorithm for computing e-close semi-discrete transport plan, as well as a data structure for the discrete OT problem. Our main results are stated in [.] and [.2 in the introduction and are summarized in the abstract.
Theorems
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Our semi-discred OT algorithm relies the existence of oracle that computes the mass of μ inside a triangle, which is common in existing algorithms for the semi-discrete optimal transport problem. We clearly state this assumption in the abstract and Theorem
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: All theorems and lemmas contain the full set of assumptions, and their proof's are either provided in the main text or the appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [NA]
Justification: This paper does not include experiments.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [NA]
Justification: This paper does not include experiments.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [NA]
Justification: This paper does not include experiments.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [NA]
Justification: This paper does not include experiments.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
30

---

Answer: [NA]
Justification: This paper does not include experiments.
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: There are no data or human subjects involved in the research process of this paper and there no known societal impacts and potentially harmful consequences associated with the research conducted in this paper.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
None of the categories deseribed by the NeurIPS Code of Ethics, namely, safety, security, discrimination, surveillance, deception and harassment, environment, human rights, and bias and fairness impacted by the research conducted in this paper.
Justification:
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: The paper poses no such risks.
Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: The paper does not use existing assets.
13. New Assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: This paper does not involve crowdsourcing nor research with human subjects.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: This paper does involve crowdsourcing nor research with human subjects.
31


Reviews:
[{"review_id": "0Hvajk9iou", "reviewer": "Reviewer_SH1u", "summary": "This paper proposes a novel combinatorial algorithm for the problem known as semi-discrete optimal transport. The proposed method constructs a residual graph by considering the cells of a $\\delta$-expanded Voronoi diagram, a relaxed concept of a weighted Voronoi diagram, as vertices. The algorithm performs augmentation on the residual graph while scaling $\\delta$. This approach significantly reduces the theoretical computational complexity compared to existing methods. Furthermore, the proposed method can be applied to discrete optimal transport with large supports, enabling sublinear time responses to OT queries with respect to support size through preprocessing.", "questions": "* Is this method an implementable and practically applicable algorithm, or is it challenging to apply in real-world scenarios due to large constant factors, thus possessing only theoretical value?", "limitations": "The authors properly address the limitations of the method within the main text.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "* The paper makes a significant contribution to the important problem of semi-discrete optimal transport in the field of machine learning.\n* The proposed method drastically reduces computational complexity compared to existing methods. Additionally, its application to query responses through preprocessing also demonstrates excellent computational efficiency.\n* The paper is exceptionally well-written, making it easy to follow the ideas despite the challenging content.", "weaknesses": "* Understanding that this is a theoretical paper, it is important to note that the lack of numerical experiments makes it difficult to assess the practical applicability of the proposed method.\n* While this is not a weakness of the paper, I am not able to fully verify the validity or the details of the proofs and methodology, and therefore cannot guarantee their correctness.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your insightful review. We answer your question below.\n\n>Is this method an implementable and practically applicable algorithm, or is it challenging to apply in real-world scenarios due to large constant factors, thus possessing only theoretical value?\n\n**Response:**\nThe main contribution of the paper is a new algorithmic framework for the semi-discrete OT problem, which we believe has much potential. Overall, our algorithm is practical and the constants hiding in the big-O notation are small. However, faithful implementation of our algorithm requires (a) constructing and dynamically maintaining arrangements of Voronoi diagrams, and (b) computing the exact mass inside any given triangle. Developing efficient and robust software that combines our algorithm with black-boxes (a) and (b) is a significant task. An interesting direction of future research is to explore which existing geometric software libraries, including the ones based on GPUs, can be used for our setting.\n\nRecollect that our algorithm does not make any assumptions on the smoothness of the continuous distribution. \nIn future work, we will investigate whether we can significantly simplify our algorithm if we are willing to make certain smoothness assumptions (similar to the ones made in existing work [1, 2]) about the continuous distribution.\n\n----\n[1] Oliker, Vladimir I and Laird D Prussner. ``On the numerical solution of the equation $\\frac{\\partial^{2}z}{\\partial x^2} \\frac{\\partial^2 z}{\\partial y^2} - \\left (\\frac{\\partial^2 z}{\\partial x \\partial y} \\right) = f$ and its discretizations, $i$.'' Numerische Mathematik, 54(3):271–293, 1989.\n\n[2] Merigot, Quentin, and Boris Thibert. \"Optimal transport: discretization and algorithms.\" In Handbook of numerical analysis, vol. 22, pp. 133-212. Elsevier, 2021.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "5CNsBwzTkj", "reviewer": "Reviewer_ZH4D", "summary": "The paper studies the semi-discrete optimal transport problem, i.e. a generalization of bipartite matching where one side of the graph is not finite, but instead a probability distribution on some subset of the R^d. \nIt gives an algorithm that computes an optimal transport plan (up to some additive $\\varepsilon) to such instances subject to an oracle that can integrate the density function of the continuous distribution in constant time over some given triangle. \nThe algorithm proposed is an adaption of the solution proposed by Agarwal et al. [SODA 24] which makes some improvements to depress the degree of the polynomial dependence of the runtime on the size of the discrete distribution in the input.\nThe authors further note that their methods yield a kind of sub-linear time online algorithm for discrete optimal transport where one (large) side of the bipartite graph is fixed and the smaller side is subject to updates. This is achieved by representing the larger side as a continuous distribution which can be efficiently sampled.\n\nNote that I was not able to check the appendix where essentially all proofs live.", "questions": "No questions", "limitations": "Limitiations and impact have been adressed fully, this is a theoretical result with no negative consequences to be expected.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is clear and easy to follow, and improves the state of the art on a very relevant topic (optimum transport).", "weaknesses": "- The work is largely a reworking of other recent results in the field, in particular the paper of Agarwal et al. at SODA24. It's not clear that improving the polynomial dependence here is highly relevant.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thoughtful review. We address your concern below.\n\n> The work is largely a reworking of other recent results in the field, in particular the paper of Agarwal et al. at SODA24. It's not clear that improving the polynomial dependence here is highly relevant.\n\n**Response:** Improvement in the run time is only one part of our contribution. What we consider the most important contribution is the new algorithmic framework we propose for the semi-discrete OT problem.\nThe algorithm from SODA 2024 discretizes a continuous distribution and then solves an instance of the discrete OT problem defined on these points in $\\tilde{O}(n^9)$ time. \nIn contrast, our algorithm extends the classical combinatorial primal-dual approach for discrete OT to the semi-discrete setting. It applies ideas such as augmenting paths on residual graphs directly to continuous regions (rather than samples from the regions). To our knowledge, our paper is the first to introduce this framework. As a result of the new combinatorial framework, we obtain a significant reduction in execution time to $\\tilde{O}(n^4)$.\n\nWe remark that polynomial improvements to discrete OT have had a significant impact on ML applications [1, 2, 3, 4]. \nSimilar improvements for the semi-discrete OT problem is an important challenge and we make significant progress in addressing this challenge.\n\n\n----\n[1] Cuturi, Marco. \"Sinkhorn distances: Lightspeed computation of optimal transport.\" NeurIPS 2013.\n\n[2] Altschuler, Jason, Jonathan Niles-Weed, and Philippe Rigollet. \"Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.\" NeurIPS 2017.\n\n[3] Lahn, Nathaniel, Deepika Mulchandani, and Sharath Raghvendra. \"A graph theoretic additive approximation of optimal transport.\" NeurIPS 2019.\n\n[4] Jambulapati, Arun, Aaron Sidford, and Kevin Tian. \"A direct $\\tilde{O}(1/\\varepsilon)$ iteration parallel algorithm for optimal transport.\" NeurIPS 2019.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ncSkbYbpuh", "reviewer": "Reviewer_Eii6", "summary": "The paper presents a primal-dual algorithm for approximately solving the semi-discrete optimal transport problem. The algorithm runs $\\log \\Delta/\\epsilon $ scales, starting with the scale $\\delta = \\Delta^2$ and halving it in each round. The idea is to maintain a $\\delta$-feasible weight function during the course of the algorithm, which in each round constructs a residual graph from the three levels of expansion of the Voronoi diagram, and augments the transport plan along augmenting paths. To ensure the invariants and that the algorithm terminates, the algorithm considers only admissible paths and performs a procedure to eliminate cycles after each round. The proposed algorithm improves the running time from $n^9$ in Agarwal et al. to $n^4$ and the size of the graph from $n^5$ to $n^3$. The algorithm also extends to any dimension $d > 2$ and any $p\\ge 1$-Wasserstein distance.", "questions": "- Could the authors provide some more background on the construction of Voronoi diagram, ie, the complexity of the construction?\n- Is there any lower bound on the runtime of algorithms for the semi-discrete OT problem?\n- How is Theorem 1.2 compared with other algorithms for the discrete OT problem?", "limitations": "The authors adequately addressed the limitations of the paper.", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper presents a set of strong results that significantly improve known results on semi-discrete optimal transport problem. \n- Although I didn't verify all the proofs, I think the algorithm and procedures and the idea of using admissibility for DFS to avoid and remove cycles are all reasonable.", "weaknesses": "- There isn't any big weakness I can see.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thorough review. We answer your questions below.\n\n----\n\n>Could the authors provide some more background on the construction of Voronoi diagram, ie, the complexity of the construction?\n\n**Response:**\nFor 2 dimensions, the weighted Voronoi diagram under the squared Euclidean distance (also known as the Laguerre diagram or the power diagram) can be constructed in $O(n\\log n)$ time [1, 2]. For higher dimensions $d>2$, the construction time would be $O(n^{\\lceil (d+1)/2\\rceil})$ [3].\n\n----\n\n[1] Fortune, Steven. \"Voronoi diagrams and Delaunay triangulations.\" In Handbook of Discrete and Computational Geometry, pp. 705-721, 2017.\n\n[2] Sugihara, Kokichi. \"Laguerre Voronoi diagram on the sphere.\" Journal for Geometry and Graphics 6, no. 1 (2002): 69-81.\n\n[3] Aurenhammer, Franz. \"Power diagrams: properties, algorithms and applications.\" SIAM Journal on Computing 16, no. 1 (1987): 78-96.\n\n----\n----\n\n>Is there any lower bound on the runtime of algorithms for the semi-discrete OT problem?\n\n**Response:** In 2 dimensions, there are no known sub-quadratic algorithms even for the discrete version of the OT problem. The semi-discrete OT problem is a generalization of the discrete OT problem. So, to obtain a runtime better than $O(n^2)$, one would expect that a subquadratic-time algorithm for the discrete OT would come first. \n\nFor higher dimensions, as stated in our paper, computing an $\\varepsilon$-close semi-discrete transport plan in a time that is polynomial in $d$ and $\\log 1/\\varepsilon$ is known to be \\#$P$-Hard [4].  \n\n----\n[4] Taşkesen, Bahar, Soroosh Shafieezadeh-Abadeh, and Daniel Kuhn. \"Semi-discrete optimal transport: Hardness, regularization and numerical solution.\" Mathematical Programming 199, no. 1 (2023): 1033-1106.\n\n----\n----\n\n> How is Theorem 1.2 compared with other algorithms for the discrete OT problem?\n\n**Response:**\nThank you for this question. Given $\\mu$ (with support size $N$) and $\\nu$ (with support size $k$), the execution times of the best-known discrete OT algorithm is almost-linear in the number of edges of the bipartite graph, i.e., $(Nk)^{1+o(1)}$ [5, 6].  In contrast, by preprocessing $\\mu$ in near-linear time, we can compute a discrete OT plan for any query $\\nu$ in $O(\\sqrt{N}k^{4}\\log 1/\\delta)$ time. Thus, when $k$ is sufficiently small (say $k < N^{1/6}$), our algorithm is faster than all existing discrete OT algorithms. We will highlight this improvement to discrete OT algorithms in the next version of our paper.  \n\n----\n\n[5] Chen, Li, Rasmus Kyng, Yang P. Liu, Richard Peng, Maximilian Probst Gutenberg, and Sushant Sachdeva. \"Maximum flow and minimum-cost flow in almost-linear time.\" FOCS 2022.\n\n[6] Agarwal, Pankaj K., Kyle Fox, Debmalya Panigrahi, Kasturi R. Varadarajan, and Allen Xiao. \"Faster Algorithms for the Geometric Transportation Problem.\" SoCG 2017.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper presents a combinatorial algorithm that solves a well-studied task, known as semi-discrete optimal transport, where the goal is to compute the cheapest way to transport the mass from a continuous distribution to a discrete one (where the transport cost scales with the distance). This work, while building on recent algorithmic results for this problem, develops novel conceptual ideas to obtain a significantly faster algorithm. The reviewers agreed that this is an interesting contribution that is appropriate for presentation at NeurIPS.
Decision: Accept (poster)