A Combinatorial Algorithm for the Semi-Discrete Optimal Transport Problem

Authors: pankaj k agarwal
sharath  raghvendra
pouyan  shirzadian
keegan  yao
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 00:20:13

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$.

Main paper: 
 introduction: 
 Optimal Transport (OT) is a 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 to transport mass from (a possibly continuous) distribution µ that is stored compactly (say, using a deep neural network) to a discrete distribution ν. In recent years, the semidiscrete OT has been used in 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 [8, 9, 19] .More formally, a semi-discrete transport plan τ between a continuous probability distribution µ defined over a compact support A ⊂ R d and a discrete distribution ν with a support set B of n points in R d is a distribution over A×B whose marginals are dominated by µ and ν, i.e., τ : A×B → R ≥0 is a transport plan between µ and ν where b∈B τ (a, b) ≤ µ(a) for all a ⊆ A and A τ (a, b) da ≤ ν(b) for all b ∈ B. A transport plan τ is complete if τ transports all mass of µ to ν. For any fixed p ≥ 1, the p-Wasserstein cost of τ is given by ¢(τ ) := A b∈B ∥a -b∥ p τ (a, b) da 1/p . In the semi-discrete p-Wasserstein problem, the goal is to compute a complete transport plan between µ and ν of the minimum cost. An optimal plan is referred to as the W p -OT plan and its cost is referred to as the W p -distance. For any parameter ε > 0, a transport plan τ is called an ε-close W p -OT plan if ¢(τ ) ≤ ¢(τ * ) + ε, where τ * is a W p -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 b ∈ B [12] . The choice of weights guarantees that mass at b is equal to the mass of µ inside the Voronoi cell (also called the Laguerre cell) of b and the optimal transport plan simply transports the mass at b 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 log n) time [24, 57] . For higher dimensions d > 2, the construction time would be O(n ⌈(d+1)/2⌉ ) [11] .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, 55] , they improve the mapping between continuous latent spaces and discrete data in Variational Autoencoders (VAEs) [9] , 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 ε-close transport plan is known to be #P -Hard with respect to d and log 1/ε, i.e., an algorithm with an execution time that is polynomial in d and log 1/ε 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 a continuous distribution µ, and computing an optimal discrete transport plan between the empirical distribution µ n defined on the n samples and ν [43] . It has been shown that this empirical p-Wasserstein distance converges to the true semi-discrete p-Wasserstein distance at a rate of n -1/2p [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, 19] .For the discrete OT problem, there are several near-optimal scalable exact and approximation algorithms [1, 4, 5, 20, 25, 26, 36, 39, 52, 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 semidiscrete transport plans remain elusive. There are algorithms to compute an ε-close semi-discrete transport plan using numerical solvers [12, 15, 17, 21, 37, 38, 41, 48] , entropic regularization [6, 16, 32] , and multiscale approaches [40, 44] . The execution time of all these algorithms is exponential in both d and log 1/ε. Furthermore, their convergence relies on a smoothness parameter of µ. For instance, a notable numerical algorithm by Oliker and Prussner [48] assumes that for a point b ∈ B, a small change in y(b) 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/ε) iterations, where each iteration requires the computation of a weighted Voronoi diagram which takes n O(d) time. Their algorithm slows down when µ 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.Recently, Agarwal et al. [5] described a cost-scaling paradigm to compute an ε-close semi-discrete transport plan. Their algorithm executes log ∆ ε scales, where ∆ is the diameter of A ∪ B. Within each scale, they create an instance of the discrete OT problem of size O(n 5 ) in R 2 (and n O(d) in R d ). 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 ε-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 O(n 9 log ∆ ε ) in R 2 and n O(d) log ∆ ε 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 a continuous probability distribution defined on a compact set A ⊂ R 2 , ν a discrete probability distribution with a support B ⊂ R 2 of size n, and ε > 0 a parameter. Suppose there exists an oracle that, given a triangle φ, returns the mass of µ inside φ in Φ time. Then, an ε-close W 2 -OT plan between µ and ν can be computed in O(n 3 (Φ + n log n) log ∆ ε ) time in the worst-case, where ∆ is the diameter of A ∪ B.Similar to the algorithm by Agarwal et al. [5] , our algorithm is also based on a cost-scaling approach and executes O(log ∆ ε ) scales. However, the algorithm within each scale is different. Agarwal et al. create a discrete OT instance with n 4 vertices and n 5 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 semi-discrete transport plan in O(n 4 log n) time in the worst case. There are several challenges in extending the combinatorial framework to semi-discrete settings, and overcoming these challenges is one of the main technical contributions of the paper.In more detail, in each scale, our algorithm maintains a (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 O(n 3 ) in the continuous space. Algorithms for discrete OT maintain dual weights for all points in A ∪ 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 shortest path algorithm has been extended to continuous space [47] , but we are not aware of any previous work that extends a combinatorial discrete OT framework to continuous space.Our algorithm extends to any dimension d ≥ 2 and any p ≥ 1 in a straightforward way. For d > 2, the algorithm in Theorem 1.1 can be shown to have an execution time of O(n d+1 (Φ+n log n) log ∆ ε ). Note that the runtime of Oliker-Prussner's algorithm has a factor 1/ε while ours has only log 1/ε. Unlike their algorithm, ours does not make any assumptions on the smoothness of µ. Furthermore, similar to Agarwal et al. [5] , our transport plan approximates an optimal weighted Voronoi diagram within a small additive factor.One consequence of our algorithm is that we can use it to design a data structure that answers ε-close optimal transport queries efficiently. More precisely, consider a large distribution µ with support of size N . We can preprocess this distribution into a data structure that can return the total weight of points inside a query triangle in O(N 1-1/d ) time [59] . By using Theorem 1.1, we can report an ε-additive approximate transport plan to any query distribution ν with support of size k in O(k d+1 N 1-1/d log ∆ ε ) time. Note that this query time is sub-linear in N . There has been some work on computing 1-Wasserstein distance approximately in sub-linear time and for answering nearest neighbor queries under 1-Wasserstein distance [13] , but these algorithms have larger error and do not extend to 2-Wasserstein distance.Theorem 1.2. Let µ be a discrete probability distribution with a support A ⊂ R d of size N . The distribution µ can be preprocessed, in O(N log N ) time, into an O(N ) size data structure so that for a discrete probability distribution ν with a support B ⊂ R d of size k, and a parameter ε > 0, an ε-close W 2 -OT plan between µ and ν can be computed inO(k d+1 (N 1-1/d + k log k) log ∆ ε ) time, where ∆ is the diameter of A ∪ B.If the points in the supports of µ and ν 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 W p -OT plan between µ and ν in O(k d+1 (N 1-1/d + k log k) log ∆ ε ) time for any fixed even value of p. We also note that in 2-dimensions, for k < N 1/4 , we can use the data structure to compute an exact discrete OT plan for any distribution ν in O(k 3 N 1/2 log ∆) time that is faster than any existing OT algorithm that takes (N k) 1+o (1) time [3, 18] .For simplicity in presentation, we restrict our presentation to d = 2 and p = 2.  
combinatorial framework: 
 In this section, we extend the combinatorial framework used by discrete OT algorithms to the semidiscrete settings. Let µ be a continuous probability distribution defined over a compact support A ⊂ R 2 and ν be a discrete distribution with a support set B of n points in R 2 . EQUATIONThe For any δ-feasible transport plan τ, y(•), if τ is a complete transport plan between µ and ν, then τ, y(•) is called a δ-optimal transport plan. Recall that any optimal transport plan transports the mass at b to its weighted Voronoi cell. In a δ-optimal transport plan, however, the mass of each point b ∈ B is transported inside a 2δ-expansion of the Voronoi cell of b. This introduces an additive increase of at most 2δ in the cost of the transport plan.Lemma 2.1. Any δ-optimal transport plan τ, y(•) between µ and ν is 2δ-close.See Appendix A for a proof. Next, given a δ-feasible transport plan τ, y(•), we define a residual graph and an augmenting path. We also introduce the process of augmenting τ along an augmenting path, which allows us to increase the mass transported by τ .Residual graph. Given a δ-feasible (possibly partial) transport plan τ, y(•), we construct a residual graph G δ by first partitioning the support A of µ into regions to form the vertex set of G δ := G(τ, y, δ) and then defining a set of directed edges.Vertex set. For each point b ∈ B, consider the Voronoi cell of b and its δand 2δ-expansions V δ b and V 2δ b . Let X δ denote the arrangement [2] of these 3n cells across all n points of B. See Figure 1 (middle). For each region φ in this arrangement, pick an arbitrary representative point r φ inside φ and assign it a mass of µ rφ := µ(φ), where µ(φ) = φ µ(a) da denotes the mass of µ inside φ. Let A δ denote the set of representative points of all regions in X δ . The vertex set of G δ is the point set A δ ∪ B along with a source vertex s and a sink vertex t. We refer to any point b ∈ B whose mass is not fully transported by τ as a free point and define its excess as the amount of mass of b that is not transported by τ , i.e., ex(b) = ν(b) -a∈A τ (a, b) da. Similarly, any point r φ ∈ A δ is free if τ does not fully transport the mass into the region φ, and its excess is defined asex(r φ ) = µ rφ -b ′ ∈B τ (φ, b ′ ). Edge set. For each pair (r φ , b) ∈ A δ × B, if τ (φ, b) > 0, we add a backward edge directed from r φ to b in the residual graph. Furthermore, if r φ ∈ V 2δb , we add a forward edge directed from b to r φ in G δ . Additionally, we add a forward edge from the source s to every free point b ∈ B and a backward edge directed from every free vertex r φ to t. This completes the description of the residual graph. Lemma 2.2. For any δ > 0 and a δ-feasible transport plan τ , y(•), the residual graph G δ has O(n 2 ) vertices and O(n 3 ) edges.While describing our algorithm, it is useful to have the definition of weighted distance for all the backward edges, including those incident on t. Therefore, we extend the definition of weighted distance to any edge (r φ , t) as follows. Let b φ denote the weighted nearest neighbor of r φ in B, i. Compressing a semi-discrete transport plan. Given a semi-discrete transport plan τ that transports mass from B to A, we construct a transport plan τ from B to A δ as follows. For each pair (r φ , b) ∈ A δ × B, let τ (r φ , b) := τ (φ, b), i.e., we assign the entire mass transported from b to φ to the pair (r φ , b) (Figure 1 (right)). We refer to the transport plan τ as the compressed transport plan. Lemma 2.3. For any δ-feasible semi-discrete transport plan τ, y(•), the compressed transport plan τ along with weights y(•) is also δ-feasible.Conversely, consider a transport plan τ , y(•) from B to A δ . One can compute a semi-discrete transport plan τ that, given a pair (r φ , b) ∈ A δ × B, arbitrarily transports a mass of τ (r φ , b) from b to φ. Lemma 2.4. Any δ-feasible transport plan τ , y(•) from B to A δ can be converted into a δ-feasible semi-discrete transport plan τ from B to A.We say that any compressed transport plan τ is a forest if the edges transporting a positive mass in τ do not create an undirected cycle.Augmentation. Given the residual graph G δ for a δ-feasible transport plan τ , y(•), an alternating path (resp. alternating cycle) is a directed path (resp. directed cycle) in G δ . 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, b 1 , r 1 , . . . , r k , t = b k+1 ⟩ is a directed path from the source vertex s to the sink vertex t in the residual graph. By construction, the vertex b 1 and the vertex r k 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{ex(b 1 ), ex(r k ), min i∈[1,k-1] τ (r i , b i+1 )}. To augment τ along P , set τ (r i , b i ) ← τ (r i , b i ) + bc(P ) for each forward edge (b i , r i ) ∈ P and τ (r i , b i+1 ) ← τ (r i , b i+1 ) -bc(P ) for each backward edge (r i , b i+1 ) ∈ P . Lemma 2.5. The transport plan obtained after augmenting a δ-feasible transport plan τ , y(•) along any augmenting path P in the residual graph is δ-feasible.Consider the following straightforward way to augment a semi-discrete transport plan. Given a δ-feasible semi-discrete transport plan τ , y(•), we can compute an augmenting path P in the residual graph and augment the compressed transport plan τ along P . From Lemmas 2.3 and 2.5, τ remains δfeasible and from Lemma 2.4, the updated δ-feasible transport plan τ can be converted to a δ-feasible semi-discrete transport plan as desired. To obtain a 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 an efficient algorithm, we iteratively compute a set of We extend the definition of admissibility to augmenting paths and alternating cycles as follows. Any augmenting path (resp. alternating cycle) P = ⟨b 1 , r 1 , b 2 , . . . , r k , b k+1 ⟩ is admissible if all triples (b i , r i , b i+1 ), i ∈ [1, k],are admissible. In Figure 2 (right), b 2 transports mass to r 1 , and b 1 is its weighted nearest neighbor. By augmenting along this admissible path, we increase τ (r 1 , b 1 ) and reduce τ (r 1 , b 2 ), thereby transporting more mass to r 1 from its weighted nearest neighbor. 
algorithm: 
 In this section, we present our cost-scaling algorithm that uses the combinatorial framework from Section 2 to compute an ε-close semi-discrete OT plan. Classical discrete OT algorithms assign weights to points in A ∪ B and use them to 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 B creates a new weighted Voronoi diagram, which in turn changes A δ , the discrete representation of A, and thus the vertex set of the residual graph. Since A δ may change significantly in each iteration during the execution of our algorithm, we cannot maintain weights for them. This creates significant challenges as 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 our algorithm to eliminate these cycles (See Section 3.4 for details). First, we present an overview of our algorithm and then present the details.The algorithm runs for O(log ∆ ε ) scales, where ∆ is the diameter of A ∪ B. In each scale δ, our algorithm maintains a transport plan τδ and a weight y(b) for every point b ∈ B. Initially, set δ = ∆ 2 and define y(b) = 0 for all b ∈ B. In each scale δ, execute the following steps.1. Initialization: Set τ δ to be an empty transport plan. Compute the residual graph G δ and the compressed transport plan τδ with respect to τ δ , y(•). 2. Iterations: While τδ is not a complete transport plan:(i) Compute a set of admissible augmenting paths in the residual graph G δ and augment τδ along these paths using the SEARCHANDAUGMENT procedure described in Section 3.2.At the end of this step, there are no admissible augmenting paths in the residual graph.(ii) Adjust the weights of all points of B that are reachable from the source by admissible paths by δ and recompute the set A δ , the residual graph G δ , and the compressed transport plan τδ using the INCREASEWEIGHTS procedure described in Section 3.3.(iii) Update the compressed transport plan τδ and the residual graph using the ACYCLIFY procedure described in Section 3.4, so that the transport plan τδ is a forest and the residual graph does not have any admissible cycles.3. Scale Update: Set δ ← δ/2.After the execution of a scale δ ≤ ε/2, our algorithm terminates by returning a complete semi-discrete transport plan τ δ obtained from the compressed complete transport plan τδ (using Lemma 2.4).Invariants. As shown in Section 4 below, our algorithm iteratively updates the weights y(•) and the transport plan τδ while maintaining the following invariants:(I1) The transport plan τδ , y(•) is δ-feasible, and (I2) At the start of each iteration, the transport plan τδ is a forest and there are no admissible cycles in the residual graphs.Remark. The algorithm by Agarwal et al. [5] creates a discrete instance in each scale of the algorithm by computing the arrangement of the δ-, 2δ-, . . ., and (4n + 1)δ-expanded Voronoi cells of each point b ∈ B. Instead of using such a fine partition to create a discrete instance with O(n 5 ) edges, we work directly with the continuous space, maintain a much smaller residual graph with O(n 3 ) edges, and use our semi-discrete combinatorial framework to find a transport plan.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 SEARCHAN-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, we 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 t. 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 SEARCHANDAUGMENT procedure terminates. We provide the details below.Let ← -G δ be the graph formed by reversing the direction of all the edges of G δ . We conduct our search starting from the sink t 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 = ⟨t = b 0 ⟩ as the search path that the procedure grows. Execute the following steps until the search path Q becomes empty. The algorithm terminates when the search path Q becomes empty, i.e., the procedure marked t as visited and removed it from Q. The following lemma shows the properties of the SEARCHANDAUG-MENT procedure. 1. If Q = ⟨t = b 0 , r 1 , b 1 , . . . , r i , b i ⟩ for some i ≥ 0, (a) If there is an edge from b i to s in ← - G δ (i.e., b i is a free point), then P = ⟨s, b i , r i , . . . ,The change in the weights of B in the INCREASEWEIGHTS procedure requires us to recompute the residual graph and the compressed transport plan. This recomputation may potentially create a cycle in the transport plan or an admissible cycle in G δ . The ACYCLIFY procedure eliminates such cycles and ensures that the invariant (I2) holds at the start of the next iteration. Converting the transport plan into a forest is critical for the efficiency of the SEARCHANDAUGMENT procedure 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 τδ a forest. Then, execute a partial DFS search from each unvisited point b ∈ B similar to the one described in the SEARCHANDAUGMENT 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 in the residual graph. Since canceling admissible cycles could have introduced new cycles in the transport plan, repeat the first step to update τδ and make it a forest. See Appendix B.3 for details. The following lemma shows the properties of the ACYCLIFY procedure. Lemma 3.3. Suppose invariant (I1) holds at the start of the ACYCLIFY procedure. Then, during the execution of the ACYCLIFY procedure, (A1) the transport plan τδ , y(•) remains δ-feasible, and (A2) the transport plan τδ is a forest and there are no admissible cycles in the residual graph.This completes the description of our algorithm. 
algorithm analysis: 
 In this section, we first prove the correctness of our algorithm and then analyze its runtime.Proof of invariants (I1) and (I2). For any scale δ, the initial transport plan τδ is empty. Therefore, τδ along with the weights y(•) is δ-feasible. By properties (S1), (W1), and (A1), the transport plan τδ , y(•) remains δ-feasible in each iteration of our algorithm, and therefore, invariant (I1) holds. The invariant (I2) is a direct consequence of property (A2) in Lemma 3.3.Proof of Correctness. From Invariant (I1), in each scale δ, our algorithm maintains a δ-feasible transport plan τδ , y(•) during its execution. The while loop in Step 2 breaks when τδ is a complete transport plan. Therefore, τδ along with weights y(•) is δ-optimal. From Lemma 2.4, one can convert τδ into a δ-optimal semi-discrete transport plan τ δ . Given that our algorithm terminates when δ ≤ ε/2, from Lemma 2.1, the transport plan returned by our algorithm is ε-close, as desired.Efficiency of the algorithm. The SEARCHANDAUGMENT procedure runs a partial DFS on the residual graph with O(n 3 ) edges. Upon finding an augmenting path P , the procedure augments the transport plan along P in O(|P |) time. In Lemma C.1, we use invariant (I2) to show that the total length of all augmenting paths found by the procedure is O(n 3 ). Hence, the running time of the SEARCHANDAUGMENT procedure is O(n 3 ). The INCREASEWEIGHTS procedure stores a sorted list of neighbors for each region and executes a DFS procedure in the residual graph, which takes O(n 3 log n) time. After increasing the weight of a subset of points of B, the procedure recomputes the residual graph G δ and the transport plan τδ , which can be done in O(n 2 (Φ + 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 a dynamic tree data structure [56] in O(n 3 log n) time. To eliminate the admissible cycles, the ACYCLIFY procedure relies on a partial DFS that runs in O(n 3 ) time. Therefore, the total time by the ACYCLIFY procedure is O(n 3 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 O(n 2 (Φ + 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 O(n). Lemma 4.1. For each scale δ, the total number of iterations of step 2 of our algorithm is O(n).Proof Sketch. Let τ 2δ , y 2δ (•) denote the 2δ-feasible semi-discrete transport plan computed by our algorithm for scale 2δ, and let τ δ , y δ (•) denote a partial semi-discrete transport plan maintained during the execution of step 2 of our algorithm. Let X 2δ (resp. X δ ) denote the partitioning of the set A with respect to weights y 2δ (•) (resp. y δ (•)). Let Y be the arrangement of all 3n cells used to construct X 2δ with all 3n cells used to construct X δ . Let τ2δ (resp. τδ ) denote the compressed transport plan for τ 2δ (resp.  
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 a large discrete distribution µ and efficiently computes an ε-close OT plan between µ and any query distribution ν in sub-linear time relative to the support size of µ. More precisely, given a discrete distribution µ with a (possibly large) support A of N points in R 2 , we design a data structure that, given a query discrete distribution ν with a support B of k points, computes an ε-close transport plan between µ and ν in O(k 3 (√ N + k log k) log ∆ ε ) time.Additionally, we show that if the support points have bounded integer coordinates and the masses are rational numbers, our data structure can efficiently compute an exact discrete OT plan.At a high level, our data structure interprets the large discrete distribution µ as a continuous distribution and uses a simplex range-searching data structure as an oracle to compute the mass of µ inside a query triangle. In this way, for any query distribution ν, one can execute the steps of our combinatorial semi-discrete algorithm to compute an ε-close transport plan between µ and ν. 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( √ N ) time [59] . Given a query discrete distribution ν, one can use our algorithm from Section 3 in conjunction with the RS-DS to compute an ε-close transport plan between µ and ν in O(k 3 (√ N + k log k) log ∆ ε) time leading to Theorem 1.2. Consider the special case where the points in A ∪ B have positive integer coordinates bounded by λ, the mass of µ (resp. ν) on each point a ∈ A (resp. b ∈ B) is a rational number of form xa T (resp. x b T ) for positive integers T and x a (resp. x b ), and p is an even number. In this case, the p-Wasserstein cost of any transport plan between µ and ν is an integer multiple of 1 T , and therefore, any 1 2T -close transport plan between µ and ν 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 a query time of O(k 3 ( √ N + k log k) log(λT )), leading to the following corollary.Corollary 5.1. When the points in the supports of the distributions µ and ν have integer coordinates bounded by λ and the mass on each point is a rational number of form x T , our data structure computes, for any even number p, an optimal solution for the p-Wasserstein problem between µ andν in O(k 3 ( √ N + k log k) log(λT )) time. 
conclusion: 
 In this paper, we designed a novel combinatorial framework for the semi-discrete optimal transport problem and used it to compute an ε-close semi-discrete transport plan. We also used this framework to design a data structure that stores a discrete distribution µ over a large support of size N and can compute ε-close OT cost between µ and a query discrete distribution ν in a time that is sub-linear in N . We conclude with the following question: Can we use our combinatorial framework to compute an ε-close semi-discrete transport plan between high dimensional distributions in poly(n, d, 1/ε) time. 
 
 Appendix: 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 τ, y(•) between µ and ν is 2δ-close. (EQUATIONLet τ * denote an optimal transport plan from µ to ν. Then,EQUATIONCombining Equations ( 3) and ( 4) and plugging τ and τ * in Equation ( 2 Therefore, ¢(τ ) ≤ ¢(τ * ) + 2δ and τ is 2δ-close.Residual graph. Next, we show that for any δ > 0 and any δ-feasible transport plan τ , y(•), the residual graph G δ has O(n 2 ) vertices and O(n 3 ) edges, leading to Lemma 2.2. To do so, we show below that the partitioning X δ consists of O(n 2 ) regions. We then conclude that the number of vertices of G δ is O(n 2 ). Furthermore, since the residual graph is a bipartite graph between set B of size n and set A δ of size O(n 2 ), the number of edges of G δ would be at most O(n 3 ).Recall that the partitioning X δ is constructed as the arrangement of all Voronoi cells, δ-expansions, and 2δ-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.e., regions) in X δ is O(|V|). Therefore, to show that X δ has O(n 2 ) regions, we show that the number of vertices of this arrangement is O(n 2 ).EQUATIONb ) denote the Voronoi cell (resp. δ-expanded Voronoi cell, 2δ-expanded Voronoi cell) of b, and let y 0 b (resp. y δ b , y 2δ b ) denote the weights of B used to compute the cell. Note that each Voronoi cell V i b for each b ∈ B and i ∈ {0, δ, 2δ} has at most n vertices. Hence, the total number of Voronoi vertices in the arrangement X δ is O(n 2 ). Next, we count the number of intersection points of these Voronoi cells. For each pair of points b 1 , b 2 ∈ B, consider the pair of cells V i b1 and V j b2 , for some i, j ∈ {0, δ, 2δ}. We show that V i b1 and V j b2 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 b and b ′ , i.e., all points x ∈ R 2 such that d y (x, b) = d y (x, b ′ ). Note that under the squared Euclidean distances, the weighted bisector of two points is a straight line. Let v denote an intersection point of V i b1 and V j b2 , and suppose v lies on the segment of V i b1 corresponding to the weighted bisector of b 1 and a point b 3 ∈ B (See Figure 4 ). Note that if b 3 = b 2 , then V i b1 and V j b2 share a segment containing v, which means thatd y (v, b 1 ) -i = d y i b 1 (v, b 1 ) = d y i b 1 (v, b 2 ) = d y (v, b 2 ),and alsod y (v, b 2 ) -j = d y j b 2 (v, b 2 ) = d y j b 2 (v, b 1 ) = d y (v, b 1 ).Therefore, in this case, i = j = 0, and the two endpoints of the shared segment is counted toward the number of Voronoi vertices of the arrangement. Hence, we assume b 3 ̸ = b 2 . In this case, EQUATIONUsing Equation ( 6), the point v also lies in the segment of V j b2 that corresponds to the weighted bisector of b 2 and b 3 with respect to weights y j b2 . Thus,d y i b 1 (v, b 1 ) = d y i b 1 (v, b 3 ) = d y j b 2 (v, b 3 ) = d y j b 2 (v, b 2 ).In other words, the point v satisfiesd(v, b 1 ) -(y(b 1 ) + i) = d(v, b 2 ) -(y(b 2 ) + j),i.e., the point v lies on the weighted bisector of b 1 and b 2 (the blue dashed line in Figure 4 ). Since the weighted bisector is a straight line, it intersects the convex polygon V i b1 in at most 2 points; hence, the two Voronoi cells V i b1 and V j b2 intersect each other in at most 2 points. Since the are O(n 2 ) 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 2 ), as claimed. Let v be an intersection point of the d Voronoi cells V i1 b1 , V i2 b2 , . . . , andV i d b d . Suppose v lies on the (d -1)-dimensional hyperplane of V i1b1 that is shared between b 1 and a point b d+1 ∈ B. In this case,EQUATION)We claim that for any k ∈ {1, . . . , d}, y i k b k (b d+1 ) = y(b d+1 ). Consider the following two cases:• If b d+1 = b j for some j ∈ {2, . . . , d}, thend y (v, b 1 ) -i 1 = d y i 1 b 1 (v, b 1 ) = d y i 1 b 1 (v, b j ) = d y (v, b j ),and alsod y (v, b j ) -i j = d y i j b j (v, b j ) = d y i j b j (v, b 1 ) = d y (v, b 1 ).Therefore, i 1 = i j = 0 and for any k ∈ {1, . . . , d}, y Therefore, from Equation (7) ,EQUATIONand for any k ∈ {1, . . . , d}, since v lies on a (d -1)-dimensional hyperplane of V i k b k , it has to lie on the hyperplane ofEQUATION)From Equation ( 9 Compressing the transport plan. Recall that for any semi-discrete transport plan τ from B to A, a set of weights y(•), and a parameter δ > 0, the compressed transport plan τ is a discrete transport plan from B to A δ . Note that the definition of δ-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 δ-feasible if τ transports the mass of each point b ∈ B to points of A δ within its 2δ-expanded Voronoi cell. Any complete transport plan τ from B to A δ that is δ-feasible along with weights y(•) is called δ-optimal.Lemma A.2. Any δ-optimal compressed transport plan τ , y(•) from B to A δ is 2δ-close. ′ from B to A δ define ¢ y (τ ′ ) := (r,b)∈A δ ×B d y (r, b)τ ′ (r, b). Then, ¢ y (τ ′ ) = (r,b)∈A δ ×B d y (r, b)τ ′ (r, b) = (r,b)∈A δ ×B d(r, b) -y(b) τ ′ (r, b) = ¢(τ ′ ) -EQUATION)Let τ * denote an optimal transport plan from B to A δ . Then,EQUATION)Combining Equations ( 11) and ( 12) and plugging τ and τ * in Equation (10) ,¢(τ ) = ¢ y (τ ) + b∈B y(b)ν(b) ≤ (r,b)∈A δ ×B d y (r, b r )τ (r, b) + 2δ + b∈B y(b)ν(b) ≤ ¢ y (τ * ) + 2δ + b∈B y(b)ν(b) = ¢(τ * ) + 2δ.Therefore, ¢(τ ) ≤ ¢(τ * ) + 2δ and τ is 2δ-close.The following observation, which is straightforward from the construction of the partitioning X δ , helps in proving Lemmas 2. Since (S3) holds for G 0 , . . . , G i-1 , for any visited point b ∈ V i-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 b ′′ that is reachable from b by an admissible path in ←--G i-1 (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.e., any point b ′′ ∈ B that has an admissible path to the visited point b in the residual graph G i-1 is also visited. By this observation, there are no admissible paths from any unvisited point in B to any visited point in G i-1 . From Lemma B.1, for any newly formed admissible triple (b, r, b ′ ), we have b ∈ V i-1 and b ′ ∈ U i-1 . Thus, all admissible triples formed after augmenting τ i-1 δ along P i 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 i as well.Property (S2). We use the property (S3) to show that (S2) holds. For any point b ∈ B that is marked as visited, as discussed above, if (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 b f ∈ B is unvisited, b f does not have an admissible path to any visited point b ∈ 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 SEARCHANDAUGMENT 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 ′ ) denote an admissible triple in G ′ that is not admissible in G. Assuming that (S3) holds prior to augmentation along P , the point b is marked as visited, and b ′ is marked as unvisited. After the execution of the SEARCHANDAUGMENT procedure, no admissible augmenting paths remain in the residual graph, i.e., there are no admissible paths from the source vertex s to the sink vertex t. The INCREASEWEIGHTS procedure increases the weights of the subset of points in B that are reachable from s by admissible paths to expand their Voronoi cells and to create new admissible triples in the residual graph. For instance, in Figure 6 (left), the path from s to t is not admissible (the triple (b 2 , r 2 , b 3 ) is not admissible, as b 3 has a lower weighted distance to r 2 than b 2 ). As shown in Figure 6 (right), the INCREASEWEIGHTS procedure then increases the weight of the points b 1 and b 2 (which are reachable from s), leading to the formation of an admissible augmenting path (note that upon updating the weights, the regions corresponding to r 1 and r 2 have slightly changed). The details of the INCREASEWEIGHTS procedure are described below. After the DFS procedure terminates, for each point b ∈ K, set y(b) ← y(b) + δ. This completes the description of the DFS step. We next describe how to recompute the residual graph and the compressed transport plan with respect to the updated weights.Let y(•) (resp. y ′ (•)) denote the weights of the points in B after (resp. before) the weight updates, let X δ (resp. X ′ δ ) denote the partitioning of the set A with respect to weights y(•) (resp. y ′ (•)), and let A δ (resp. A ′ δ ) denote the set of representative points of the regions in X δ (resp. X ′ δ ). Furthermore, let τ ′ δ denote the transport plan maintained by the algorithm for partitioning X ′ δ . To compute the new transport plan τδ for the point set A δ , the INCREASEWEIGHTS procedure first computes the arrangement Y of all 3n cells used to construct X ′ δ with all 3n cells used to construct X δ , i.e., Y is the arrangement of Voronoi cell, δ-expanded Voronoi cell, and 2δ-expanded Voronoi cell of each point b ∈ B both before and after weight updates. See Figure 7 . For each region φ ∈ X ′ δ ∪ X δ , let C(φ) ⊆ Y denote the set of regions of Y that lie inside φ. For each region ϱ ∈ Y, pick an arbitrary representative point r ϱ inside ϱ. We denote the set of all representative points of the regions in Y by Y .The INCREASEWEIGHTS procedure first converts τ ′ δ to a transport plan τ over the finer partitioning Y × B by simply splitting each region φ ′ in X ′ δ (and its mass transportation) to fine regions of Y inside φ ′ . The procedure then uses τ to construct a transport plan τδ over A δ × B by merging the regions of Y (and their mass transportation) to regions of X δ . The details are provided next. The goal of the ACYCLIFY procedure is to ensure that the invariant (I2) holds, i.e., that at the beginning of each iteration of our algorithm, the transport plan is a forest and there are no admissible cycles in the residual graph. The procedure runs in three steps: (1) make τδ a forest, as described in Section B.3.1, (2) cancel any admissible cycles from the residual graph, as described in Section B.3.2, 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 a forest in the third step is essential for invariant (I2), as canceling admissible cycles might introduce cycles in the transport plan. To remove all admissible cycles from the residual graph, we use a partial DFS similar to the one described in the SEARCHANDAUGMENT procedure. Let ← -G δ be the graph formed by reversing the direction of all the edges of G δ . The procedure first marks all points of B and all backward edges as unvisited and defines U := B as the set of unvisited points. While there exists an unvisited point b ∈ B, the procedure initializes a partial DFS by setting Q = ⟨b = b 1 ⟩ and searches as follows until Q becomes empty. 1. If Q = ⟨b 1 , r 1 , . . . , b i ⟩ for some i ≥ 1,(r t , b t )) in C, set τδ (r t , b t+1 ) ← τδ (r t , b t+1 ) + bc(C) (resp. τδ (r t , b t ) ← τδ (r t , b t ) -bc(C)).In this way, at least one of the backward edges of the cycle C is removed from the residual graph and the cycle has vanished. Lemma 3.3. Suppose invariant (I1) holds at the start of the ACYCLIFY procedure. Then, during the execution of the ACYCLIFY procedure, (A1) the transport plan τδ , y(•) remains δ-feasible, and (A2) the transport plan τδ is a forest and there are no admissible cycles in 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, b) ∈ A δ × B, the point b is a 2δ-weighted nearest neighbor of r, i.e., the transport plan obtained after the first step of the ACYCLIFY procedure is δ-feasible.Next, we show property (A1) and (A2) for the second step of the procedure. By construction, all triples on the search path maintained by the procedure are admissible, and therefore, any cycleThe SEARCHANDAUGMENT 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 ∈ B (resp. r ∈ A δ ), marks the point b (resp. the backward edge (b ′ , r) used to reach r) as visited and does not visit the point b (resp. the backward edge (b ′ , r)) again in the same execution. Upon finding an augmenting path P , the procedure augments the transport plan along P in O(|P |) time. Let ⟨P 1 , . . . , P k ⟩ denote the set of all augmenting paths found by the SEARCHANDAUGMENT procedure. In 2 dimensions (resp. d dimensions), since the residual graph has O(n 3 ) edges (resp. O(n d+1 )), the running time of the procedure would be O Proof. Let τ 0 δ denote the transport plan maintained by the algorithm at the beginning of execution of the SEARCHANDAUGMENT 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 v consisting of augmenting paths whose bottleneck capacity is determined based on the residual capacity of its endpoints, and (2) set P e consisting of augmenting paths whose bottleneck capacity is determined based on mass transportation over its backward edges. We first show that |P v | = O(n 2 ) and then show the same bound for P e . Since each augmenting path has a length of at most 2n, we then conclude that the total length of all augmenting paths is O(n 3 ). Proof. Let τ 2δ , y 2δ (•) denote the 2δ-feasible transport plan computed by our algorithm for scale 2δ, and let τ δ , y δ (•) denote a partial transport plan maintained during the execution of step 2 of our algorithm. Let X 2δ (resp. X δ ) denote the partitioning of the set A with respect to weights y 2δ (•) (resp. y δ (•)). Let Y be the arrangement of all 3n cells used to construct X 2δ with all 3n cells used to construct X δ , i.e., Y is the arrangement of Voronoi cell, δ-expanded Voronoi cell, and 2δ-expanded Voronoi cell of each point b ∈ B with respect to weights y δ (•) along with the Voronoi cell, 2δ-expanded Voronoi cell, and 4δ-expanded Voronoi cell of each point b ∈ B with respect to weights y 2δ (•). For each region ϱ ∈ Y, pick an arbitrary representative point r ϱ inside ϱ. We denote the set of all representative points of the regions in Y by Y . Let τ2δ (resp. τδ ) denote the compressed transport plan for τ 2δ (resp. τ δ ) using the partitioning Y. Note that the partitioning Y is a refinement of both partitionings X δ and X 2δ . Define τ ′ := τ2δ -τδ . Equivalently, y δ (b k ) -y 2δ (b k ) ≤ 6nδ. By property (W2) in Lemma 3.3, our algorithm increases the weight of the point b k by δ in each iteration while its mass is not fully transported by τδ ; therefore, the point b k cannot remain free after 6n iterations, i.e., after O(n) iterations, there cannot be any remaining free points in B and the step 2 of our algorithm terminates after O(n) iterations.(n 3 + k i=1 |P i |) (resp. O(n d+1 + k i=1 |P i |)). In Lemma C.1, we show that k i=1 |P i | = O(n 3 ) (resp.Discrete OT algorithms use the weights assigned to A ∪ B to define admissible edges. Since A δ evolves during the execution of our algorithm, we cannot maintain weights for them, forcing us to define admissibility on a sequence of two consecutive edges (b1, rφ) and (rφ, b2) rather than per edge.To perform this step efficiently, in the construction of the residual graph, for each r ∈ A δ , we store the list of all neighbors of r sorted in increasing order of their weighted distance.Here, the constant is O(n d ) hides a factor of 3 d . 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: * 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.
* 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.
questions: 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: limitations: The authors properly address the limitations of the method within the main text.
paper_summary: paper_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. 
It 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. 
The 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.
The 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.

Note that I was not able to check the appendix where essentially all proofs live.
strength_weakness: strength_weakness: - 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.
questions: questions: No questions
limitations: limitations: Limitiations and impact have been adressed fully, this is a theoretical result with no negative consequences to be expected.
paper_summary: paper_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.
strength_weakness: strength_weakness: - There isn't any big weakness I can see.
questions: questions: - Could the authors provide some more background on the construction of Voronoi diagram, ie, the complexity of the construction?
- Is there any lower bound on the runtime of algorithms for the semi-discrete OT problem?
- How is Theorem 1.2 compared with other algorithms for the discrete OT problem?
limitations: limitations: The authors adequately addressed the limitations of the paper.

Comments: title: Comments on the rebuttal
comment: I thank the authors for the response. I maintain my current score.
I have reviewed the authors' response. My concerns have been resolved. I will maintain my current score.