Efficient Algorithms for Robust and Partial Semi-Discrete Optimal Transport
Authors: Pankaj K Agarwal
Sharath Raghvendra
Pouyan Shirzadian
Keegan Yao
Venue: NeurIPS 2025 poster ()
Keywords: Semi-Discrete Optimal Transport
Approximation Algorithms
Robust Optimal Transport
Primary area: optimization
TL;DR: We relate robust and partial variants of optimal transport (OT), introduce novel characterizations of their solutions in the semi-discrete setting, and develop two algorithms for these variants of semi-discrete OT.

Abstract: The sensitivity of optimal transport (OT) to noise has motivated the study of robust variants. In this paper, we study two such formulations of semi-discrete OT in $\mathbb{R}^d$: (i) the $\alpha$-optimal partial transport, which minimizes the cost of transporting a mass of $\alpha$; and (ii) the $\lambda$-robust optimal transport, which regularizes the OT problem using the total variation (TV) distance. First, we provide a novel characterization of the optimal solutions in these settings, showing they can be represented as a restricted Laguerre diagram. Second, we exploit this characterization to establish a strong algorithmic connection between the two problems, showing that any solver for one can be adapted to solve the other with comparable precision. Third, we overcome key challenges posed in extending the cost-scaling paradigm to compute these variants of OT and present an algorithm that computes the exact solution up to $\log (1/\varepsilon)$ bits of precision in $n^{O(d)}\log (1/\varepsilon)$ time, where $n$ is the support size of the discrete distribution.  Finally, we present an $n^{1+o(1)}\varepsilon^{-O(d)}$ time approximation algorithm for the above variants of OT.

Efficient Algorithms for Robust and Partial Semi-Discrete Optimal Transport
Pankaj K. Agarwal
Sharath Raghvendra
Department of Computer Science Duke University
Department of Computer Science North Carolina State University
Pouyan Shirzadian
Keegan Yao
Department of Computer Science Virginia Tech
Department of Computer Science Duke University
Abstract
The sensitivity of optimal transport (OT) noise has motivated the study of robust variants. In this paper, we study two such formulations of semi-discrete OT in RR4: (i) the a-optimal partial transport, which minimizes the cost of transporting a mass of a; and (ii) the A-robust optimal transport, which regularizes the OT problem using the total variation (TV) distance. First, we provide a novel characterization of the optimal solutions in these settings, showing they be represented as restricted Laguerre diagram. Second, we exploit this characterization to establish a strong algorithmic connection between the two problems, showing that any solver for one can be adapted to solve the other with comparable precision. Third, we overcome key challenges posed in extending the cost-scaling paradigm compute these variants of OT and present an algorithm that computes the exact solution up to (og(1/z) bits of precision in n ((d) log(1/E) time, where n the support size the discrete distribution. Finally, we present an n [+(()g -O(()) time approximation algorithm for the above variants of OT.
Introduction
Given two distributions μ and v supported on sets A and B in Rd and a parameter p 2 1, the cost of transporting units of mass from E A to b B is m|la - bll.. The optimal transport cost (also called the p- Wasserstein distance) quantifies the minimum total effort required move mass from to v. In many machine learning applications, one wishes to compute the OT cost between model and target distributions that are continuous. Due to lack of efficient algorithms that work directly on continuous distributions, practitioners often approximate the OT cost by computing it between two discrete empirical distributions formed by samples [5] 17 39]. This OT between two discrete distributions can be computed by using minimum-cost flow algorithm (38
more computationally feasible variant is the semi-discrete optimal transport (SDOT) problem, where one distribution is continuous and the other is discrete. In this case, the OT cost and plan can be characterized by set of weights y: + R=o such that, for any point E B, the additively weighted Voronoi cell Vú (b) of (also called the Laguerre cell) has the same probability mass as b. The OT plan assigns all mass in V, (b) 88. The semi-discrete OT has found numerous applications in
machine learning 6 20 27
Despite its usefulness, OT cost may be unduly influenced outliers in either distribution. Two robust variants help reduce the impact of noise. The first variant, called the a-optimal partial transport (a-OPT), minimizes the cost of transporting mass of E [0, 1] from μ to v, potentially discarding
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

outlier mass [14. The second variant, called the X-robust OT (A-ROT), introduces a regularization via total variation distance. Given a parameter À O. m units of mass can either be transported from E A to b E B incurring a cost of mila - b or be burnt incurring a cost of Xm [4 The X-ROT minimumizes the cost under this tun-orrtranssort trade-off. In this paper, we study both a-OPT and X-ROT in the semi-discrete setting. We provide characterization of them and design both exact and approximation algorithms for them.
Related Work.
Several algorithms have been presented for additively approximating the semi-discrete
37 and entropic regularization 10 24 These algorithms take ( poly(1//5 time compute this approximation.
OT using numerical solvers 8 13
20 29
Recently, Agarwal et al. [ adapted the cost-scaling paradigm that is used for the network-flow algorithm to the semi-discrets settings, obtaining the exact solution with log(1/s) bits of precision in ne ne log 1/c time. Thus, this algorithm achieves polynomial time exact solution with (nly(n) bits of precision, while prior algorithms required exponential (in n)) time to achieve the same accuracy. Agarwal et al. [3] also developed
near-linear-ting Monte Carlo algorithm that provides a relative approximation. It is worth noting that weighted Voronoi diagrams under éy metrics can have complexity as high as n ((d)) in the worst case. This inherent geometric complexity suggests that substantially improving upon the current n
bounds for exact solutions may be difficult.
The o-ootimal partial transport (P-OPT) problem for discrete distributions was introduced by Chapel et al. [4 as a means detect and remove outliers, particularly when the fraction of outliers is known in advance; see also [2 35 for useful properties of a-OPT. Phatak et al. [40] later proposed computing the a-OPT cost as a function of a E [0. 1], allowing automatic outlier detection and removal without prior knowledge of the proportion of outliers. Although the a- a-OPT distance is not a metric, Raghvendra al. 42 introduced an a-OPT based metric, called the RPW distance, that enjoys improved statistical robustness compared to the p-Wasserstein distance. In the semi-discrete settings, it known that the a-OPT routes mass from each discrete point to a subset of its weighted Voronoi cell [. and the support of a partial transport plan satisfies the interior ball condition [2] that is, the transported continuous mass lies within the union of balls of a fixed radius centered at the support of the discrete distribution.
Mukherjee er al. 34 introduced the X-ROT (originally referred to by ROBOT distance). There is substantial work in establishing mathematical properties of X-ROT distance as well as designing algorithms for it [33 36 Most of the existing algorithms for both a-OPT and X-ROT distances are restricted to discrete settings and no prior work has extended these robust formulations to the semi-discrete settings. We address the gap in this paper.
Our Results. We introduce the semi-discrete versions of a-OPT and X-ROT and make four novel contributions related to them.
First, we show that the optimal solution for the semi-discrete a-OPT (a- SDOPT) can be characterized using a restricted Voronoi diagram (see Lemma (.1). The restricted Voronoi diagram combines the weighted Voronoi diagram structure [11] with the interior ball condition 12]. We further establish a stronger characterization than was previously known: in the optimal solution, for each point b B, all continuous mass within its restricted Voronoi cell is transported from b and either (i) all mass of b transported, or (ii) b attains the maximum weight among all points of B. We also establish a similar characterization for the semi-discrete X-ROT (A-SDROT) problem; see Lemma 2.2 in Section 2 To our knowledge, this is the first paper to give such a characterization of partial robust OT in the semi-discrete setting.
Second, we demonstrate that any approxàmate solver for X-SDROT can be used to obtain an approxi- mate solution to the aDSDPPT problem. Until now, these two problems have been studied separately, each with its own family of exact and approximation algorithms. The standard reduction used to establish duality between them [12 does not, however, preserve approximation guarantees, as we highlight in Section 3. We provide refined reduction showing that an approximate solution one problem can still be transformed into an approximate solution to the other with nearly matching accuracy.
Third, using our characterization of partial and robust OT the semi-discrete setting, we design a cost-scaling algorithm for the SDROT problem. Extending cost-scaling approaches to solve optimal partial transport is challenging, even in the fully discrete case [43], and these challenges only increase
2

---

in the semi- setting. We detail this challenge as well as our approach for resolving them in Section[] Roughly speaking. the cost-scaling algorithm maintains an approximate restricted Voronoi diagram, which is successively refined at each scale. The algorithm maintains the invariants that the untransported mass of the continuous distribution lies outside the approximate Voronoi cells. Ensuring this invariant as we move from one scale to the next is challenging and requires new ideas. Our algorithm runs in n (dd) log( A/E) time and computes a solution to the SDROT problem with (og(1/E) bits of precision: here n the support size of the discrete distribution and A is the diameter of the two distributions. Through our reduction, we can also compute a solution for the a-SDOPT problem with a similar performance. We develop a preliminary implementation of this algorithm and experimentally demonstrate that (i) the algorithm executes within the theoretically predicted number of iterations, and (ii) both SDOPT and SDROT predominantly transport inlier mass. See Section 4
Finally, we describe an e-approximation algorithm for the SDROT problem running in near- linerr time n. Our algorithm first reduces SDROT the discrete robust OT problem, then reduces the problem to computing a transport plan on a sparse graph, and finally either uses the recent minimum-com cost flow algorithm [6] or uses the Although the algorithm builds on the one in [1, new ideas are needed at both reduction steps to make it work for partial OT. For p- Wasserstein distance, i.e. the cost is ||a bll", the algorithm takes n 1 4d log A log e time for any constant value of p. For p = 1, we present a faster Monte Carlo algorithm that runs in O(ne -dd---d--1o log A log n log $ = time. Because of the lack of space. this algorithm is presented in Appendix A. See Theorem A.I.
sertipliiaatveewiiilee update method as in 22 28
Characterizing Partial and Robust Semi-Discrete Optimal Transport
In this section, we provide a characterization of partial and robust semi-discrete OT using weighted Voronoi diagrams. We begin by giving a formal definition of the variants of OT.
Partial and Robust OT. Suppose is a continuous probability distribution defined over compact support A C Rd and v is discrete distribution with support set B of points in F For parameter 2 1, define the distance of each pair of points (a, b) E A x B as d(a. b) := ||a bll". A transport plan 5 A x + Ezo is a function whose marginals are dominated by μ and v, i.e. EvEB ((a,)) b) S p(a) for all a C A and JA T(a, b) da S vbb) for all b E B. The cost of transport plan between μ and is defined as ((T) ;= Eoen d(a, b)r b) da. Define M(T) as the total amount of mass transported by transport plan t, i.e., M(T)= = SA EBEB T(a,
For a parameter a E [0. 1]. 7 is an a-partial transport plan if M(T) = a. An a-optimal partial transport (a-OPT) plan is minimum-cost a-partial transport plan.
For a parameter À >> 0, define the X-robust cost of a transport plan T as Wx(T) := <(r)+ (()) A X-robust optimal transport (X-ROT) plan is transport plan with the minimum X-robust cost.
Voronoi cells. For a set of weights y(-) for the points in B, the weighted distance between each point a E A and E B is defined as (o, b) == d(a, b) - y(b). For any point a E A, define its (weighted) nearest neighbor as the point 6 B with the smallest weighted distance to a. For any point b B, the Voronoi cell of b is the set of points of A with b as their weighted nearest neighbor, i.e.,
Vy( == {a E A d b) S dd Vb' B).
See Figure I left). is well known that for any set of weights v(-) for B such that vVV,, = v(b) for all points E B, the transport plan that transports the mass of each point b E B to the mass of μ inside V, is an OT plan. Furthermore, one can prove the existence of such weights using the LP formulation of semi-discrete OT and the strong LP duality; see S 19 23
Restricted Voronoi cells. Given a set of non-negative weights v(-) for points in B, for any point 6 E B, the restricted Voronoi cell of b captures the points of the Voronoi cell of that are within distance of y(b) from 6, i.c.,
RV2(b) == {a E Vy(b) d(a, b) S b(b)).
See Figure rright). We refer to any point E B as balanced (resp. surplus, deficit) if (RVV, (b)V (b)) = (b) (resp. n (b) < (b), μ (RV, (b)) > v(b)).
Valid weights, cover, and cap. A weight function y : B + R2oo is called valid if all points b E B either balanced or surplus and the surplus points have the maximum weight, i.c. for all surplus
3

---

Figure 1: (left) Voronoi cells, and (right) restricted Voronoi cells
points b E B, (b) max y(b).. For a valid weight function y. we define the transport plan induced by y, denoted by Ty : A B - Rzo. as the one that for each point B, transports the mass of b to the continuous mass of inside the restricted Voronoi cell (b). Since each point b E B is surplus or balanced, there is enough mass at transport to (b). We define the cover of v-) to be Cover(y) = Eve B μ RV, y (b)), and its cap to be Cap(y) == maxse v(b).
Characterization of partial and robust OT. The following lemma characterizes a-OPT plan using a set of valid weights with cover of a.
Lemma 2.1 (a-OPT Characterization). Let μ be continuous distribution over compact support A C. Rd a discrete distribution over a set B of n points in R and E [0, 1] parameter. There exists valid weight function y - Rzo with Cover(y) = a such that the transport plan Ty induced by is an a-OPT plan. Furthermore, this property holds for any valid weight function with Cover(y) a.
Proof Sketch. We provide brief sketch of the proof below. The complete proof is provided in Appendix C The a-partial transport problem is linear optimization problem with the linear primal objective function ((T) == =SA EbeB T(a, b)d(a, b)da and mass constraints SA r(a, b)da S v(b), We consider its corresponding dual optimization
EnER (a, b) S μ(a) and Sa LIE T(a b)da problem, which can be written as:
at Σ b(b)v(b) 1 (a))
max
A.
BEB
subject to
((b) (aa) S d(a,b), v(a,))
Φ + 2
(1)
Set (b) = X - @(b). By complementary slackness conditions the optimal solutions for the primal and dual problems, one can argue that if primal constraint JA (a, b)da < v(b) is not tight then y(b) = X, and that if primal variable T(a, b) > 0 then y(b) (aa)) = d(a, b). By feasibility of the optimal dual solution, both of these conditions imply that (i) (b) max y(b') for every surplus point E B, and (ii) T(a, b) > for pairs (a, b) only if a E RV,
Next, we differentiate the dual objective with respect to the unconstrained variable X, argue that the dual objective is maximized when Cover a, and use strong duality and conclude the existence of optimal weights We finally use the feasibility constraints of the dual problem argue that any such weights y satisfying these conditions induce an optimal transport plan.
Next, we characterize a A-robust optimal transport plan using a set of valid weights with cap of X. Lemma 2.2 (A-ROT Characterization). Let be continuous distribution over compact support A C Dd discrete distribution over set B in points in Rd and X > 0 parameter. There exists a valid weight function : + Rş0 with Cap(y) = X such that the transport plan induced by is X-ROT plan. Furthermore, this property holds for any valid weight function y with Cap(y) X.
Suppose y(-) is set of valid weights with a cap of À and cover of a. A transport plan 7 transporting the mass of each point b B to the mass of μ inside (V, would be both c-OPT plan well

---

as A-ROT plan. In the next section, we show strong relationship between the two robust variants of the remi-discrete OT problem.
Relating the two characterizations. For any value À 2 O, let MaxCover(2) denote the largest value [0, such that there exists set of valid weights y() for B with cover of and cap of A. We note that MaxCover(A) gives the maximum value of mass that can be transported if we restrict any mass to be transported at most À distance. Similarly, for any value a 2 let MinCap(a) denote the smallest value À 2 0 such that there exists set of valid weights y(-) for B with Cover(y) = a and Cap(y) X. That is, the maximum distance the mass has be transported in any a- a-OPT plan has to be at least MinCap(a). We show in the following lemma that both functions are monotonically non non-decreasing.
Lemma 2.3. The mappings MaxCover: Rzo - [0, 1] and MinCap: [0, 1] + Rzo are monotoni- cally non-decreasing functions.
From Lemmas 2.2 and 2.3] given an exact black-box solver A for the robust OT problem, one can compute highly accurate a-OPT plan using simple binary search on the value of MaxCover(a). Hence, we obtain the following:
Proposition 2.4. Given procedure for computing semi-discrete robust OT, for any parameter [0, 1] that requires m bits of representation, semi-discrete a- OPT can be computed by making (m) calls to A.
3 Relating Approximate Partial Transport Models
In this section, show that an approximate solver for X-robust OT can be used to compute near-optimal transport plan. We begin by formally defining near-optimal transport plans.
Approximate transport plans. We call an o-partial transport plan (er: -apprxxiiie if c(r) S (1 + ,,))(r")) ) + Ea. i.c., both relative and additive error terms are allowed. Analogously, define transport plan where M(T) = as (Er, aappproxxiie X-robust optimal transport plan if (())+ + (1 a) S (1+ +e,)
+ + a²] }}] Ea. where r* is a A-ROT and a = M(T"). When e, = O, we refer to any (0, caprroxapio a-partial transport plan any (0, S.) -apprrxxmmte X- robust transport plan as Ea close We finally define approximate partial transport plans that also allow for some slack in the amount of mass transported. Define a partial transport plan 7 as (er 6)- approximate a a-partial transport plan if (1) M(T) 2 (1 - d)a, and (2) c(T) I (1 (1+ + Ea. where r" is an optimal o-partial transport plan. The first condition says that the partial transport plan т routes at least (1 5) mass, while the second condition says that the partial transport plan is comparable in cost to the best possible o-partial plan.
Approximation algorithm. We now describe an algorithm ApT for partial OT using a X-robust OT solver. Let E (0, 1) and E, 2 0 be parameters, and let 0 < 8 < as be some constant. Suppose we are given an
algorithm Arr for the X-SDROT problem that besides returning an (E,, aapprrxxiittttt X-robust transport plan (A(A) also returns the amount of mass (A) transports. We perform a binary search on to compute an approximate a-OPT plan.
First initialize values ' 1 minn,, E d(bı, by) and Au = (1- de A diameter of A U and is finite since is assumed that μ has compact support. Run algorithm ArT with inputs År. and ÅR. obtaining approximate partial transport plans Tr. and TR. respectively. If M(TR) < ax we send the untrasported mass in a greedy manner and return the resulting transport plan (As we will see below, if a is close to then ArT might transport less than mass even for A == A). If TL routes more than a mass, then one can greedily route all a mass within a tiny ball of radius most λ around each point of the discrete distribution. Otherwise, we perform a simple binary search on A. Set À = 1 (Ar + AR) and use algorithm AgT to compute an approximate partial transport plan Tx given input A. If M(Tx) 2 a, then assign Ag f- À and Tn + Tx Otherwise, assign AL < À and TL < TX- This maintains the invariant M(TL) S a < M((R) while decreasing the gap between \ and AR. The algorithm stops when dR (1 E, E, + Ea- If M((R) > a, we compute an c- a-partial transport plan T, from Ta by reducing the mass being transported in a greedy manner, i.e. not sending the mass the farthest points of A. See the appendix for details.
where A is the
Pathological example. Before we analyze the performance of this algorithm, we remark that one cannot guarantee Lemmas 2.1 and 2.2 to hold for an approximate partial transport plan. As
5

---

consequence, one cannot immediately use the nice structure of Lemma 2.3 to guarantee that the above binary search returns a near-optimal Q-OPT. Roughly speaking, suppose X MMiCCap(a). Even if we set À = (1 + =)\", an (e, 0)approxximate X-ROT solver may not transport the mass at distance close to X and thus it may return a transport plan T.A.A with M(TAA)) < a. The problem becomes more acute as a approaches 1, as M(TA,N) < even for À > A. For instance, consider the following simple example where A = {a}, B = {b], d(a, = and a = 1. Then a transport plan тλ with (Tx)) = 0 is an (s, -aaapproxiemtee A-robust OT plan for X = 1 1+ Furthermore, for À = 10 a transport plan T with M(Ta) = 1- J is an (e, -aapproxxmate X-robust OT plan (as its cost is 10 5 + (1 - 5) 1 + e). We, however, show that the binary search returns transport plan that has an additive error in the cost or that guarantees to transport at least (1 5)a mass in a near-optimal manner.
The example above highlights the key phenomenon that allows us to provide a guarantee on the quality of the transport plan returned by the above algorithm: as À increases, the maximum amount of mass that is not transported by an approximate X-robust transport plan r decreases. We make this phenomenon precise in the following lemma.
Analysis.
Lemma 3.1. Suppose is an arbitrary (e,, Ea) -approximate X-robust transport plan and ō E 1).
(i) For any a E (0, if M(T) S (1-b0 (1-5)a, then X I (1 8a MMinCap(a) +
&a 80
(ii) Equivalently, for any a E (0, - -). M() S a, then S (1 ;:)) MinCap (1-8) + Co.
The following lemma proves that when the algorithm terminates with interval [}z, AR). the most expensive mass of TR has a cost comparable to Ag-
Lemma 3.2. If AR S (1+E,)AL t then e ((Tn) -((TL) 2 [M(Tn) - M(TL) An (4 45.).
Using Lemma 3.2 we prove that APT is also an approximation algorithm for the Q-OPT problem. Theorem 3.3. Let p. be continuous distribution with compact support A C Rd, and let be discrete distribution with support B C IR and minimum pairwise distance 1. Define A to be the diameter of A AUT B and let (0,1) be a parameter. Suppose Er 2 0 and ART an arbitrary (Er,Ea)-appe approxxmaatii algorithm for the X-robust OT problem. For any ó (0, min { (Ee 1-2 2)). the algorithm APT is an (E,, 11) approoxiiitno algorithm for a-partial transport, where n = O (%t 5r (c.MinCap(a(1 + 8)) + Ea)). For E, = O, ApT is an (0, 55.) aap))opproiniation algoritim. The algorithm ApT makes O (log A) calls to ART-
a
The challenge for the approximate solver arises when some mass needs to be transported much farther than the cost of the transport plan (which is the expected distance of mass being transported), as in cases where has a long tail and a is such that some mass from this tail has to be transported. If is in the range for which the cost of an a-OPT plan Ta is comparable to the maximum distance the mass is transported in Tater, a stronger claim can be made on the quality of the plan returned by Apr.
Corollary 3.4. Let y be a continuous distribution with compuct support C Rd, and let v be discrete distribution with support B C RRd. Suppose En 2 0, a E (0, 1 - Ex) is parameter, and ART is an urbitrary (ar,,a)-approximme --apprxximation algorithm for the X-robust OT problem. Additionally, suppose e(r") 2 MinCaapaa + E.) for some constant E (0, 1). where T" is a-OPT plan. Then, ApT is an (O(E)), (llaa)--aproxiimaton algorithm for a-partial transport.
Finally, if we allow some flexibility in the amount of mass transported by the returned partial transport plan, then one can still obtain similar cost approximation guarantee for arbitrary distributions by running the above binary search but with input being al1 6).
Theorem 3.5. Let μ be a continuous distribution with compact support A C Rd, v a discrete
distribution with support B C Ted, and a E (0, a parameter. Suppose Er, 2 0, E, < a, and 8 5e. Then algorithm Apr is a (E,, 5 8) aaprroximmaiion algorithm for a-partial transport given any arbitrary (at-aaaaappoxiimatio algorithm ART for the X-robust OT problem.
a
4 A Highly Accurate Algorithm for the Robust Optimal Transport Problem
In this section, present algorithm that computes E-close A-robust transport plan, for X, s > 0. between continuous distribution and discrete distribution v in Ö(n log(A/s)) time in R²² and in ((d) log( (A/z) time in Rd, where A is the diameter of A U B.

---

Figure 2: (left) The 6-expanded restricted Voronoi cell of RV, (b), (middle) in a 8-feasible transport plan, the point b E B can transport mass to its 28-expanded restricted Voronoi cell, (right) the decomposition Xz..
4.1 Our Combinatorial Framework
We begin by introducing the notion of a ô-feasible transport plan. Given a weight function 3(-) on points of B and a parameter 8 2 O, the ô-expanded restricted Voronoi cell of is defined
RV (b) == {a E A b(a,b) S v(b) + 8 and (a, S dyla, 8') 5, AV B).
for any pair (a, b) A x B with ((a, b) > 0, we have E RV22 (b), i.e., each point b E B transports mass only to the regions within the 26-expansions of its restricted Voronoi cell (Figure (iggure[mmidlee)) middle)). We refer to any point b E B (resp. a A) as surplus (resp. deficit)
See Figure left). Any transport plan 7 along with weights y(-) 6-feasible
T(a, b) da < v(b) (resp. u(a)) and as balanced if ʃA r(a, b) da vbb (resp. ΣbEB T(a μ(a)).
EBEB ba,
Given weight function v(-) for the points in B, derive weight ū(a) for each point A as y(a) == max(0, max max( d(a, b))}.
BEB
Note that any 0-feasible transport plan transports the of each point b B to regions inside its restricted Voronoi cell. From Lemma 2.2 a X-ROT plan a O-feasible transport plan where Cap(y) A, all surplus points b B have y(b) = X, and all deficit points A have Ø(a) The following lemma shows that any d-feasible transport plan with similar restrictions on surplus and deficit points is 28-close X-robust transport plan.
Lemma 4.1. For any parameter 8 > O, suppose T, T.D(.) denotes a 8-feasüble transport plan, and ler À Cap(y). Suppose the following two conditions hold: (FI) for every deficit point A, T(a) 0. and for all surplus points b E B, y(b) = X. Then, T is a 28-close X-robust transport plan.
Challenges in extending cost-scaling technique. Our algorithm builds on the cost-scaling framework of Agarwal et al. [1. which operates over multiple scales. Their algorithm begins coarse scale 5 = A with all weights ybb = O. In each scale 8, weights are inherited from the previous scale 28. The algorithm discretizes the continuous support using an arrangement of 0-, 5-, and 26-expansions of all Voronoi cells and iteratively computes 8-feasible transport plan by adjusting the weights and increasing the amount of transported mass using augmenting paths. When the plan transports all the mass, the resulting weights are passed on to the next scale 5/2. Their algorithm terminates when = e/2 and returns an s-close transport plan.
Extending their approach to our setting creates new challenges. Unlike in the standard setting, our algorithm computes a d-feasible transport plan that may leave some mass untransported. We must, therefore, cap the weights assigned to each point E B by À (Lemma 4.1 condition (F2)) and restrict the weights p(a) of all deficit regions in A 0 (Lemma 4.1 condition (FT)). A naïve implementation of the scaling paradigm risks generating deficit regions within the restricted Voronoi cells deficit points of A with positive weights violating Lemma 4.1 For example, in Figure 3 the algorithm creates deficit region (green region) inside the restricted Voronoi cells, illustrating such violation. To resolve this, we introduce "consolidating paths", which carefully restructure the transport plan to relocate deficit regions from inside the restricted Voronoi cells to their exterior (Figure 3 rightt).
7

---

Figure 3: (left) In scale 28, the point b transports mass to a region yp inside its 45-expanded Voronoi cell, and (right) in scale ö, the oo becomes violating deficit region. A consolidating path (ro, ri, b2, r2) can be used to bring this deficit to the exterior of the restricted Voronoi cells.
Each scale of our algorithm consists of two steps. Step updates the weights and the transport plan along set of consolidating paths to restore condition (FI). Step 2 of our algorithm maintains õ-feasibility and condition (F1) as invariants, and carefully raises the weight of the surplus points to
x, enforcing condition (F2). From Lemma 4.1 the resulting transport plan is 26-close. We provide non-trivial proof that these adjustments do not impact the asymptotic execution time of our algorithm.
4.2
Our Algorithm
We begin by describing the basic terminology required to describe our algorithm.
Decomposition of A. Given d-feasible transport plan т, y()), let X's denote the arrangement of the restricted Voronoi cell (b) and the â- and 26-expanded restricted Voronoi cells RV (b) and 2V (b), for all points b E B, namely, the subdivision of IR induced by overlaying these cells. See Figure right... For each region E Ks. pick an arbitrary representative point rhp inside and assign it a mass of pr, := #(p) J&p p(a) da. Let Ag denote the set of representative points of all regions Xs. We refer to r E As as a deficit point if the corresponding region yp E Xg deficit. The point rip is violating if it is deficit and p(rp)) > 0.
store and represent a ò-feasible semi-discrete transport plan t, y()) using (discrete) transport plan f over As x B, where for any point b E B and any region yp E Xs. we set f(r,, 6) := (lp,b) Note that any ô-feasible transport plan f, v()) over As x B can be converted back to a d-feasible semi-discrete transport plan that arbitrarily transports a mass of f(r,,. b) from b to yp for any pair (1p, b) E X3 B. We refer to / as an implicit representation of
In our algorithm,
Residual Graph. Given a 8-feasible transport plan T, y-)), define a residual graph Go over the point set As U B, where for each pair (rip. b) As x B, we add forward edge directed from b to rip if rig E RV22(()) (b). Additionally, if ( b) > 0, we add backward edge directed from r.p b. For any triple (b,, b2) formed by forward edge followed by backward edge, the triple (bi, Tep b2) is admissible if (rp, b1) < (r.p b2). Any directed path P in the residual graph is an alternating path. An alternating path P is admissible if all triples (b,, r2. bb+1) E P are admissible. We update along P with a mass of /3 by increasing (resp. decreasing) the mass transported along each forward (resp. backward) edge of P by 3. When P is a path from a surplus point b E B to a deficit point ri E Ag, we refer to P as an augmenting path and to updating i along P as the augment process.
Consolidating paths. Any alternating path P = (ro, Ti....... b rk) in the residual graph is called a consolidating path if ro is a point with ((ro) = and Thk is a deficit point with ((rx) > 0, i.c., P is path from point outside of all 26-expanded cells to a violating deficit point. The consolidating path P is admissible all triples (b,, ri, b,+1) E P are admissible. Note that updating & along P with mass of increases the total mass transported to the violating deficit point r by 3 (potentially making it balanced) and decreases the amount of mass transported to the zero-weight point ro by B.
The scaling algorithm. Our algorithm works in O(log( Ae- scales, where A is the diameter of B. Each scale is associated with a parameter 8 >> and computes a 8-feasible X-robust transport plan. Our algorithm maintains set of weights y(-) for B and uses few subroutines that
8

---

are described in the appendix. At the beginning of the first scale, set d = A and y(b) 0 for all b E B. Execute the following steps while ò > e/2.
1. Removing violating deficit points: For any point b E B, set Tz(a, b) (a, b) for each point (b) and Ta(a, b) for any point a & RV (b). Compute Go and ts with respect to Ta, y()) While there exists violating deficit point E As:
(i) Execute the SEARCCHANDCOOSSOLIDNET procedure, which computes set of admissible augmenting and consolidating paths in 98 and updates i. along these paths. At the end of this step, there are no admissible augmenting and consolidating paths to the violating deficit points in the residual graph.
(ii) Execute the REDUCE WEigHTS procedure, which reduces the weights of all points of B that have admissible paths to the violating deficit points in Aş by 5 and recomputes G3 and řz. iiii) Execute the ACYCLIFY procedure, which updates fa and Gs so that the transport plan řiz is forest and the residual graph does not have any admissible cycles.
Raising surplus weights to X: Set all points b B with y(b) < À as active. While there are active surplus points in B:
(i) Execute the SEARCHANDAUGMENN procedure, which computes set of admissible aug- menting paths and admissible alternating paths from surplus active points to inactive points in G3 and updates îg along these paths. At the end of this step, there are no admissible augmenting paths in the residual graph.
(ii) Exccute the INCREASEW WEIGHTS procedure, which increases the weights of all active points of B which have admissible paths from the active surplus points by 6, marks any point 6 E B with (b) = À as inactive, and recomputes G5 and ř..
(iii) Execute the ACYCLIFY procedure, which updates fa and Go so that the transport plan fg is forest and the residual graph does not have any admissible cycles.
3. Scale Update: Set +- 8/2.
The SEARCIAANDAUMMETT adaptations of the procedures outlined in [I] and similarly maintain the 8-feasibility conditions. The following lemma states the properties of the two new procedures.
EERREASEWEIGITES and ACYCLIFY procedures are straightforward
Lemma 4.2. During the execulion of the SEARCHAND and REDUCEWEIGHTS procedure, the transport plan remains ô-feasible, each balanced point b B remains balanced, the weight of each point E B containing deficit regions inside RV, decreases by 5, and the weight of each surplus point remains unchanged.
4.3 Analysis
Correctness. In each scale 8, our algorithm begins with ô-feasible transport plan, restores the condition (F1) in the step 1, and restores (F2) while maintaining the 8-feasibility condition and (FI) step 2. Since both (FI) and (F2) are satisfied, from Lemma the computed transport plan at the end of each scale 5 is 28-close X-ROT plan.
Efficiency. As shown in the appendix, the residual graph contains O(n²) nodes and (nn²) edges in R Using this. we show that each iteration of steps I and 2 of our algorithm takes (m22((++> log n)) time, where & denotes the query time of an oracle that computes the continuous mass μu inside query triangle. We also show that both steps run ((n) iterations. Summing over all O(log( A/e)) scales, we get total running time of (n 3 ( + n log n) log(A/s) A/))) for our algorithm, as claimed. In the remainder of this section, provide sketch of the proof of the number of iterations.
Suppose V28(-) denotes the of weights of B maintained by our algorithm at the end of scale 26. For any iteration of step 1, let T', (-) denote the 5-feasible transport plan maintained by our algorithm after the ith iteration and define (b) = y2s y² (b).
Lemma 4.3. For any subset S C B, suppose in an iteration of step of our algorithm, the reduction in the weight of any point in S is more than G6 greater than the reduction in the weight of any point B A S, Le., minges (b) > maxyeBy (V)) + 68. Then, there are no deficit regions inside the restricted Voronoi cells of the points in S.

---

Iterations
path and cycle lengths
Error
Regions
Step
Step 2
Step
Step
10
0,009
1
4
41
361
198
20
0.018
2
4
320
991
529
30
0.021
3
5
688
1719
862
40
0.015
5
979
2385
1212
50
0.018
4
5
1603
3061
1544
60
0.015
6
1973
3673
1855
70
0.017
5
6
2512
4129
2130
80
0.016
6
7
2845
4785
2421
Table The results of our experiments empirically verifying our theoretical claims.
Assume that for some iteration i >> 6n, there still remains some violating deficit regions. Using Lemma 4.2 we show that for each surplus point b, E B, the weight of b, remains unchanged during step (i.e., (b.) = and for each point b B with deficit regions in inside (b), the weight of b decreases by in all iterations (i.e., (b) 16). Therefore, for 6n,
compute set B with deficit regions in RV,(b), b, É S for each surplus point
where E S for all points b, B, and that satisfies the conditions of Lemma 4.3 However, from Lemma 4.3 there are no deficit regions in the restricted Voronoi cells of points in S, which is contradiction to our assumption. Hence, the number of iterations of step (n)) Using a similar approach, we show that the number of iterations of step 2 of our algorithm is also O(n). We conclude with the following theorem.
Theorem 4.4. Let μ be continuous distribution with compact support A C Rd. discrete distribution with support of n points S Rd, and > O a parameter. Suppose there exists an oracle that returns the mass of inside any query triangle in t time. Then, there is algorithm that computes an E-close X-robust optimal transport plan in n O(d) log A/e time, where A is the diameter of AUU B.
Experimental evaluation. We implemented a prototype of our algorithm in Python to empirically verify its performance. The implementation omits two ffficiency-criticee components: (i) the compu- tation of continuous mass within regions defined by restricted Voronoi cells, which we approximate using a fine uniform grid over the domain by summing the mass grid points lying inside each region, and (ii) the dynamic tree data structure, which we simulate by naively iterating over all edges due to the lack of efficient Python implementations. We evaluated our implementation on Gaussian distribution with mean [0.5, 0.5] and covariance 0.15/ bounded within the unit square, and added 10% noise to the bottom-left corner by mixing with an exponential distribution with rate parameter 3. The discrete distribution consists of n samples drawn from the same Gaussian with an additional 10% noise sampled from an exponential distribution with rate parameter 0.2 and 0.02, and varied from 10 80. For each value of n, program 10 times and reported the averaged statistics Table I We observe that as n increases, the additive error remains within the target threshold 0.02. Moreover, the number of iterations in the two main steps of the algorithm consistently stays below the theoretical upper bounds (6n and 12n; see Lemma and E.9). Finally, the total lengths of paths and cycles computed in these steps grow sub-quadratically in n., matching the asymptotic behavior predicted by our complexity analysis (Section 4.3). Additional experimental results are provided in Section for the sake of space. The Partial and Sobbust_SDOT
the top-right the
corner. We fixed
code is available at htttps//gghhbb.oom//uuuua//////t
5
Discussion
In this paper, we present new characterizations of the semi cemi-discrete partial and robust OT problems, and design approximation algorithms based on these formulations. Our techniques are stated for -p costs and extend to any cost function whose weighted Voronoi diagram has bounded complexity. The presented algorithms assume access to an oracle that, a polyhedron of constant complexity, com- putes the continuous mass contained in P in & time. For many natural distributions- with sample access or simple closed-form densities-on one can typically assume (poly log n). Achieving polynomial dependence d and 1/e remains interesting direction for future work.
those
10

---

References
[1] Pankaj Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Keegan Yao. A combinatorial algorithm for the semi-discrete optimal transport problem. Advances Neural Information Processing Systems, 37:25857-25887, 2024.
[2] Pankaj K. Agarwal, Kyle Fox, Debmalya Panigrahi, Kasturi Varadarajan, and Allen Xiao. Faster algorithms for the geometric transportation problem. In Proc. 33rd International Sympo- sium on Computational Geometry, pages 7:1-7:16, 2017.
[3] Pankaj K. Agarwal, Sharath Raghvendra, Pouyan Shirzadian, and Keegan Yao. Fast and accurate approximations of the optimal transport in semi-discrete and discrete settings. Proc. 35th Annual ACM-SIAM Symposium Discrete Algorithms,, pages 4514-4529, 2024.
[4] Jason M Altschuler, Jonathan Niles- Weed, and Austin Stromme, Asymptotics for semidiscrete entropic optimal transport. SIAM Journal of Mathematical Analysis, 54(2): 2022.
[5] Luca Ambrogioni, Umut Guclu, and Marcel van Gerven. Wasserstein variational gradient descent: From semi-discrete optimal transport ensemble variational inference. arXiv preprint arXiv: 1.02827, 2018.
[6] Dongsheng An. Yang Guo, Na Lei, Zhongxuan Luo, Shing-Tung Yau, and Xianfeng Gu. AE- optimal transport. In Proc. 8th
OT: A new generative model based on extended International Conference on Learning Representations, 2019.
[7] Sanjeev Arora, Elad Hazan, and Satyen Kale. The multiplicative weights update method: A meta-algorithm and applications. Theory of Computing, 8(1):121-164 2012.
[8] Franz Aurenhammer, Friedrich Hoffmann, and Boris Aronov. Minkowski-type theorems and leaat-squares clustering. Algorithmica, 20( 1066666776, 1998.
[9] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numerische Mathematik, 84(3):375-393, 2000.
[10] 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-987, 2021.
[11] David P Bourne, Bernhard Schmitzer, and Benedikt Wirth. Semi-discrete unbalanced optimal transport and quantization. arXiv preprint arXiv:1808.0/ 0696,, 2018.
[12] Luis A Caffarelli and Robert McCann. Free boundaries in optimal transport and Monge- Ampere obstacle problems. Annals of Mathematics, 171 2010.
[13] Paul B Callahan and S Rao Kosaraju. A decomposition of multidimensional point sets with applications to k-nnarest--eiighbors and n-body potential fields. Journal of the ACM, 42(1): 67-90, 1995.
[14] Laetitia Chapel, Mokhtar Alaya, and Gilles Gasso. Partial optimal tranport with applications on unltbble learning. Advances in Neural Information Processing Systems, 33:2903- 2913, 2020.
[15] Rick Chartrand, Brendt Wohlberg, Kevin Vixie, and Erik Bollt. A gradient solutin the Monge-Kantorovich problem. Applied Mathematiaal Sciences, 3(22) 2009.
[16] Li Chen, Rasmus Kyng, Yang P Liu, Richard Peng, Maximilian Probst Gutenberg, and Sushant Sachdeva. Maximum flow and minimum connmum--oottttt flow in almost-linear time. In Proc. 63rd IEEE
Annual Symposium on Foundations of Computer Science, pages 612-623, 2022.
[17] Nicolas Courty, Rémi Flamary, Devis Tuia, and Alain Rakotomamoniy. Optimal transport for domain adaptation. IEFE Transactions on Pattern Analysis and Machine Intelligence, 39(9): 1853-1865, 2016.
11

---

[18] Nicolas Courty, Rémi Flamary, Amaury Habrard, and Alain Rakotomamonjy. Joint distribution optimal transportation for domain adaptation. Advances in Neural Information Processing Systems, 30:3730-3739, 2017.
[19] Gianluca Crippa, Chloé Jimenez, and Aldo Pratelli. Optimum and equilibrium in transport problem with queue penalization effect. Advances in Calculus of Variations, 2(3):207-244 2009.
[20] Fernando De Goes, Katherine Breeden, Victor Ostromoukhov, and Mathieu Desbrun. Blue noise through optimal transport. ACM Transactions of Graphics, 31(6):1-11, 2012.
[21] Alessio Figalli. The optimal partial transport problem. Archive for Rational Mechanics and Analysis, 195(2):533-560, 2010.
[22] Emily Fox and Jiashuai Lu. deterministic near-linear time approximation scheme for geomet- ric transportation. In Proc. 64th IEEE Annual Symposium Foundations of Computer Science, pages 1301-1315, 2023.
[23] Darius Geiß, Rolf Klein, Rainer Penninger, and Günter Rote. Optimally solving transportation problem using Voronoi diagrams. Computational Geometry, 1608) 2013.
[24] Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach. Stochastic optimization for large-scale optimal transport. Advances in Neural Information Processing Systems, 29: 3432-3440, 2016.
[25] Sariel Har-peled. Geometric Approximation Algorithms. American Mathematical Society, 2011.
[26] Sariel Har-Peled, Piotr Indyk, and Rajeev Motwani. Approximate nearest neighbor: Towards removing the curse of dimensionality. Theory of Computing, 8(1):321-350, 2012.
[27] 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-22 28, 2023.
[28] Andrey Boris Khesin, Aleksandar Nikolov, and Dmitry Paramonov. Preconditioning for the geometric transportation problem. Proc. 35th International Symposium on Computational Geometry, pages 15:1-15:14 2019.
[29] Jun Kitagawa. An iterative scheme for solving the optimal transportation problem. Calculus of Variations and Partial Differential Equations, 51(1):243-263, 2014
[30] Jun Kitagawa, Quentin Mérigot, and Boris Thibert. Convergence of a Newton algorithm for semi-discrete optimal transport. Journal of European Math. Society, 219)):263--2611 2019.
[31] Arthur Leclaire and Julien Rabin. 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):282-308, 2021.
[32] Bruno Lévy and Erica L Schwindt. Notions of optimal transport theory and how to implement them on computer. Computers & Graphics, 1488 2018.
[33] Yiming Ma, Hang Liu, Davide La Vecchia, and Metthieu Lerasle. Inference via robust optimal transportation: theory and methods. International Statistical Review, 2025.
[34] Debarghya Mukherjee, Aritra Guha, Justin M Solomon, Yuckai Sun, and Mikhail Yurochkin. Outlier-robust optimal transport. In Proc. 38th International Conference on Machine Learning, pages 7850-7860 2021.
[35] Sloan Nietert, Ziv Goldfeld, and Rachel Cummings. Outlier-robust optimal transport: Du- ality, structure, and statistical analysis. In Proc. 25th International Conference on Artificial Intelligence and Statistics, pages 11691-11719, 2022.
[36] Sloan Nietert, Rachel Cummings, and Ziv Goldfeld. Robust estimation under the wasserstein distance. arXiv preprint
2023.
12

---

0= [37] Vladimir Oliker and Laird D Prussner. On the numerical solution of the equation 8x By²
8 Uzoy
= and its discretizations, I. Numerische Mathematik, 54(3):271-293, 1989.
[38]
James Orlin. A faster strongly polynomial minimum cost flow algorithm. In Proc. 20th Annual ACM Symposium on Theory Computing, pages 377-387, 1988.
[39] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1 2021.
[40] Abhijeet Phatak, Sharath Raghvendra, Chittaranjan Tripathy, and Kaiyi Zhang. Computing all optimal partial transports. In Proc. 11th International Conference on Learning Representations, 2022.
[41] Hongxing Qin, Yi Chen, Jinlong He, and Baoquan Chen. Wasserstein blue noise sampling. ACM Transactions Graphics, 36(5):1-13, 2017.
[42] Sharath Raghvendra, Pouyan Shirzadian, and Kaiyi Zhang. A new robust partial p- Wassersttinn- based metric for comparing distributions. In Proc. 41st International Conference on Machine Learning, pages 41 41867-41885, 2024.
[43] Lyle Ramshaw and Robert E. Tarjan. A weight-scaling algorithm for min-cost imperfect matchings in bipartite graphs. In Proc. 53rd IEEE Annual Symposium on Foundations of Computer Science, pages 581-590, 2012.
[44] Jonah Sherman. Generalized preconditioning and undirected minimum-cost flow, In Proc. 28th Annual ACM-SIAM Symposium Discrete Algorithms, pages 772-780, 2017.
[45] Daniel D. Sleator and Robert E. Tarjan. data structure for dynamic trees. In Proc. 13th
Anmual ACM Symposium on Theory of Computing, pages 114-122, 1981.
[46] Marc van Kreveld, Frank Staals, Amir Vaxman, and Jordi Vermeulen. Approximating the earth mover's distance between sets geometric objects. arXiv preprint arXiv:2104.08136 2021.
[47] Goran Zuzic. simple boosting framework for transshipment. In 31st Annual European Symposium on Algorithms, pages 104-1, 2023.
13

---

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: We provided the novel characterizations of a-SDOPT and A-SDROT in Section 2 the algorithmic connection between the two problems in Section 3 the high- precision algorithm in Section 4 and the near-linear time approximation algorithm in the appendix.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Our algorithms require the existence of an oracle that returns the amount of continuous mass inside any given triangle. This assumption stated clearly in Theorem [4.4 and in Section Another limitation is the exponential dependence of the running time of our algorithm in Section on the dimension, which we overcame by presenting algorithm (with a lower precision) our final contribution. The exponential dependence on the dimension and the assumption of the existence of the oracle are also present previous work.
linear
3. Theory assumptions and proof'
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: All of the lemmas and theorems are stated including the full set assumptions, and the complete proof of all lemmas is included in the appendix.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of 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: [Yes].
Justification: The experimental setup discussed in Appendix B
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]..
Justification:
As mentioned in Section 4.3 the code is available at https://gihbbb.omm/ooyya/////ftt Partial and Robust SDOT
Robust_SDOT
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: [Yes]..
Justification: The experimental setup discussed in Appendix
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: the paper does not include results that require statistical analysis.
14

---

8. Experiments compute resources
Question:
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?
Answer: [NA].
Justification: the experimental sections does not include any resource specific results, such as running times.
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics /tticsssttt////eerrr http///rrrrpppeeoooooo/o//aaacaa
Answer: [Yes]
Justification: the research conducted the paper does not involve human participants, does not include any datasets, and there are known potential harmful consequences.
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA].
Justification: None of the categories described by the NeurlPS Code of Ethics, namely, safety, security, discrimination, surveillance, deception and harassment, environment, human rights, and bias and fairness is impacted by the research conducted in this paper.
11. Safeguards
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)?
Question:
Answer: [NA].
Justification: the paper poses such risks.
12. 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 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 is 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 participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA].
Justification: the paper does not involve crowdsourcing nor research with human subjects. 15. Institutional review board (IRB) approvals 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 an equivalent approval/review based on the requirements of your country or institution) were obtained?
15

---

Answer: [NA]
Justification: the paper does not involve crowdsourcing nor research with human subjects.
16. Declaration LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA].
Justification: the method development in this research does not involve LLMs.
16

---

A Fast Approximation Algorithms for Robust Optimal Transport
In this appendix, we describe algorithms that compute an (E,, -apppr--ppoo A-robust transport plan between discrete distribution v and continuous distribution μa in near-linear time with respect to the support size of the discrete distribution v.We keep the description of the algorithm at high level, and describe technical details in Appendix Our low-precision algorithms can be summarized with the following high-level sequence of steps.
T
T
Figure 4: (left) Cover of A by hypercubes, and (right) graph formed by two trees with cross edges
1. Discretize continuous distribution. First, we decompose the continuous domain A into a near- linear size set P of hypercubes (see Figure 4. Each hypercube U is guaranteed to satisfy one of two properties: (1) is within distance + /|R - bll to one of the input points b E and has sidelength roughly e² minyysb yb ||b' b ). or (2) the sidelength of is at most roughly times the minimum distance from its center to every point in B. We query the mass of each hypercube and contract this mass to a single point. If the hypercube very close to a point of B, we contract its mass to b. Otherwise, support B also has near linear size.
contract the mass to its center. This induces a discrete distribution whose
2. Reduce to X-capped transport. Mukherjee et al. 34 showed that for the discrete optimal transport problem, computing a X-robust OT plan is equivalent to computing complete transport plan under capped distance function (a, b) == min (d(a,))) b), A}, which denote transport plan. We compute (e, (a.00)aapp -ppprooxie X-capped transport plan i between ja and v using some discrete OT solver and discard all mass transported at cost more than X. The best choice of near-linear time algorithm for X-capped OT depends on the choice of cost function d; we provide efficient algorithms under squared Euclidean distance (Section F3 and é, metrics (Section F.4).
X-capped
(a) Discrete X-capped p- Wasserstein. We construct an E-well separated pair decomposition (s- WSPD) [13 to approximate the set of n² edges with a set of O ") pairs of regions in Rd. We then construct two quadtrees, an up-tree Tr and a down-tree TJ, and set of cross edges connecting each WSPD region in the up-tree to each WSPD region in the down-tree (see Figure[p. This graph satisfies the property that every a in the support of j has a unique directed path to each 6 in B. We cap the edge costs in natural manner and guarantee the path length
from to b is at most (1 + (a) b). then run minimum-cost flow algorithm [I6] on the
directed graph and shortcut the resulting flow using standard techniques.
(b) Faster X-capped 1-Wasserstein. A modification to the graph construction above allows for faster, specialized machinery for the mininum-com flow problem via the multiplicative weights method [47]. First, we employ a random shift on the quadtree to incur bounded distortion in expectation. Instead of making cross edge between each WSPD pair, we take each WSPD pair and add a cross edge between the lowest ancestors whose sidelengths are at least Jog & times their tree distance. Then cap the edge costs as above. One can argue in standard fashion that the shortest path metric de satisfies E (d b)] < (1 + =)d b) for all a, b in the supports of f and v. We make all edges undirected and run the multiplicative weights method.
Unlike prior works multiplicative weight update algorithms for the geometric transportation problem [3] 12 28), our new approach uses a well-separated pair decomposition to inform what cross edges are good. This type of technique is necessary if one wants to simultaneously (1) bound the number of edges by o(d) log A, (2) approximate the capped Euclidean distance with a shortest path metric, and (3) argue the shortest path metric embeds into a tree metric with deterministic Olog A) distortion.
17

---

3. Extract semi-discrete plan. Finally, we design simple scheme to disperse the discrete transport plan from the contracted point throughout its corresponding region. The scheme disperses + ( 6) uniformly throughout each hypercube [1 that is farther than distance ming *b ||b - from b, where en denotes its center. By conditions (1) and (2) on the set of hypercubes, we argue that every point E O is approximately equidistant from b, and therefore any arbitrary dispersion gives comparable cost. Additionally for each hypercube U within distance miny œe 1 bll from 6, we greedily assign the (approximate) closest ((, b) mass from within small ball of to b. The remainder is uniformly distributed among the remaining b' E B.
We make the following conclusion about the algorithm.
Theorem A.1. Let I be a continuous distribution with compact support A C. R v be a discrete distribution with support B S IRd X> O be a parameter and d. p 2 1 be constants. Suppose Q1 is the time complexity to compute Sg u(a)da for any hypercube O [ A and Q2 is the time complexity to compute, given a point b E IR' and constant 2 0, the radius r 2 0 for which the Euclidean ball B (b, r) of radius r centered at satisfies SESD.r) u(a)da c. Then (e, apppprxxiitee X- robust semi-discrete transport plan can be computed in ne- log loge-1 (no(1) log A + Q1 + Q2) time if d(a, b) = ||a - bllP and in (ne -2d logg" (s-2d-4 log" A log n + Q1 + Q2)) time with probability at least _ if, b) = I||== b/|-
For simplicity of exposition in the introduction, we assumed Q1 Q2 = Olog A). If one can represent the pdf function μ in way that allows for fast integration over balls and hypercubes, then Q1 Q2 = (log A) amortized complexity holds. Some examples of such distributions include histograms and Guassian mixtures.
A.1 Comparison with prior work
Before describing the algorithm in full detail, we make some remarks about our algorithm in comparison with existing work for semi-discrete OT and more generally the geometric transportation problem. The structure of our algorithm has many similarities with existing algorithms for full semi- discrete OT 33 46): we first discretize the Euclidean space by dense sampling of the continuous distribution with hypercubes, then run fast discrete OT solver using this geometric sampling, and finally disperse the mass from the discrete OT plan throughout each hypercube. There are a few major bottlenecks when using their algorithms for A-capped OT that we overcome in this paper.
First, the algorithms of 33 [46] as described only work for p = 1. To reduce semi-discrete A-robust p- p-Wasserstein to the discrete A-robust p- Wasserstein problem for p > 1, one cannot simply perform greedy routing of the mass within an E ball of each point b B as done in 3 46 This is because the triangle inequality does not hold for p > 1, and therefore one cannot charge an e error to the greedy routing. We take slightly different approach, where one greedily routes mass within the E ball of point only the discrete transport plan routes mass from to its e ball. This change sufficient to construct an (e, -ppapoooxiittoon for X-capped p- y-Wasserstein distances when p > 1, as proved in Appendix
Second, existing algorithms for the discrete geometric transport problem (c.g. 2 3 22 consider X-capped distances. Approximating this distance with a shortest path metric in a graph is nontrivial. If one chooses to cap edge costs to guarantee particular path has length at most X, then it necessary to guarantee that other paths in the graph are not shortened by too much. too many edges in the graph have their costs capped, then extraneous paths in the graph may develop undesirably low costs. To circumvent this challenge, we design geometric graphs whose shortest paths provably have a bounded number of edges with capped edge costs. This difference is highlighted in Appendices
28 do not
R3 and
B
Additional Experimental Results
We first provide figures for our experiments described in Section 4.3 showing the shape of the transport plan as well as the untransported parts of the mass of the distributions. We then present results using the Beta distribution.
18

---

2.0
10001
-.00003
125
C--
-C0000
-400004
100002
1003
02
00
3.005
Figure 5: (left) The continuous distribution, and (right) the discrete distribution.
À = 0.12
X = 0.22
x= 0.32
X = 0.42
Figure 6: The parts of the mass of the distributions that are not transported by the X-ROT plan computed by our algorithm.
B.1 Experiments on Gaussian Distribution
As described Section 4.3 in our experiments, we used continuous distribution as a Gaussian distribution bounded inside the unit square that is contaminated with 10% noise in the bottom left corner (See Figure 5 (left)). Our discrete support consists of 0.9n samples from the same Gaussian distribution (the red points in Figure (right)) along with 0.1n samples from exponential distribution centered at the top right corner of the square (the blue points in Figure 5 (right)), each assigned mass of 1/n.
Figure 6 shows the parts of the mass of the continuous and discrete distributions that are not transported by the X-ROT plans computed by our algorithm, the value of À increases from 0.13 to 0.42. Here,
ground distance the squared Euclidean distance. We observe that by increasing the value of X, the X-ROT plan transports the majority of the inlier mass of the distributions and leaves the noise parts of the mass untransported.
We extended our experiments to 3 dimensions using a similar framework and constructed continuous and discrete distributions using the Gaussian distribution inside the unit cube with 10% noise mass added to the different corners of the cube for continuous and discrete distributions. The average statistics of our experiments on different values of are shown Table 2
B.2 Results using
a Mixture of Beta Distributions
mixtures of Beta distributions that induce skewed and dispersed mass. Specifically, the discrete distribution consists of n samples drawn from a Beta(a = 2, 8) distribution, with 3 varying in the range [3, 10). As 3 increases, the distribution becomes increasingly concentrated near [o, 0). The continuous distribution is constructed by contaminating the same base distribution with 15% noise drawn from Beta(a 5, B 2) distribution (a distribution concentrated close to [0.9. 0.9).. For smaller values 3, the two components overlap significantly: for larger values, they become well-separated. Each experiment was repeated 10 times, and report the averaged results.
In this experiment,
T the costs and the percentage of inliers and outliers that are transported by a 0.15-ROT plan computed by algorithm. We observe that as & increases (and the inlier and outlier parts
19

---

Iterations
path and cycle lengths
Error
Regions
Step
Step
Step
Step
10
0,002
2
5
25
454
389
20
0.016
5
6
303
2514
1776
30
0.013
5
6
500
4664
3292
40
0.016
11
8
1587
7057
4807
50
0.028
12
11
2820
10954
6681
60
0.024
16
12
3644
13141
8072
70
0.020
15
10
4735
15310
9008
Table 2: The results of our experiments in dimensions.
% Inliers transported
% Outliers transported
Total mass transported
Cost
Robust cost
3
0.022
0.041
0..80
0.499
0.907
4
0.020
0.041
0,992
0.350
0.896
0,019 5
0,042
0.995
0.245
0.883
6
0.017
0.041
0.998
0.191
0.877
7
0,015
0.042
0.998
0.091
0.862
8
0.013
0.041
0.999
0.085
0,862
9
0,012
0.040
0.999
0.053
0.857
10
0.011
0.040
1.000
0.048
0.857
Table 3: The results our experiments on mixtures of Beta distribution.
become more distant from each other), the amount of outlier mass transported by the 0.15-ROT plan computed by algorithm decreases, while such transport plans transport almost all of the inlier mass for all values of 3.
These experiments confirm that our algorithm remains reliable even when inliers and outliers overlap. In particular, the computed partial and robust transport plans consistently focus on the inlier mass while effectively disregarding much of the outlier mass, especially as the main and noise distributions become more dissimilar. Performance metrics (shown in Table 4 remain consistent with those observed in the Gaussian case.
C Missing Proofs and Details from Section 2
Lemma 2.1 (a-OPT Characterization). Let μ be continuous distribution over compact support A C. Rdd, a discrete distribution exists valid weight function y B + Rz0 with Cover(y) = a such that the transport plan Ty induced by is a-OPT plan. Furthermore, this property holds for any valid weight function with Cover(y) = aa
set B of points in R. and [0, parameter. There
Proof. We rewrite the primal linear optimization problem for the a-partial transport problem:
min
Σ
/ r(a, b)d(a, b)da
BEB
s.t.
(a, b)da S v(b) vb B,
T(a,b) b) S m(a) Va E A,
b
Σ
r(a, b)
BEB
T 720 0.
20

---

Iterations
path and cycle lengths
B
Error
Regions
Step
Step 2
Step
Step
3
0,012
4
6
420
993
495
4 0.016
7
8
332
1026
495
0.005 5
6
7
190
986
462
6 0.007
6
8
131
971
435
0.010 7
9
8
229
886
421
8 0.008
6
161
810
382
9 0.007
6
8
89
849
388
10
0.007
3
5
55
766
351
Table 4: The performance of our algorithm in our experiments on mixtures Beta distribution.
Then its corresponding dual linear optimization problem is
max ad A.
Σ (bb)((b) ((a))(a)da
(3)
BEB
s.t. A (b) (a) s d(a,b) A.EA, Φ, v
Note for any fixed value of and o, the objective is maximized when - as large as possible. The largest feasible - one choose is -a) min (minse (dd(a, - X o(b)), 0} for all a E A. After making this substitution, using the definition the restricted Voronoi cells, we obtain the following optimization problem
min min(d(a, - A d(b)), o) BEB
max A.P
at +
Σ BEB 1VV--00
s.t. Φ 2
This problem is concave and unconstrained on \; therefore, it is maximized when 0 (A), where D&: R the above dual objective for fixed values and ÖD, denotes the superderivative of Dg. For values of & where D is smooth, we note
dD (A)
Σ (((R BEB
Cover(A Φ).
Therefore, dD (A) 0 Cover(> - ) -) a. Similar observation extends this to the case when D& dA
not smooth, where one gets E (D) if Cover(() 8) x) S a S Cover((A + 8) Φ) for all 5 > O. For choice of weights (b)
o(b) we conclude Cover = as desired.
Now suppose make the substitution yb) problem 33 to obtain the equivalent dual optimization problem
Φ(b) into the original dual linear optimization
Σ y(b)u(b) - /1 ((a)u(()da - (1-a)A
max X.9.v
BEB
s.t. y(b) ((a) S b(a,b) Va AbbE, (bb) vb e B,
+ 2 0.
By complementary slackness, for any optimal primal solution and dual solution y, Y, if T(a,b) b) > 0 for any a, b, then y(b) 15(a) = b(a,b). But the objective minimized when conclude that if b(a > 0 then E RV, Similarly, if SA T(a, b)da < v(b) then Φ(b) 0 or equivalently y(b) = A- o(b) = A. But by feasibility of the optimal dual solution, we note that À 2 y(b)) for all T(a, b)da v(b) then y(b) maxve y(V/).
-v((a min mnineeB(d(a) b) (b)), as observed above. Hence,
B' B. Therefore,
This proves that y exists, then y satisfies the desired properties. then note the optimal solutions to both primal and dual linear optimization problems exist by strong duality.
21

---

We finally show that for any where Cover(y) and y(b) = max ((6") for all free b E B. the transport plan induced by y is an optimal a-transport plan. Let y be some arbitrary set of weights that satisfy these properties. Define weight function 1 A R that assigns vh(a) == maxse B y(b) - d(a, b) to each point a A. Additionally define 1z} == v(b) - SA ( b)da and (a) == pla) - -ESEB T(a, b) as the free mass of b and a with respect to transport plan т, respectively. We can rewrite the cost of as:
<(T) Σ BEB*
/
d(a, T(a, da
E
//00) (1()) - ((a))
A
DEB
Σ ((b)u(b) / dp(a) Σ
u(b)uz(b)+ + S,, v(a) (17())
BEB
bEB,
Σ v(b) / (p))d((a) BEB
du(a) (1 a)
where Vmax = maxtEB y(b). Let T denote any optimal a-partiil transport plan between μe and v. Then, by the definition of dual weights for points A,
<(7*) = Σ{ /
/ (u(b) - (ba)) da
d(a, b) T" (a, da
E
EEB
Σ (bb)u(b) |1 (1()) du(a)
s(b)rj (b) + S (b)) (a)
Σ
BEB
bEB,
(b))((b) S1 (aaa (u(()) - (1 a N))nnxx
2 Σ
DEB
where the last inequality holds since vb(a) 2 0 and y(b) S ymax- Combining the two bounds,
<(T) = Σ (bb) - / (a) du((a) - (1 a)ymax S c((*).
BEB
Since 7* an a-partial OT plan, we conclude that T is also an a-partial OT plan.
Lemma 2.2 (A-ROT Characterization). Let j. be continuous distribution over compact support A C Rdd discrete distribution over set B in points in Rd and X > 0 parameter. There exists a valid weight function B > Rz0 with Cap(y) = X such that the transport plan induced by is X-ROT plan. Furthermore, this property holds for any valid weight function y with Cap(y) X.
Proof. The primal linear optimization problem for X-robust optimal transport is as follows.
Σ / baab)) da A ( 1 1--(
min
r(a, b) da
BEB
BEB
A
s.t.
b(da sv(b) vb B,
Σ b(a,b) S m(a) A
bEB
2 O.
The corresponding dual linear optimization problem
Σ ((b)((b) -f o(o))(a)da
max 0,
BEB
S.t.
((b) ((a) S d(a, b) E.E B,
6 + 2 0.
22

---

Define y(b) := À @(.). Then, since À is a given parameter and the total mass at the points of B is the dual problem can be rewritten as
Σ ((b)u(b) - - ((a))(a)da - >
max M.V
BEB
y(6) vaa) d(a,b) Vae A,bEB, y(b) - A, Ab B, v 2 O.
By complementary slackness, we observe that if SA b)da < v(6) for some b B then y(6) A. Then one would have y(b) maXNER v(b) and Cap(y) maxsen u(b) X as desired by dual feasibility constraint y(b) S X for
E B. Otherwise if J r(a, b)da = u(b) for all b E B, then let Ymax == maxbe B y(b) and define the new dual weights y/ (b) == y(b) + (A ymmax).
(a) := (aa) + (X Vmnx)- It follows that y' v are feasible and have an equal dual objective value since both ys and v are probability distributions. Furthermore, Cap(y') == max B s X,
as desired.
Additionally note that for any choice of weights y(-) for B, the optimal choice for the function would be vb(a) == max(0, maxbe y(6) - d(a, b)) for each point E A. By complementary slackness, for any pair (a, b) E AA x with T(a. b) > O, we have y(b) (a) d(a, b). Therefore, the point has to be inside (b)(
We argue by strong duality that choice of y optimizing this dual problem exists. Therefore, conclude that there exists a valid weight function satisfying the desired properties that induces a X-robust OT plan.
We finally show that for any y where Cap(y) = and y(6) = maxy EB y(b) for all free E B, the transport plan t induced by is A-robust optimal transport plan. Let y be some arbitrary set of weights that satisfy these properties. Define a weight function vi A + that assigns (aa) max max y(b) d(a.b)} to each point E A. Let
Σ beB SA T T(a, b)
Additionally define v6(6) == vbb) Sa r(a, b) and (7((a) = μ(aa ELEB T(a, b) as the free mass of b and with respect transport plan T, respectively. Using the assumption that is induced by y, we can rewrite the robust cost of as:
(x(T)=
Σ{
/ d(a, d(a.b) +1 11(1II (I-I/
/ B
r(a, b)da
BEB
Σ
(y(b) ((a))T(a (1
Σ (bb)((b) / 5(a) dp(a)
Σ ((b) / JA, +(a)d + aa)-a)
ben
Σ ((b) / A
X(1 a) - A(1 a.)
BEB
/ ((a) du(a).
Σ (b))u())
BEB
Let 7* denote any optimal X-robust optimal transport plan between a and v. Define a" Eren LA T (a, b)da. Then, by the definition of dual weights for points in A.
((a (a, da + > (1 (1-
wx(7*) = Σ
Σ
(aa, b)da
BEB A
BEB A
2 Σ
(y(() - ((a))>>(() da + 2(1 aa")
BEB
A
Σ ((b) / du(a)= Σ
u(b)ut (b) (6)+ / vi(a) duz (a)+ A(1 (1-a")
BEB
BEB,
23

---

Σ v(6)a(b)
S, ((())d(a))) X(1 - a²) X(1 - -a")
BEB
S (())d)(a), du(a),
Σ BEB
v(6)((6)
where the second last line above holds since vb(a) 2 0 and y(b) ymax. Combining the two bounds,
(T) = Ev(b)u(b) - (ba) du(aa s wx(+").
BEB
Since x" a X-robust OT plan, we conclude that r is also a X-robust OT plan.
Lemma 2.3. The mappings MaxCover: Rz0 - [0, 1] and MinCap: [0, 1] + Rz0 are monotoni-
cally non-decreasing functions.
Proof. Consider any two values 0 < λ < A2. Define ai := (lxC (resp. :=
MaxCover(M)) (A2)) and let 7i (resp. 72) denote the corresponding 01 ar-buutst (resp. aegroobust) OT plan. Since is an optimal λy Ar-rottst OT plan,
<(T)) + -1 01) S <((()) + +1 - a2))
Consequently,
- x(T2) (a1 a2))1.
(4)
Similarly, since is optimal A2-robust OT plan,
((T2) + a1)a <(()) + (1-a1)
Therefore,
<(T) <(72) (ar 022)22
(5)
Combining Equations 4 and (),
(0) o2) x2) 2 0.
In words, since Å₁ < A2, then an S 02. This proves MaxCover is an increasing function.
We now show MinCap is increasing function. Let 0 S on < 02 S 1 be two constants. Define A == MinCap(a;) and x2 := MinCap(ag). Let ys denote the corresponding A;-capped weight functions for each a;, i E {1,2} known to exist by definition of MinCap. Additionally let T2 denote the corresponding Ay-robust transport plans for E {1, induced by weights y by
Lemma 2.2 Then we analogously conclude by Lemma 2.2 that each 7; is optimal A3-robust
transport plan, and therefore f(71) <(T2) S (o1 02 022) and d(T)) ((T2) 2 (a1 - a2) Combining these inequalities again gives (a1 (2) X2) Therefore if an a2 then MinCap(a,) < MinCap(a2).
D Missing Proofs and Details from Section
3
We give the subroutine to return transport plan with precisely a mass in this section, then allocate the rest of the section to proving various lemmas and theorems from the main text.
D.1 Routing Exactly Mass
In the binary search procedure as described from the main text, we obtain two transport plans 7y and which satisfy M(TL) I and M(TR) 2 From these two transport plans, we require simple procedure return precisely a mass, mentioned in the main text.
We assume it is possible to return a linear combination of two transport plans. Then return the linear - c)TR + CT₁, for the constant > 0 which guarantees 7PT transports a mass. interpolation T,4' M(ra)-o only depends on M(TL). M(TR) and a, and assumes M(TL) S a The choice of c = M(rz)-M(rL)
and M(TR) 2 This invariant indeed maintained by construction of the binary search procedure.

---

We note that if one possesses stronger oracles on the distributions TAn and μ. then it is possible to greedily remove mass from Tn until Tn transports a mass while guaranteeing comparable precision. Let f, denote the transport plan defined by fr b) [a(a, b) for all a, such that d(a, b) S r and (r b) otherwise. If one can compute the value r* 2 where f,- routes precisely mass, then we can return f.- and guarantee comparable approximation.
D.2 Proofs of Lemmas
Lemma 3.1. Suppose is an arbitrary
Eu (ar apppoxxiiitt X-robust transport plan and 5 E (0, 1).
(i) For any a E (0, if M(T) $ (1 - 5)a, then S (1 - MMinCap(a) + Sa.
(i)) Equivalently, for any a E (0, 1 - 6). ((T) & a, then (1 aa MinCap
1-6
+
Proof. We prove result (i). Then (ii) follows from substituting a' = 1-8 into the bound from (i) and noting 1-8 1.
Let at JA EoEB T(a, b) be the total mass routed by the partial transport plan Then since is a (Er, Ea appproxmmaee X-robust transport plan, we note that
<(T) + (1 a, = (1 + e,) (e(T2) + (1 - a) + Ea
where T& is an optimal a-partial transport plan. To get a bound on the cost transforming (1+E,) into Tā we apply some algebra on the above incquality and observe
<(T)
(a
ar E+
Er {
Ea
++
Er
1 1
1
a+
Now let rà denote optimal utr-partil transport plan. By optimality of Tá we first note
<(Tā,) c((T I+ E S 0.
Then by Lemma 2.3 well optimality of r. (a, b) and ra, (a, b), observe that
c(T
c(Tāa)
MinCap(a).
1 s
e
E+
1+
Finally, by Lemma 2.1 we observe
c(Tx)
<(T.))
<
MinCap(a).
E
1 Er Er
Adding the three above inequalities gives
<(T.))
c(()
Ur 1 eter
MinCap(a).
I I+
Hence,
ar
Er
Eo
S
a
ar
MinCap(a).
It Er The result follows from dividing this inequality by assumption of the Lemma statement, and isolating Lemma 3.2. If An (1+e,)A, + then e(TL) 2 [M(Tn) - M(TL)) Ar
1 + Er
1+ E,
1 It
a- Iter
which bounded above by da by
Proof. Since Txx is a (E,, =.) (a approxximte ro.uustst OT plan, we observe
<(TA.) + (1 ((TAA))) S + Additionally, since TAa is (er, Ea (a approxximat ro-robtst OT plan, we observe {(TAR) + (1 M(TAa))AR < (1+E,)AR Ea,
+1+ (e(TAR) + M((xa)))el) + Ea.
25

---

and therefore e(TAa) - (1 + Er AR + Eu since (Tnn) - 1. Conclude
<(TAR) - x(TA.) 2 MM(AR) - M(TAa) - M(TAM dy - Er((TAR) 6a
2 M((AR] M(TAL) - Ar. e, [11 + +rr)A + Ea) Ea
M [M(TAR) -((TA.) -EolAs -E,,AR - 2Ea
Now we use the assumption AR 2 & ta, observation -1 -1 S -(1 + ex and bound 0 s M(TAa) M(Tx2) S 1 to obtain the inequality
(M((Ta)) M(TA.)) c A 2 (1+ ((M(TAR) -M(xA.) - 6r) An Er)c.
z [M((AR)) (M(TAA)))M(MT)))) - 28r) An 28a
Combining inequalities gives the desired result.
Theorem 3.3. Let p. be continuous distribution with compact support A C RRd. and let be discrete distribution with support B _ Tid and minimum pairwise distance Define A to be the diameter of A U B and let E (0,1) be a parameter. Suppose Er, 2 0 and ART is an arbitrary (Er,Eal -aprrooxxnnntooon algorithm for the X-robust OT problem. For any ó E (0, min { ( 1-a 2)). the algorithm ApT is an (Er, -p-approximmatin algorithm for a-partial transport, where 1J = O (5r MinaaaCa + 8)) + Ea)). For & = O, Apr is an (0, (eaa)-apprximmaaiii algorithm. The algorithm ApT makes O (log A) calls to ART.
Proof. By the condition AR = O initially and end criteria AR S (1 + EE)AA + Eas we Aι know that ApT terminates within Olog A) iterations. Additionally by construction of the plan TAPT. M(TAAT) = SA EbEB TApT (a, da a as desired. Therefore, suffices bound the cost of TAPT
Since TAR is
An-robtst OT plan, observe c(TAR) (1 - M(Txa))AR I (1 E,) (<(7") ) + +1 aa)a +
(6)
M(ran)-a Using ((fx))-M(Tx.))
Now we note that the precise choice of constant c in algorithm ArT is e := the inequality from Lemma 3.2 and multiplying by
observe that
c <(<(TAz) ((Txa)) (a M(TA,)) ))RR S 4E,AR 4Ea.
By definition, we note TArr = Tan + c - -An). Adding Inequalities 66) and gives c(TApT) + (1- o)An = (((xx)) + ce((TT), e((xx)) + a
< <1+ + [1++,) aa) a)+ 4 Ag + 5Ea = =(1+e))((7") + [1 - a- - E,a -5Ey] An + SEa-
Subtract (1 a))n from both sides of the inequality to get
c(TAPT) < (1 +,) <(7") + (5s, E, AR + 3Ea
(1 <(7*) )+++)+ + 55,AR SEa-
(8)
This immediately implies the result when &r 0. For the remainder of the proof, we assume
E, > 0. Then using Lemma 3.1 we note that A < (1 [=) MinCap (a(1-5)-1) + da for
any 0 < < 1. 18 I, then we can use < (1 5) - 1 + 28 to simplify
dr. < (1 (1-
Er δa
MinCap (a<(1 + 26))
Ea
50
additionally 5 25c, then we use this same trick to obtain
2,,
6
d <
1
MinCap (a(1 + 25)) +
5a
&a
26

---

Substituting into the condition An $ Er + gives
[1+
25r dix MinCap (a((1 + 26)) + Ea
AR S (1 (1+
+ Ea
da
4Er
3Ea
1 +
da MinCap (a(1 + 26)) +
δa
3
MinCap (a(1 + 26)) + 8a
8a
We conclude the Theorem statement after substituting this bound AR into Incquality (8.
Corollary 3.4. Let μ be continuous distribution with compact support A C R and let be discrete distribution with support B _ RR. Suppose Er 2 0. a E (0, 1 -+) is parameter, and ART is arbitrary (Er,Ea)-au (app)aiappooiiition algorithm for the X-robust OT problem. Additionally, suppose e(r") 2 e MinCap(a + Ex) for some constant E (0, 1), where r" is an a-OPT plan. Then, ApT is an (O(E)), approximaaiin cappppppmatt algorithm for a-partial transport.
Proof. By the proof of Theorem note
4 (TAAT) S (r) <((")) + 55,An + 58a-
and
AR MinCap (a((1 + 26)) + 35a
8a
8a
Substituting this bound on An gives
65, MinCap (a(1 + 26)) +6)) +
35a
<(TAPT) S ((+ <(T") + 5E,
+ 55a
óa
8ix
30e
158,ea
= =1++ <(r") +
MinCap (a:11 + 25)) 1-a
+ 5E.-
da
8a
Now choose 5 fr which is indeed within the range (0, min < a 2n
2)) for &, small enough.
Substituting this precise choice of 8 gives
30s MinCap (a (1+2
Er
15
c(TAPT) S (1+ sx) <(7") +
+
+ 5
2a
2a a
(1 E,) <(T") + 60z, MinCap (a Ex) 35E..
Substituting the assumption the Corollary gives the desired result:
((TApT) - (1 + E,) <(T*) 60e, MinCap (a + e,) + 355a
$ (1+ E,) ((T") + 60e, <((7") + 35ea
11 + 60(1 + c) Er c(T")
) + 35E
Theorem 3.5. Let μ be continuous distribution with compact support A C IRd, v a discrete distribution with support B C IR, and a E (0, parameter. Suppose E,, 2 0, E, < ja, and
8 5c. Then algorithm ApT is a (E,, 5E. appprooxinntionnn algorithm for a-partial transport given any arbitrary (aprrxattaoooiiimatti algorithm AgT for the X-robust OT problem.
Proof. By the condition of As = 0 initially and stopping criteria AR S (1+ (1+ + Ea- X
we know that ApT terminates within Ollog A) iterations. Additionally by construction of the plan TAPT M(TApr) := SA Eben TApz (a, b) = (1 8)a as desired. Therefore, suffices to bound the cost of TApr'
Since TAn is an (Er,ea) (Gr,,a)-pppooxt approximate ro-robtst OT plan, we observe
c(TAR) (1 M(TAa))AR I + (1+ (e(7") (1 a)ARl t
(9)
27

---

Now note that the precise choice of constant in algorithm Apr is (((ag)-a(1-0) Using the inequality from Lemma 3.2 and multiplying by S 1, we observe that
M(TAn) M(TAn)-M(TA,))
<(F(TA)) - <(TAa)) (a(1 6) - M(TXa) AR S 48, + 4Ea-
(10)
By definition, we note TAPT TAx + +TTA, -AR) Then adding Inequalities p and 10 gives f(TArT) + (1 a(1 - 6))AB ((TAR) + c(((TA))) ((Tx))) + (1 -aa1 8))AR
< (1 ((++,))((7") + [1+ (1++ a) + 4s,r] AR 5Ea
(1+ =1++,,))(7)) + - a - a a+ + 55,] AR + 5Ea
Subtract (1- a(1 - A) from both sides of the inequality to get ((TApT) ) 5 (1 + e+) <(7") + (5s, (e, + 0)() XR + 5ea S (1 + ex) <((") + 58a,
where the second inequality above true 5 2 Ge- This concludes the Theorem statement.
E Missing Details and Proofs of Section
In this section, we present the missing details and proofs of the lemmas in Section 4
We begin by presenting missing proofs our combinatorial framework in Section E. We then provide the details of the implementation of the procedures in Section E.2 E.6 Finally, we provide the proofs of
correctness and efficiency of our algorithm in Section
E.1 Missing proofs of the combinatorial framework
Lemma 4.1. For any parameter 8 > O, suppose T. T.y(.) denotes a 8-feasible transport plan, and ler À Cap(y). Suppose the following two conditions hold: ( every deficit point A, T(a) O, and (F2) all surplus points b E B, y(b) = X. Then, T is 28-close X-robust transport plan.
Proof. For any point E A, suppose ba E B denotes the weighted nearest neighbor of with respect to weights v(-). i.c., 6, := arg minse (y b). For any transport plan r', let μ²' (resp. s²'') denote the sub-measure of μ (resp. v) that transported by 7', i.c. 1e" (resp. v is the part of the mass of (resp. v) that is transported by r'. Let
and
of :
denote the sub- uub-measures of μ and / that not transported by T', i.c., noj and of are the untrans- ported parts of the mass of μ and respectively. We can rewrite the X-robust cost of the transport plan 7 as
= Σ/ / d(a () b) + 2(1 M(T))
(X(T)= = ((T) + A(1 - -(())= -
BEB
=Σ
()(a,b)-vll y(b) + 5(a))T((a)))) da + X(1 M(T))
DEB
+
Σ u(b)uu (b) - (a)) da.
(11)
BEB
Note that from property (F2), any point E with (b) > (i.e., any surplus point e B) has weight y(b) = >; therefore,
Σ (b) Σ (b)
Σ [b)) = Σ (bb)u(b) X(1 M(T)).
(12)
BEB
bEB
BEB
BEB
Similarly, from property (F1), any point a E A with (7((a > 0 (i.e., any deficit point a E A) has zero weight. Hence,
/ ((a)) = ((a)m(a) da - ((()11((a) da = o(a))(a)
(13)
28

---

Finally, due to the 5-feasibility of T, y(-), the point lies inside the 26- 26-expanded restricted Voronoi cell of b, and therefore, by increasing the weights of b by 26, the point becomes the weighted nearest neighbor of a and also contains a in its dual disc, i.e., m(a) s (u(b) + 26) d(a, b). By rearranging this inequality,
d(a,b) b) v(b) ((a) S 26.
(14)
Plugging Equations
12 and 14 into Equation
/. M(a)'' (a) da + A(1 M(T))
(x(T) = Σ
da +
Σ BEB
y(b) (b)
BEB
A
(d ((a)) + v())u(b) - Sx ((a))(a) da
Σ
BEB" A
BEB
28,(a, b) da + Σ (()) - ((a)p(a) da
BEB"
A
BEB
Σ (bb)u(b) J. M(a))(aada BEB
$ 26
+
(15)
where the last inequality holds since t transports at most 1 unit of mass, i.c. M(T) I 1. Next, let 5* denote any X-robust OT plan. We rewrite the X-robust cost of T*
wx(T*) <(7') + X(1 M(7*))
Σ
d(a,b) (a, X(1 100-MM -M(7))
BEB
Σ DEB
(d(a, b) v(b) v(a)))((((())) (a, da 1(1 101-MMn -M('))
Σ b(b)u'' (b) / "(a)u²" (a)
(16)
BEB
Since À Cap(y).
Σ ybbb)² (b) (b)= - (bb)
Σ ubbby (b) 2 Σ v(b)u(b) - A(1 101--((11))) M(()).
(17)
BER
BEB
Similarly, since all weights 7(-) are positive,
J glo)y' (a) (a) /, a(a)ula)
J, o(a)) (a) S g(a))(a) da.
(18)
Finally, by the definition of the derived weights, for any pair of points (a. E Ax d(a,b) y(b) + ((a)= 2 0.
(19)
Plugging Equations
and into Equation 16
wx(7")
(dy 6) + ((a))a"(aa da + X(1 - M(7"))
BEB
+ Σ u()))²"' (b) - -1 sla)ji" (a) da
BEB
(dg(a, + + (a, b) da + Σ Eb(b)u(b) - v(a)u(a) da
BEB
BEB
(bb) / ((a))(a) da.
(20)
BEB
Combining Equations 15 and (00p.
(bb)v(6) - f ((a)u(a) da S wx((") + 28,
(x(T) < 28 +
and is a 2ô-close X-ROT plan.
29

---

Figure 7: (left) the Voronoi-basse boundaries (blue solid lines) and disk-based boundaries (red dashed lines), (middle) the subset of boundaries of the expansions that lie inside the restricted Voronoi cells, and (right) the subset of boundaries of the expansions that lies outside of all restricted Voronoi cells.
The following lemma shows that for any d-feasible transport plan т, v(-), one interchangeably use the implicit representation & without violating the 5-fessibilityy conditions.
Lemma E.1. For any s-feasible semi- discreee transport plan T, v.),,, the implicit representation
i, y())) is also ô-feasible. Similarly, for any set of weights y(-), any ò-feasible transport plan i over
Ag x B can be converted into a ō-feasible semi-discrete transport plan T, v() over A x B.
Proof. Recall that for any semi-discrete transport plan T between μ and v, the implicit representation of r is defined as transport plan where f(r,,, 6) := S T Ta, b) da for each pair (4p, b) E X, x B. For each region yp E Xgs the representative point ri is inside yp. Furthermore, since RV,2 (b) included in the set of cells used in the construction of Xs. the region is either completely inside /VV22 (b) or is completely outside of it. Hence, due to the 5-feasibility of T.y()). for any pair (a, E A x B such that r(a, 6) > (since a E RV (b)), the region ip E Xx containing a also completely lies inside 2V (b) and rip 2VV22(b). Consequently, for any pair (r, b) E Ag x B with f(r, b) > 0, the point r lies inside the 26-expanded restricted Voronoi cell of 6, and f, v-) is 8-fcasible.
Next, suppose & is 8-feasible transport plan over As X B, i.e., for any pair (r, b) As x Bf f(r, b) > 0, then E RVV2 As discussed above, for any region &p E Xz. if ri E (V)) then oo would be completely inside RV Hence, if we define a transport plan 7 such that y Xss B, then for any pair (a,b) A x B such that r(a, b) 0, then a E RV (b), and T, v-) is d-feasible.
T(p, b) f(rp.b) for each pair (1,,b)
Size of residual graph. Next, we show that the residual graph has O(r)) points and O(n edges for d-dimensional distributions, which will be used in analyzing the running time of different steps of our algorithm.
Fix any set of weights v(-) for B. For any point b E B, for simplicity in presentation, we omit the weights y(-) from notation and let (b) = RV; (b), (b) == RV and RV)(b) = 2V225()) to denote the restricted Voronoi cell, 6-expanded, and 26-expanded Voronoi cells of b respectively. For any cell (b) for b E B and
{0, 1, let (b) (resp. (b) denote the part of the boundary of (b) that is formed by the Voronoi cell V7 (b) (resp. disk with radius y(6) + i6 centered 6). We refer to this subset of boundaries as Voronoi- based (resp. disk-based) boundaries. See Figure lefft..
The following lemma provides an important property of these boundaries, which helps us in bounding the size of the residual graph.
Lemma E.2. Suppose R == UsEB (V,(b) denotes the union (of the interior) of all restricted Voronoi cells. For any point E B and any E [1, all interior points of the foronoi-based boundaries RV (b) lies inside R and all interior points the biske boundaries RV; RV) lies outside R.
30

---

Proof. Consider any interior point p RV2 (b). Since this point is determined by the i ed-exaadd Voronoi cell of b. then d(p, b) < y(6) + 16. Suppose p lies on weighted bisector between b and another point V/ E B. In this case.
d(p, 8') v(b) d(p, b) (b(b) + 16) > O.
Since is a weighted nearest neighbor p and has positive weighted distance to p. the point p lies inside the restricted Voronoi cell of b'. See Figure mmiddle).
Similarly, consider any interior point p E (b). By construction, d(p,b) = y(b) + id. Consider
weight function v/(()) that assigns (b) y(b) + iō and vτ(66) == v(b) for all V' + b. By the definition of restricted Voronoi cells, the weighted nearest neighbor of p with weights v)(() is b, and any other point 8/ E B would have non-positive weighted distance to p. Since the weighted distance of each point 8' * to p is the same with respect to y(-) and vy(+) the point p lies outside of the restricted Voronoi cell of b'. See Figure 7 right).
From Lemma H.2 any intersection point in the arrangement is formed using either d Voronoi-based boundaries or disk-based boundaries. As shown in D. there are ((n)) intersection points between the Voronoi-based boundaries. Furthermore, any subset of d disks intersect in most 2 points. Therefore, the total number of such intersection points is ((()) intersection points, and consequently the number of regions in the arrangement is O(n)). Since the residual graph is a bipartite graph between set As (with size O(n²)) )) and set B (with size n), the number of its edges would be O(nd11), leading to the following lemma.
( Hence, the number of
Lemma E.3. For d-dimensional distributions, for any d 2 2, the residual graph has O(n nodes and (n4+1) edges.
E.2 SEARCHANDCONSOLIIO
Procedure
Define V as the set of all violating deficit points of As. Mark all points of B and all backward edges as unvisited. Let U
denote the set of unvisited points, and at any point during the execution of the algorithm, let Q denote the search path. The
procedure initiates a partial DFS from each point in V in the graph To find admissible consolidating and augmenting paths. For any point r E V, set Q = (r = r1) and execute the following steps until Q becomes empty.
1. Q = (r= bi,bi............
bz).
(a) If b, is surplus point, then the reverse path P (br, T's, bi-1...., 71) is an admissible augmenting path in Go. Augment ř8 along P. If remains a (violating) deficit point, then set Q = (r = r1) and continue. Otherwise, stop the search from r.
(b) If b is not a surplus point, for any unvisited backward edge (b,, in 5.. add r' as
T++1 to Q. If there are no unvisited backward edges from by in Gs, then mark b, as visited and remove b, from U and Q. If Q becomes empty, stop the search from b.
2. If i=Q=ch = (r
rı, b..b....,he bi, ri41), let by arg minseu (r((r+1, D') denote the unvisited point with the minimum weighted distance to r++1-
(a) If g(r++1) = 0, then the reverse path consolidating path in Gs. Update fz along P. If remains a (violating) deficit point, set Q = (r = r1) and continue. Otherwise, stop the search from b.
(re+1, bi, .,, ri) is an admissible
(b) Otherwise, if (b', r++1.b)) b.) is admissible, i.e., dy(r++1,b") < (ri+1.,b) b;), then add D* as b,+1 to Q. (c) Otherwise, remove T1+1 from Q and mark the backward edge (b;, r2411 as visited.
The procedure terminates when all partial search procedures from all points in V are terminated. Lemma E.4. Given 6-feasible transport plan Ta, v(-), after the execution of the SEARCHANDDCNN SOLIDATE procedure,
(SCI) the transport plan remains 8-feasible,
(SC2) each balanced point b E B remains balanced, and
31

---

(SC3) there are no admissible consolidating paths and no admissible augmenting paths to violating deficit regions in the residual graph.
Proof. Suppose To (resp. T/)) denotes the transport plan maintained by the algorithm after (resp. before) the execution of the SEARCHANDCONSOLLIAAI
procedure.
For any edge br) E A B with (x(r; b) > 0, if (r, > O, then by the 5-feasibility of f 5()). forward edge
the point/ lies inside RVZ (b). Otherwise, the edge (r, b)
augmenting or consolidating path P so that after updating the transport plan along P, it has become a backward edge, In this case, there a forward edge from b in the residual graph, which means that E RV (b) and (SC1) holds.
Next, note that any augmenting path P from a surplus point b E B to point E As. Furthermore, any consolidating path P is from a zero-weight point r' E As to a violating deficit point rie E Ag- For any point 6 B other than the two endpoints of path P. the amount of mass transported from b remains unchanged after updating the transport plan along P. Therefore, updating a transport plan along any augmenting consolidating path does not change the amount of mass transported from any balanced point, leading to (SC2).
Finally, to prove property (SC3), we begin by showing the following properties:
(SC4) there are no admissible cycles in the residual graph, and
(SC5) any point E B and any backward edge that are marked as visited does not form any admissible augmenting paths to violating deficit points consolidating paths in the same execution of the SEARCHANDCONSOLIDATT procedure.
Suppose pl ..... pk denotes the sequence of augmenting and consolidating paths computed by the SEARCIJANDCONSOLIDIXT procedure, and for any E [0, kj, let F and G3 denote the transport plan and residual graph maintained after updating the transport plan along Pt, respectively. Define V (resp. U') to be the set of visited (resp. unvisited) points of B with respect to fg.
We prove (SC4) using an induction on i. Since 7 is obtained from the transport plan computed by the ACYCLIFY procedure, there are no cycles of admissible triples in the residual graph, and (SC4) holds for g.. For any i [1, kj), assuming (SC4) on ,
. - show that (SC4) also holds for g'. For any admissible triple (b, r, b') in G' formed after updating # along pi, we show below that the triple (b, r, b') does not form admissible cycles and conclude (SC4) on Gj-
For any visited point b E Vi-1 \ Vi-2 the search from b did not lead to the computation of an augmenting or consolidating path; therefore,, any point 5/ that is reachable from b by an admissible path in the reverse of residual graph gi (and therefore, reachable by our partial DFS procedure from b) would have also been marked as visited. In other words, any point D' E B that has an admissible path to the visited point b in the residual graph G also visited. By this observation, there are no admissible paths from any unvisited point in U to any visited point in Vi-1.
Below, we show that for any newly formed admissible triple (b, (,,r,b) 8'), b E Vi-1 and b' E U Assuming this, there are no admissible triples from any visited point to any unvisited point (in the reversed residual graph), and therefore. the newly formed admissible triples (which are from an unvisited point to a visited point in the reversed residual graph) do not form any admissible cycles and (SC4) holds for 9) as well. Consider any triple (b, 6') that is admissible in G but not in procedure does not change the weights of points,
g/-1. Recall that the SEARCHANDCONSOOLIDAR and therefore, since (b, r, b') is not admissible in 93- (r, b') is not a backward edge in G3 i... (r, (r,6)= E P' as a forward edge. From step 2(b) of the weighted nearest unvisited neighbor of r; therefore, since (r, 6) > d, (r, 5') (from the admissibility of (b, r, b')) and the procedure added 1' to the search path (instead of b), the point b was marked as visited by the procedure. Therefore, for any newly formed admissible triple (b, r', b E Vi-1 and E -i
procedure, b' is the
Next, we show that (SC5) holds. For any visited point b E Vi, as discussed above, all vertices that are reachable from by an admissible path in the reversed residual graph are also visited. Since all surplus points b. E B (resp. points by E B transporting mass to a zero-weight point As) are
32

---

unvisited, by (resp. b.) is not reachable from any visited point b E B by our search algorithm and therefore, the visited points do not participate in an admissible augmenting or consolidating path. Furthermore, the procedure marks backward edge (r,b) b) as visited if, for each admissible triple (V,r, r, b). the point b' is visited. Since the point 1/ cannot be included in an admissible augmenting path to a violating deficit point or consolidating path, the visited backward edge (r; b) also does not form an admissible augmenting path to a violating deficit point or a consolidating path.
For any violating deficit point T E As with respect to řss y/-), the search from The has been terminated since Q was empty; hence, all forward edges incident on r'p are to visited points of B, which from (SC5) do not form admissible nngmentig/coosooidatiinn form any admissible augmenting paths to violating deficit points or consolidating paths during the same execution of the SEARCHANDCONSOLIDATT procedure, leading to (SC3).
paths; therefore, Thp will not
E.3 REDUCEWEIGHTT Procedure
The procedure reduces the weights of all points in B that have admissible alternating paths to some deficit points AA with ū(r) > 0. leading to the formation of new admissible augmenting and consolidating paths to the violating free points in the residual graph. The procedure performs partial DFS from all violating deficit points in the reverse of the residual graph to compute a set IC' of all points of B that have admissible alternating paths to the violating deficit points. The procedure then decreases the weights of all points in I. by 8 and shrinks their restricted Voronoi cells. It finally recomputes As. Gs. and fz. We describe the steps of the DFS below.
Define V c A& as the set of violating deficit points. Given the reversed residual graph To. mark all points b B and all backward edges (b, r) as unvisited and set IC and let U = B denote the set of unvisited points of B. For any point E V, start a partial DFS in the reverse residual graph To by setting Q- == (r r1) as the search path that the algorithm grows. Execute the following steps.
1. If Q- (r = ri....,b))) b,).
(a) If there exists an unvisited backward edge (bs, r") ) in Go. add r* to Q as Ti+1- (b) Otherwise, mark 6, as visited and remove b; from U and Q. Add b, to K.
2. If Q = r........ r.). let 6* == arg minser ev dr,r,, 6").
(a) If (r,, < (ris, b;), i... (b", Tr, b.) is admissible, add 6* as 6,41 to Q.
(b) Otherwise, remove r, from Q and mark (b,, rs) visited.
After all DFS procedures from all points in V terminate, for each point b K, set y(b) + y(b) 6. Lemma F.S. Given a ô-feasible transport plan Ts, v(-) such that there are no admissible consolidating paths and no admissible augmenting path to a violating deficit region, after the execution of the LEDUCEWEIGHTS procedure.
(RI) the transport plan remains 5-feasible,
(R2) the weight of each point b E B containing violating deficit regions inside RV, (b) decreases by 8, and
(R3) the weight of each surplus point b B remains unchanged.
Proof. Let y(-), As (resp. y'(+). Xy A/) denote the set of weights of B, decomposition of A, and the set of representative points before (resp. after) the execution of the REDUCEWEIGHN procedure.
To prove property (RI), first note that for each point b K. there are no zero- eeroht points r E RV such that 78 (r, b) O, since otherwise, there would have been an admissible consolidating path in the residual graph, which is in contrast with Lemma E.4 (SC3). Below, we show that for any pair of points (a, b) E x B such that T3 Ta b) > O, we have E 2V (b), and hence (R1) holds. Let b (resp. (/) denote the weighted nearest neighbor of with respect to weights y(-) (resp. y'(-)) From the d-feasibility of Ta. v(-). dl (a, b) < da by)) + 25. Let (resp. r.)) denote the representative point of the region containing a. As discussed above, ((ra) > and ra lies inside the restricted Voronoi cell of b..
33

---

If E K,
- if E KC, then
dy = d)(a, + dy + 30 dy b2) + 28.
Ta, b) would be admissible (since otherwise b/ K.), and Otherwise, the triple (6/ Ta:
dT(Ta, 6) S (Ta, ba). Hence, dy (a, b) S d, by) + and
(a, = dy b) S dy b.) + 20 dy b'g) 28.
Otherwise, b & K. In this case.
d. b) d,(a, b) - da (a k2) + 20 S dy (a, b',) + 25.
Therefore. for any pair b) with 7a b) > 0, we have a E RV (b). and To, v-) would be d-feasible, i.c., (RI) holds.
Since the BEDUCEWEIGIITY procedure initiates partial DFS from each violating deficit point r E V. all points of B with violating deficit regions in their restricted Voronoi cell will be added to the search path and to the set IC in step 1(b) of the procedure; hence, their weights will be reduced by (leading to (R2)). Furthermore, for any surplus point b E B, since there are admissible augmenting paths from the violating deficit regions (from Lemma E.4), the surplus point will not be added to the search path and hence, cannot be added the set K. leading to (R3).
E.4 SEARCHANDAUGMENT Procedure
The SEARCIIANDAUGMENT procedure initiates a partial DFS from deficit points of As in the reverse of the residual graph to find admissible augmenting paths and admissible alternating paths to inactive points. The procedure then updates the tramsport plan along those paths. The SEARCHANDAUGMES procedure in this paper is similar to the SEARCHÂNDAUGMENT procedure introduced in M1 and the only difference is that when the search path reaches an inactive point of B, i..., a point E with (b) X, the procedure updates the transport plan along the alternating path and restarts the search.
Let D C Ag denote the subset of deficit points of A with respect to Tg. Mark all points of and all backward edges unvisited. Define U B as the set of unvisited points of B. For any point r E initiate a partial DFS Q (r rı) in the reverse residual graph To as the search path that the procedure grows. Execute the following steps.
If Q (r = rı, br,,b,,..,,7 bz),
(a) If is surplus, then the reverse path P = re..... ri) path in Go- Augment is along P. r still a deficit point, set Q = (r vi). Otherwise, stop the current search.
-r,...,r1,)......... is admissible augmenting
(b) Otherwise, if b, is an inactive point (i.e., y(b,) A), then the reverse path P (b,, Ti..... T1) is an admissible alternating path to an inactive point in Gg. Update ig along P. If is still a deficit point, set Q (r ri)- Otherwise, stop the current search.
(c) Otherwise, if there is an unvisited backward edge (b;, r') in To then add r' as T:41 to Q. Otherwise, mark b, as visited and remove b, from and Q.
2. If Q (r = rı, b..,,,.., bb, r++1), let == arg mingeu (r.
(a) If b) < (r(4r,,, by), i.e., (b, r:+1, bi) is admissible, add b as d:41 Q.
(b) Otherwise, remove 1 from Q and mark the backward edge (b;, r241) visited.
The procedure terminates when all partial DFS executions from all points in D terminate.
Lemma E.6. Given ô-feasible transport plan Ta, T&,y() y(-) satisfying condition (FI) in Lemma 4.1 after the execution of the SEARCHANDAUGMENN procedure,
(SAI) the transport plan 78, u(.) remains ô-feasible and (FI) remains satisfies, and
(SA2) there are no admissible cycles, admissible augmenting paths, and admissible alternating paths from inactive points deficit regions in the residual graph.
34

---

Proof. Suppose TO (resp. 7/)) denotes the transport plan maintained by the algorithm after (resp.
before) the execution of the SEARCHANDAUGMENT procedure.
For any edge (r,b) b) E AB x B with is (r) > 0, ff b) > 0, then by the d-feasibility of fzz ((),, the point r lies inside RV22 (b). Otherwise, the edge (r, b) was a forward edge on an augmenting path P (resp. alternating path P' from an inactive point) so that after updating the transport plan along P (resp. P'), it has become backward edge. In this case, there is a forward edge from b to in the residual graph, which means that r & RV (b), i.e., f8, v()) is also d-feasible. Furthermore, since the SEARCHANDAUGMERS procedure does not change the weights, (F1) continues to be true, and I holds.
Suppose Pl ,?..... pk denotes the sequence of admissible augmenting paths and alternating paths to inactive points computed by the AARRCIIANAAGGEET procedure, and for any E [0, ki. let f3 and 9j denote the transport plan and residual graph maintained after updating the transport plan along Pt, respectively. Define V (resp. U*) to be the set of visited (resp. unvisited) points of B with respect to f3..
We use an induction on i to show that there are no admissible cycles in the residual graph. Since ig is obtained from the transport plan computed by the ACYCLIFY procedure, there are no admissible cycles in the residual graph G.- For any E [1, k),, assuming there are no admissible cycles in Gg. g - we show that the same holds for g For any admissible triple (b, r. B') in g' formed after updating 2/ z along Pt, we show below that the triple (b, b)) does not form admissible cycles and conclude that there are no admissible cycles in G3.
For any visited point ò E Vi-1 \ Vi-2, the search from b did not lead to the computation of an augmenting path or an alternating path to an inactive point; therefore, any point 5' that is reachable from 6 by an admissible path in the reverse of residual graph gi (and therefore, is reachable by our partial DFS procedure from b) would have also been marked as visited. In other words, any point 8' 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 ( to any visited point in Vi-1.
Below, we show that for any newly formed admissible triple (b, (6,r,8), B)), b E Vi-1 and 8' E U Assuming this, there are no admissible triples from any visited point to any unvisited point (in the reversed residual graph), and therefore, the newly formed admissible triples (which are from an unvisited point to a visited point in the reversed residual graph) do not form any admissible cycles, as claimed. Consider any triple (b, r, 8') that is admissible in 93 but not in 93- Recall that the SEARCH ANDAUGMENT procedure does not change the weights of points, and therefore, since (b, b') is not admissible in G (r, b') is not a backward edge in G.- i..., (r, b²) E pi forward edge. From step of the SEARCHANDAUGMENN unvisited neighbor of r: therefore, since dy (r. b) > dy (r. b') (from the admissibility of (b, r. 8))) and the procedure added 5' 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, B'), b E Vi-1 and D' E Ur-- To summarize, any newly formed admissible triple is from a visited point to an unvisited point (in the residual graph), while there are no admissible triples from unvisited point to visited points; hence, the newly formed admissible triples do not form admissible cycles, and using the inductive hypothesis, the residual graph remains free of any admissible cycles.
procedure, b' is the weighted nearest
To prove property (SA3), we next the following properties:
(SA4) any point b E B and any backward edge that is marked as visited does not form any admissible augmenting paths or alternating paths to inactive points in the same execution of the SEARCIIANDAUGMER procedure.
For any visited point b E V as discussed above, all vertices that are reachable from b by an admissible path in the reversed residual graph are also visited. Since all inactive points b, E B (resp. points b- E B with pero-weight points r E As in (V22(b)) are unvisited, by (resp. b,) is not reachable from any visited point & B by our search algorithm and therefore, the visited points do not participate in an admissible augmenting path or alternating path to inactive points. Furthermore, the procedure marks backward edge (r, b) as visited if, for each admissible triple (b', r, b), the point
35

---

8' is visited. Since the point d' cannot be included in an admissible augmenting path or alternating path to inactive points, the visited backward edge (r, b) also does not form an admissible augmenting path or alternating path to inactive points.
For any deficit point 'D with respect to F8. v(), the search from has been terminated since Q was empty (otherwise, would have been balanced); hence, all forward edges incident on are to visited points of B, which from (SA4) do not form admissible augmenting paths or alternating paths to inactive points; therefore, r will not form any admissible augmenting paths or alternating paths to inactive points during the same execution of the SEARCHANDCONSOLIDD.E procedure, leading to (SA2).
E,5 INCREASEWEIGHTS Procedure
The
procedure is also similar to the INCREASEWEIGHTT procedure introduced in M.. with the slight difference that it inactivates the points whose weights reach À after increasing their weights. In particular, the INCREASE W EIGHTS procedure computes set K. of points that are reachable from the surplus points by admissible alternating paths and increases the weights of those points by 8. It then marks all points of B whose weights have reached À as inactive and recomputes As. G8. and Fz. For completeness, we include the description of the INCREASEWEIGHTS procedure below.
WEIGHTS
For each point r E As. let (r)) S B denote the set of points b E B with fa > 0. Let S C B denote the set of active surplus points of B. Mark all points b E B and all forward edges (b, r) in the residual graph as unvisited and set K. = Let = B denote the set of unvisited points of B. Initiate a partial DFS from any point b E S in the residual graph Gs by setting Q == (6 = bo) as the search path that the algorithm grows. Execute the following steps.
1. If Q (b- bu, P1, . by),
(a) If there exists an unvisited forward edge (b:, r) in Gs. add r to Q as Ti+1-
(b) Otherwise, mark b, as visited and remove b; from U and Q. Add b, to K.
2. If Q- = (b (b=60 bu .,...., b, r;), let 6* := arg maxye eUrW(r)) dy (ra. b).
(a) If dr b*) > (r,, by), i.e. (bs, r;, b²) is admissible, add d* as d,+1 to Q.
(b) Otherwise, remove r; from Q and mark (b,, r;) as visited.
After all DFS procedures from all points in S terminate, for each point b E K. set v(b) + y(b) 8. We next describe how to recompute the residual graph and the compressed transport plan with respect to the updated weights.
Lemma E.7. Given 8-feasible transport plan Ts, v(-) satisfying condition (FI) in Lemma 4.1 after the execution of the INCREASEWEIGHS procedure,
(IW1) the transport plan fs. v(-) remain 6-feasible and condition (FI) remains satisfies,
(IW2) the weight of each active surplus point E B increases by 8, and
(IW3) the weights of all points b B with deficit regions inside RV2*(b) (b) and all inactive points remain unchunged.
Proof. Let X', As (resp. X' A3) denote the set of weights of B, decomposition of A, and the set of representative points before (resp. after) the execution of the INCREASEWEIGHTS procedure.
To prove property (IW1), we show that for any pair of points (a. b) E A x B such that (, b) > 0, we have a E 2V (b). Let ba (resp. bg) denote the weighted nearest neighbor of a with respect to weights v(-) (resp. y (-)). From the 5-feasibility of Ta, v(-), we have dy b) S dy K.) + 26. Let ra denote the representative point of the region containing
If b E K, then by also is included K., since either b b, or the triple (b,, ray b) would be admissible. Hence, bn remains the weighted nearest neighbor of with respect to weights by and
(a b) b) 8 < dy bw) d da bw) + 28.
36

---

Otherwise, b E K. In this case,
- if 6/ & K, then
(a,,,)= dy(a, R/) + da, + 26.
- otherwise, we would have dy 1n) 2 dr b), since otherwise (Ka ra. b) would be admissible and since 1/2 K, then b É K. well. In this case, since 1/ be- came the weighted nearest neighbor of rar should have To RVV and since dy (Ta b/2) 2 dr( b). we should also have that Ta E (b) and consequently a E RV, (b). Then,
(y = dy(a, b) & dy +.) +0= dy (a, b() + 26.
Therefore, for any pair (a,b) with ra(a, > 0, we have a RV (b), and Ts, v()) would be d-feasible, i.e., (IW1) holds.
Note that the INCREASEWEIGHTY procedure initiates partial DFS from each active surplus point, which terminates when the search path becomes empty; hence, each active surplus point E will be added to the set KC in step I(a) of the procedure, and its weight will be increased by 6, proving (IW2).
Finally, note that from Lemma E.7 there are no admissible augmenting paths or alternating paths from active surplus points to inactive points the residual graph. Hence, the partial DFS procedure from the active surplus points cannot reach any inactive point or any point b E B with deficit regions inside RV26 (b); therefore, all such points would not be in the set IC and their weights remain unchanged, leading to (IW3).
E.6 ACYCLIFY Procedure
Given a d-feasible transport plan Ts: v(-) obtained after the execution of REDUCEWEIGUTS or INCREASEW WEGGHTS procedure, the ACYCLIFY procedure updates the transport plan such that, for the implicit representation ig of Tg. the transport plan is a forest, i.e., the graph containing all edges (r, b) E As B does not contain any edges. Furthermore, the ACYCLIFY procedure makes sure that there are no admissible alternating cycles in the residual graph. This procedure is identical to the ACYCLIFY procedure presented in [1] The procedure first uses a dynamic tree structure Sleator and Tarjan
45 to make the transport plan forest (using phase 1). The procedure then uses partial DFS
to follow the admissible triples and remove any admissible alternating cycles (using phase 2). The procedure finally makes the transport plan a forest again (using phase 1). We summarize the two main phases of the procedure next.
Phase 1: Acyclifying the Transport Plan. To make the transport plan a forest, the procedure
initializes an empty transport plan i and iteratively adds the edges transporting mass in řs to f. The procedure also maintains a dynamic tree data structure and for any new edge, if the edge creates a eycle of transportation, the ACYCLIFY procedure cancels the cycle right away. More precisely, after adding an edge (r1, b1) to the transport plan, if a cycle (b), T1, b2..... bR, Tk) is formed such that for all [1, k] (assuming To = тх), then the procedure increases
F(r., 6.) > 0 and f(b,, ri-1) (resp. decreases) the amount of mass transported on the even-indexed (resp. odd-indexed) edges of this cycle until at least one of the odd-indexed edges gets removed from the transport plan. Using the dynamic tree structure by Sleator and Tarjan [4]. each of the above operations takes Olog n) amortized time; therefore, using Lemma E.3 since the number of edges of the residual graph is (nn this phase takes total of ( ( log n) time.
Phase 2: Acyclifying the Admissible Triples. This step uses partial DFS on the reverse of the residual graph Eo to remove any admissible cycles from the residual graph. For any admissible cycle
C= (bi, ri...., b,, rs, di+1 = bi in Gs. define the boutleneck capacity of C as min f(r, by). The procedure cancels C' by increasing (resp. decreasing) the amount of mass transported along the forward (resp. backward) edges of C' by the bottleneck capacity. The update will remove at least one of the transporting edges, hence removing the admissible cycle from the residual graph.
37

---

Mark all points of B and all backward edges as unvisited. Define U := B as the set of unvisited points. For any unvisited point b B, initiate a partial DFS in the reversed residual graph To by setting Q = (6 = b1) and search as follows until Q becomes empty.
1. If Q- (b1, ri.....,be).
(a) If there exists an unvisited backward edge (b:: rip,, add r as T++1 Q.
(b) Otherwise, mark bs as visited and remove by from Q and U.
2. If Q- (b), .r.................... be, rt), let == arg minseu gu dy (r,,
(a) If dr b) < (ry i... (b, ri, b.) is admissible, if b lies the search path Q as by. then the cycle C = (bg. Ti.r....-. bi, Ti, be41 = = bz) is admissible. Cancel the cycle C' and set Q = (b,, by). Otherwise, add as bi+1 to Q.
i.e.
(b) Otherwise, (b, rs. b,) not admissible. Remove r; from Q and mark the backward edge (bi, r.) in 98 as visited.
Lemma E.8. Given 6-feasible transport plan To, v(-) satisfying condition (F1) in Lemma 4.1 after the execution of the ACYCLIFY procedure,
(AI) the transport plan 7s; v remains ō-feasible and condition (FI) remains satisfies, and (A2) the transport plan ts is a forest and there are admissible cycles in the residual graph.
Proof. We begin by proving (A1). The ACYCLIFY procedure first makes the transportation graph acyelic. Since the edges carrying mass the new transport plan a subset of the original transporta- tion network prior to applying ACYCLIFY for each edge (r, b) E As B that transports mass in the new transport plan, r E IVV22(b) (b), which implies that the transport plan produced after this phase remains d-feasible. In the second phase of the procedure, at any time, all triples in the search path are admissible. Thus, updating the transport plan along any one of such cycles does not violate the ö f-fasiibllity conditions.
Furthermore, since the amount of mass transported to each region remains the same at all steps of the procedure, any point E A, that is deficit in the new transport plan is also deficit with respect to the transport plan prior to this phase of ACYCLIFY procedure. Since the procedure does not change the weights of the points, if (FI) holds prior to the execution of the ACYCLIFY procedure, any deficit region would have a zero weight, and condition (FI) continues to be true, i.e., (Al) holds.
To argue that property (A2) also holds, we aim to demonstrate:
(A3) any point b E B (and any backward edge b)) that marked as visited does not participate in any admissible cycle in the same execution of the ACYCLIFY procedure, and
(A4) the portion of the residual graph consisting of visited points and their adjacent regions does not include any admissible cycles. Moreover, there are no admissible paths from an unvisited point to any visited point.
Assuming (A3) is true, we conclude that no visited point b E B becomes part of an admissible cycle. Since the ACYCLIFY algorithm terminates only when all elements of B have been marked as visited, this implies that no admissible cycles exist in the residual graph at termination. Additionally, since the procedure runs phase once more after the execution of phase 2, the transportation network becomes acyclic, leading to (A2).
Furthermore, as discussed next, property (A3) naturally follows from (A4): For any visited point b E B, all points 1/ E B from which there is admissible path to b must also have already been visited (since the algorithm performs the search in the reversed residual graph). If (A4) holds, then no admissible cycle can be formed solely from visited nodes. Similarly, backward edge (r, 6) is marked as visited only if, for every admissible triple (b', r,, b). the point 6' has already been marked as visited. If t/ does not participate in any admissible cycles, then neither can the triple (b'r b). Accordingly, we prove (A4) below, which proves (A3) and consequently, proves (A2).
We induction show that (A4) holds throughout the second phase. Let C² i [1, k]] denote the sequence of admissible cycles computed by the procedure, and let #8. f....., # denote the sequence
38

---

of transport plans, where rg is the transport plan maintained at the beginning of phase 2, and fi is derived by updating # along C*. Initially, at the start of the second phase, property (A4) clearly holds.
A point b E B is marked as visited when search from b does not lead to the formation of an admissible cycle. This means that for all admissible triples (b', r. b), the point b' and backward edge (r. b) must already be marked as visited. Consequently, all points that can reach b via an admissible path (all points that are reachable from b by our backward DFS) are also visited, and b has no admissible paths to any visited points in the residual graph. Therefore, if (A4) holds prior to marking b as visited, it will remain valid afterward. We next consider the case where cycle C* is eliminated.
Any admissible cycle computed by the procedure is made up only of unvisited points of B. After canceling such cycle, any newly created admissible triple (6, r, b') will consist of visited point b and an unvisited point t'. This means that cycle cancellation can only produce new admissible paths from visited unvisited points. Provided that (A4) is held before canceling C*, canceling C" does not introduce any admissible cycles among visited points. Additionally, since cancellation does not generate any admissible triples directed from unvisited to visited points, there will be no admissible paths from any unvisited point to any visited point. Therefore, property (A4) is preserved after canceling any cycle C
E.7 Missing Proofs of Efficiency
Number of iterations of step 1. We begin by showing that the number of iterations of step 1 of our
algorithm is O()).
Lemma 4.3. For any subset S C B, suppose in an iteration i step I of our algorithm, the reduction in the weight of any point in S more than 6ô greater than the reduction in the weight of any point in B \ S, i.e minses (b) > max (b) (10)+++ + 68. Then, there are no deficit regions inside the restricted Voronoi cells of the points in S.
Proof. Let T25, (-) denote the 26-feasible transport plan computed by our algorithm at the end of scale 26. Recall that 7', (-) denotes the 8-feasible transport plan maintained after iteration of the initialization step and (b) (b) y (b) for each point b E B. Define R Ubes RV2 (b) y
as the union of all 25-expanded restricted Voronoi cells of the points in S.
To prove this lemma, we show that all the continuous mass inside the set R$. is transported to the points in S in the transport plan 725 Given that there are no deficit points with respect to 725. the total continuous mass inside the set RS would be no more than the discrete mass at the points in S. Due to the ò-feasibility conditions, 7¹ transports the mass of the point S only to the regions in Rg. i.e., the total mass transported from S in 7 is no more than vSS)= -bes Etes v(b). Since the REDUCE WEIGHTS procedure does not reduce the weight of any surplus points, no points in can be surplus; hence, points in S can be deficit, and all points in S are balanced. In the following, show that the total continuous mass inside Rgg is at most vSS).
Consider any point a E RS, and let b, B be any 45-weighted nearest neighbor of a with respect to weights (-s(.). Below, we show that b. E S. For the sake of contradiction, suppose b, & S, and let b' E S denote any point in the set S. Note that dges (a, ba) S dyza (a, b') + 45. By the lemma's assumption, (OV) < (baa) + 68, since ba & S. Therefore,
dy (a, b') = dyss (a, b') + (b) (b') 2 dyza (a, by) + 7 - 46
= (a, b..) - (bb.) (6') 46 > dy b..) + 26.
In words, since ba & 5, with respect to weights y (-), the weighted distance of b2 to a would be more than 28 less than the weighted distance of any point D' E S: consequently, the point cannot lie inside the union of the 26 expanded restricted Voronoi cells of S (i.c., & R$). which is a contradiction. Hence, any 45-weighted nearest neighbor of a with respect to weights y25() has be in the set S.
Next, note that any point a E R$ lies inside the 25-expanded Voronoi cell of at least one point in S with respect to weights y²(.). Since all points in S have their weights reduced by at least 65 (from the lemma's s assumption), the derived weight (2)(a) with respect to weights y22(.) would be positive, and from the 26-feasibility of T25, (-)(.), the point a cannot be free point.
39

---

Consequently, all continuous mass inside Ri is transported to the points in S in T26, as claimed.
For any iteration i of step 1, let A'p C As (resp. B'g C B) denote the subset of violating deficit points (resp. surplus points) with respect to r¹. Note that step 1 terminates as soon A becomes empty. Since the total amount of continuous mass transported by r' is no more than the total discrete mass at the points in B, if A + 0. then there exists at least one surplus point b, E B'g. Since step does not create new surplus points (from Lemma E.4 balanced points remain balanced), the point b, was surplus in all previous iterations as well, and 7 (b.) = O. Furthermore, from Lemma 4.2 the weight of the point by E B containing any deficit point E A reduces by 5 in each iteration of step 1, and therefore, (b,,) ) = 6ř. Consequently. if > 6n, then (b,) (bb)) > 6n1,, and there will be a subset S B such that b, E S for all points b, B containing violating deficit point
A'p. b. & S for all surplus points b, E B's. and that satisfies the conditions of Lemma 4.3 In this case, from Lemma 4.3 there are no deficit regions inside the restricted Voronoi cells of the points in S, which is contradiction. Hence, there will be violating deficit points in As after 6n iterations of step 1, i.e., step 1 will terminate after at most 6n iterations.
Number of iterations of step 2.
Next, we show that step 2 of our algorithm in each scale runs ((n) iterations. For any iteration i, support "step2 Visep (-) denote the 5-feasible transport plan maintained after after iteration i of step 2. Define "ihepp
Viiep(()) y22(b))
Lemma E.9. For any subset S C. B, suppose in an iteration i > 12n of step 2. the increase in the weight of any point b E S is more than 66 greater than the increase in the weight of any point E B S, i.e. millies (b) > max Vsupp + 6ò. Then, all points in S are balanced in Supp
Proof. First, recall that step 2 terminates when there are no active surplus points. Hence, we assume that there exists an active surplus point by E B in iteration i. Since each iteration, the weight of by increases by 6, Tsepp (b.) = i5. Also, b, has the highest increase in weight and therefore, b, E S. Note that any point b E B that is surplus with respect to the transport plan Tsp was also surplus point in all previous iterations of steps and 2. Hence, the weight of has not been reduced in step 1. Recall that if was a surplus point in T28. then b has been dcactivated by our algorithm, which means (2s(b) X. However, step has not reduced the weight of b, and therefore, the point b would have been deactivated in step 2, i.e., the surplus point by is not balanced in T26 as is balanced.
Second, note that any surplus point b E B with respect to T25- U26(.) has weight (b) = x, and since the number of iterations of step at most 6n, 7aupp2 (b.) Z Thep (b) + 6nd, and cannot be in S. Therefore, all points b E S are balanced in T28-
Define R$ = Uses RV-26(b) as the union of all 25-shrunken restricted Voronoi cells of the points
in S. From the d-feasibility of Taepp Vasepp (-), the continuous mass inside Rg. is transported only to the point in S. To prove this lemma, we show that the amount of continuous mass inside the set R is no less than the discrete mass at the points in S. Assuming this, since there are no deficit points in B in step 2, no points in S can be surplus, and all points in S are balanced. In the following, we show that the total continuous mass inside R3 is equal to v(S).
Consider any point b E S and any point a E A such that (25 b) > 0. Below, we show that a E Rg. Let b' & S be any other point. Since T26 b) >> 0, then is 46-weighted nearest neighbor of a with weights y2s(.) and therefore, dyzo (a, S dyzs (a, + 48. Furthermore, since D' S, Tssop(()) (b) > Tsspp2 (b') + 65. Therefore,
y (a, b) dyze (a, + (b)(b) Vatepe (b) dyzs (a, b)) + 45 T5mp2
duyu (a, b¹) Nappa (b') + 48 T/ppp(( (b) < Vang (a, b') 28.
In words, with respect weights B (-), the weighted distance of b a would be more than 28 less than the weighted distance of a to any point 1/ & S.i.e., a R3.
Consequently, all continuous mass inside R$ is transported to the points in in T28, as claimed. Since all points in S are balanced in 725, the total discrete mass at the points of S are no more than
40

---

the continuous mass inside R2, and since there are no deficit points in B, there are also no surplus points in S and all points in S are balanced.
From Lemma E.7 the weight of each active surplus point b E B increases by ó in each iteration of step 2. Furthermore, note that for any point 6 B that had a weight of X in the previous scale (e.g., the surplus points of B in T28). the weight of has been reduced in step by at most 6nd, and in step 2. it can increase by at most 6nò. Assume, for the sake of contradiction, that the number of iterations of step 2 has exceeded 12n, and consider any iteration i >> 12n. In this case, for any surplus point by E B, since by was a surplus point in all previous iterations, including the iterations of step 1, we have
(bsp)(bs) - Viuue (b.) v2s(b,) - (Vsnepe (b.) - y (b.)) (U2s(b.) - v (b.) 12nd 0 0 12nd.
Furthermore, for any inactive point 6' B with respect to 725. 1/2(()), after increasing the weight of 6' by at most 6nd,, it becomes inactive; hence, the difference in the increases in the weights of b, and D' is more than 6nō, and there exists subset S of B where all active surplus points are in S and all inactive points are in B S. However, from Lemma E.9 all points in S have to be balanced, which is a contradiction. Therefore, the step 2 of our algorithm terminates in most 12n iterations.
Efficiency of the SEARCHAND execution of the BAARCHANDCONSOLLIIE DFS procedures on the residual graph to find a set of admissible augmenting and consolidating (resp. alternating and augmenting) paths. The partial DFS procedure, upon backtracking from a point b E B (resp. r E As), marks the point b (resp. the backward edge r)) used to reach r) as visited and does not add it to the search path again in the same execution. Upon finding an admissible path P, the procedure updates the transport plan along P in (( time. Let (Pi,-... Pk) denote the set of all paths found by the DEARCHANDCONSOLEE the residual graph has On edges (from Lemma E.3) the rumning time of the procedure would be O(na+1 1 + Ek=1 (P.)). Next, we bound the total length of all paths found by the procedures.
and SEARCHA AUAMMMETT Procedures. Each (resp. BEARCHANDAUGMEERS procedure runs partial
(resp. SEARCHANDAUGMEL
) procedure. Since
Lemma E.10. The lotal length of augmenting and consolidating paths computed during the execution of the SEARCHHANDCONSOLIDAATT procedure is ( in dimensions.
Proof. Let F8 denote the transport plan maintained by the algorithm the beginning of the execution of the SEARCHA COOOSHHAANCOOOODDATT procedure. To prove this lemma, we categorize the augmenting and consolidating paths found by the procedure based on the source of their bottleneck capacity, namely (1) set P., consisting of augmenting and consolidating paths whose bottleneck capacity is determined based on the residual capacity of its endpoints, and (2) set Pa consisting of paths whose bottleneck capacity is determined based on mass transportation over its backward edges. We begin by showing that [P,] = O(n and then show [Pel = O(n)). Since each path has length of at most 2n (since the number of points in B is n),, we then conclude that the total length of all paths computed by the procedure is O(n²).
Let P be a consolidating path in P.,. If the bottleneck capacity of P is determined by its free endpoint r E As.. then the mass of r will be fully transported after updating the transport plan along P: similarly, for any augmenting path P E P.. if the bottleneck capacity of P is determined by its free endpoint b B (resp. r E As), then b (resp. r) will be balanced after augmentation. Therefore, [P,l O(n²) since |As U B O Next, let P be path in Pa; in this case, the backward edge (r, b) determining the bottleneck capacity of P will be removed from the transport plan after augmentation.
Consider any triple (6, b') that is admissible in G' but not in G. Recall that by the definition of the admissible triples, (r, 5') is a backward edge in G' and (r, 6) > (r, 6). Since the procedure does not change the weights y(-). the only case where (b, r, b') is not admissible in G is when (r. is not backward edge in G. i.e., the pair (b', r) is in P as forward edge, and updating à along P results in transporting mass from b' to r. On the other hand, by step 2(b) of the added to the search path only if b' is the weighted nearest unvisited neighbor of r; in other words, since (r, b) > (r: 1/) and the procedure added I/ to the search path (instead of b), the point b was marked as visited by the procedure. Therefore, for any newly formed admissible triple (b, r, b'), point b (resp. 6') is marked as visited (resp. unvisited).
SEARCIANDCONSOLDDA
procedure, a forward edge (b r) will be
41

---

By property (SC2) in Lemma E.4 the point b' cannot form an admissible augmenting path during the same execution of the SEARCHANDCONSOLIDAAN the bottleneck capacity of P was a backward edge of the initial transport plan rg and updating the transport plan along each path P E Pe removes one of the transporting edges of the transport plan ig Since the transport 1g is obtained from the ACYCLIFY procedure, using property (AC2) in Lemma E.8 the transport plan is a forest and the number of backward edges is O(n"); hence, |Pel = O(n)) ) claimed. The total number of augmenting paths found by the procedure, therefore, O(n)), and since each augmenting path has a length of at most 2n, their total length is O(n+1)).
procedure. Hence, the edge (r, b) determining
Using an identical argument, one can show that the total length of all alternating and augmenting paths found by the SEARCIANDAUGMENT procedure is O(n Therefore, each exccution of the SEARGHANDCONSOLIDAR and SEARCHANDAUGMEN procedures takes O(z time.
Efficiency of the REDUCEWEIGHTS and INCREASE WEIGHTS Procedures. The RE- DUCEWEIGHTS and INCREASEWEIGHT procedures run DFS that visits each edge of the residual graph at most once and has a total running time of O d dimensions. Furthermore, in the arrangement used to construct partitioning V, each point E has at most restricted Voronoi cells (three cells that are used in the construction of X3 and three that are used in the construction of X²). Similar to our proof for the size of the graph, one can show that the total number of vertices in the arrangement used to construct Y' is Onn)), and the number of regions in y is at most O(z)). The construction of the transport plan f, therefore, can be done in ( + n)) time since
(1) the mass of all regions in y can be determined using the oracle in ((n time, and (2) the mass transported on each pair (e) E y x B can be determined in ((1) time.
Converting f to fs. as is done in the merge step, also takes (nd+)) time, since there are 1( pairs of points in A, x B, and therefore, the total complexity of 4 is O(n++). Finally, storing a sorted list of neighbors for each region r E As takes 1 log n time in total. Hence, the executions of the REDUCEWEIGHT and INCREASEWEIGHTS procedures take (n + n log n))
time.
Efficiency of the ACYCLIFY Procedure. In the first step, the ACYCLIFY procedure uses a dynamic tree structure to remove all cycles of transportation from the transport plan î.. Using the dynamic tree structure by Sleator and Tarjan 45,, since the total number of edges of the graph is O(n in d dimensions, the running time of this step would be O(rd+1 log n). In the second step, the procedure runs a partial DFS procedure on the residual graph and cancels the admissible cycles. The partial DFS procedure, upon backtracking from point E B (resp. backward edge r)) used to reach r) as visited and does not add it again to the search path in the same execution. Furthermore, upon finding an admissible cycle C. it cancels the cycle in O(CC) time. Let (Ci,..., C& } denote the set of all cycles found in the execution of the ACYCLIFY procedure. Given that the size of the residual graph is at most ((d the second step of the ACYCLIFY procedure takes a total of On + Ei |C,I) time. The following lemma bounds the total length of all cycles found by the ACYCLIFY procedure.
As), marks the point b (resp. the
Lemma E.11. The total length of admissible cycles computed during the execution of the second step of the ACYCLIFY procedure is O(n²)) in 2 dimensions and ( time in d dimensions.
Proof. Let #9 denote the transport plan maintained by the algorithm after the first step of the ACYCLIFY procedure, i.e., 7 is forest. To prove this lemma, we show that the ACYCLIFY procedure finds O(nd) admissible cycles, where each cycle has a length of 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 in Lemma E. 10 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 [F.8, the point does not form an admissible cycle in the same execution of the ACYCLIFY procedure and therefore, the newly formed backward edge
42

---

(r. 8)) cannot be included in any admissible cycles. Therefore, each cycle cancellation removes one of the backward edges of fg. where Fg is forest and the number of its transporting edges is O(n"). Therefore, the total number of cycles found by the ACYCLIFY procedure is O(n"), and their total length is O(n*+)), as claimed.
From Lemma E.11 the total execution time of the ACYCLIFY procedure is O(n*+)) in d dimensions.
F Full Algorithm Details from Section A
We first make some necessary definitions and useful notations for the remainder of the section.
Let B (c, r) denote the Euclidean ball of radius r centered at c, and let G5 = { [0, 5) + Z: 8 8z4) be the set of axis-aligned grid cells of side length 5. For any hypercube D. let Bn == B n denote the set of points of B in D. cn denote the center of D, and En denote the sidelength of D.
Given a hypercube grid cell D I + 8) from Gs, we define the set of its child cells
as C[]] = - a; + L + 5 : T E {0, } Additionally, given a tree T = (V, E),, we use pa(v) and Ciu) to denote the parent of v and the set of children of for each v E V.
We define any pair of sets P. I2d as s-well separated if
Ilp - p2ll, max 119: 92ll S 8 min q/1-
max
max PI. PI.PSEP
qu.4aEQ
p-qEPxQ
Then given a set S C Rd, a e-well separated pair decomposition of S is a set of pairs of subsets W QI),
(Px, Qs)} of S such that (1.) the pair (P:, Q:) is s-well separated for all Isis k, and (2.) for every p; q E S, there exists a unique j such that E Pj and E Qj. It is known that an e-well separated pair decomposition of a set of n points in Rd of size (nn--d) can be constructed in O(n(e- + logn)) n)) time, e.g.
13.
We divide the remainder of this section to describe each of the components of the algorithm in Section A
F1 Construction of discrete distribution
We first describe the construction of the collection of hypercubes, Compute an e-well separated pair decomposition W of the discrete point set B. For each pair (B) B2) E W, choose arbitrary representative points d E Bj and by E Bz. Define & 2'el - byl for all E Z. Then let be the set of axis-aligned grid cells O from Ge& where on B (b), 6,) B (b2, 5,). and
PB, B₂)
initially
(o) U P,(BI, B2). For every O E such that there exists some C also in P, we replace D with its 2d child cells in P. For convenience, we will denote P Uuep as the set of points covered by P. We conclude with set of hypercubes P with the following desired properties.
Lemma F.I. [P| = O(ns 2) and can be constructed in O(n(s 2dd loge 1 + log n)) time. log
Moreover, every O in P satisfies one of the following two properties: (1.) tg < 2e minke B ||b elll
or (2.) b aol || < e min b lbb- b'|1 and fn v E ming
IV -all for some b E B.
Proof. We first prove the bound IP| O 2d logs Note that |W| O(ne d)
and can be constructed in O(n(s- + log n)) time. Then by initial construction of P UB U log P,(Bi, B2) and potential replacement of each hypercube D in P with at most its 2d child cells, we observe that
P 2d "loge-t)
max (P((B,, B2) Bg.3
Then one can argue that [Pr B2)l = O(E-d) for all and all (B1, B₂) E W since at most O (de interior-disjoint hypereubes of sidelength e r can fit inside the ball of radius r. We require ( B2)l) time to find the grid cells in each set P,(Bi, (2).
43

---

We now prove that one of the two properties holds for each hypereube D in P. Let O E P be arbitrarily chosen. Suppose ||6 - --lll
y ||8 bl] for some b B. Let 6* be an arbitrary element of arg ming #b || - bll, and let (P, Q) be the pair in well-separated pair decomposition
W such that b E P and b* E Q- Then must be the case that [ is an element of Pu(P..) Q). which
implies LD I =²||b= - bll. Let D' * b be an arbitrary element of B. Then observe
|lb* - bbll S 1/
by definition of 6*. But by triangle inequality, ||bl - coll 2 | - bll ||b coll 2 (1 - - b1|. Combining gives
E² s - bll S 1-E
11 coll-
We note 6' + b was arbitrarily chosen.
Now suppose there does not exist any b E B such that || - cnll < E ming 1| -bl- Then for every b1. b2 E B, there exists an E 2 such that > 0 and en B (b1, 6,) B (b). 6;-1 where again 5, = 2 |b| b2ll- If there exists a b2 E B such that < 2 log 8 then by construction of P we
have that {g S εō; < 28 |bb roll-
2 log== for all by then cn & B (b1: # 1162 ball). This implies ||6 call > A since by triangle inequality one can argue that for every b1 B there exists a by B with ||| - -b2l 2 A.. But since E P. we note that there must exist a by E B such that
If instead
en E B (bs," maxb libs - ball) and lo I maxbu lbb - ball < A. Conclude that tn S <
2a|l6 - lll.
Once we have constructed the set of hypercubes P. we create a discrete distribution in the following way. For each b E B, define the s-neighborhood of b as the set
(ebb)= (D P: 11b- coll S
min I|b #b
8'113
of hypercubes in P that are much closer to b than any other point in B. We will use the notation N₂ (b) = (b) D to denote the set of points covered by (b) for convenience. Then define j(b) = SN.WS) daa)daa for each b B and (((=)
u(a)da for each hypercube O P UbeB-N(()) Finally choose arbitrary point from IR \ P and define ju(r) = 1 fr pla)da as the leftover mass from μ not covered by P. This completes the construction of the discrete distribution Jt.
F.2 Reduction to >-capped OT
Mukherjee et al. (34] showed that for the discrete optimal transport problem, computing a A-robust OT plan is equivalent to computing a complete transport plan under a capped distance function. We extend the same result to the semi-discrete OT problem.
Given a value A > 0, define the X-capped distance function (+.-)) d' b) == min(d(a, b), A}- For any complete ammi-discrete transport plan yd between and v with respect to the distance function -,),,, we derive partial semi-discrete transport plan т, referred to as X-capped transport plan, by removing all mass that is transported on edges with distance at least X. Formally, we define 7 for each pair of points (a, b) E A x B as r(a, b) (a, b) if d(a, b) < X and r(a, b) = O otherwise. We refer to the transport plan derived from a complete OT plan under the X-capped distances by a X-derived OT plan. The following lemma relates the X-derived OT plans and X-robust OT plans.
Lemma F.2. Given parameter X O. let TX be an (E, appproxximat X-capped transport plan, and let T be a partial transport plan derived from TX Then, T is (c, 0) approximme X-robust OT plan between μ and v.
Proof. Let dx b) == min/d(a, b), A} denote the X-capped distance between and b. Since 7> is an (c, apppproxiee X-capped transport plan, we note that for any arbitrary transport plan TA where M(Tī) = 1,
Σ BEB
(alab)dd b)da S (1+s)
Σ BEB
T (a.ddaa,)ddd
44

---

Now let r" denote the X-derived transport plan from xx i.e.
T (a, b), d(a, b) S otherwise.
T'(a, b) -0, 0,
Then follows from the definition of X-capped distances and derived transport plans that
Σ BEB"
ra b)dx b)da
Σ /
=(a, b)d((,
(1-I)
(1
Σ/ BE
T(a, b)daa =x(())
BEB*
and
/ ri b)dx(a, b)da = Σ (a, b)d(a, b)da + a (1-Σ// T (a, b)da = wx(+"). BE
Σ BEB*
BEB
Therefore, substituting into the first inequality at the beginning of the proof, we conclude
(T) (1 (T")
The result follows from the fact that the transport plan T was chosen arbitrarily among transport plans where M(Tx)) =1.
F33 Algorithm for X-Capped p- Wasserstein distance
To compute approximate discrete X-capped p- p-Wasserstein distances (we equivalently consider d(a,b) a blll and do not take the nd root the plan cost), we use a compressed directed graph construction of 22) cap distances at X, run directed min-cost flow in near-linear time using [16], and shorteut flows to produce transport plan. We describe the construction of the instance of minimum cost flow for completeness, and defer procedures for minimum cost flow [16] and shortcutting the flow to construct transport plan 22.
Sparse directed graph. Suppose we are given two discrete point sets Å, B C Rd of at most m points, with respective distributions and v. We construct an e-well separated pair decomposition W of pairs of quadtree cells from quadtree T on Å U B as in e.g. (25] The well-separated pair decomposition W guarantees that for every pair of cells (01,0g E W, 24. S lo, < 260, We replace each pair (D,, (2) in W where &, 26m, with the set Utrequ, D2) of pairs of O2 with each child E' of D. This guarantees that each pair of cells in W is at the same level of the quadtree T while only increasing the number of pairs in W by factor of 2d
We now use the e-well separated pair decomposition W to construct a sparse graph. Make two copies of the quadtree T = (Vr, ET), called 7i (V1, E1) and T2 (122 E2), where En consists of only directed edges going up the tree and Eg consists of only directed edges going down the tree.
We note that W consists of pairs of quadtree cells (C1, 02) where |An, > 0 and |Bn > by
( |An, 1,
B >
construction. For each (2) E W, let 01
a, Åo = {a} and 0/2 b.
Ba, = {b}
be representative points from Dj and 02, respectively. Let Fgz {m - 02 (0,,0) E W] be the set of directed edges from U1 to 12. Then the sparse graph we construct is G = (Ve, FG), where VG Vi U V22 and EG = Es U Eg U Ez.
Lemma F3. The graph G constructed satisfies IVal Onns log A loge-1) and |Egl O(ne-3a 3d 3doog log A log==1). Additionally, G can be constructed in O(ne- log A llg time.
Proof. The vertex set Va of our graph consists of two copies of centers of all quadtree cells the quadtree T. We note that the quadtree 7 on points with spread (or diameter if minimum pairwise distance is assumed to be 1) has size at most Om log A). Then the quadtree built on the support of f (contains B) has size (nngdd log A log by IP| + B (nns 22 log e-1) by Lemma F.1
-11 since the size of the support of á bounded above
The edge set consists of two copies of each tree edge (one up and one down) plus one cross edge for each pair in the well-separated pair decomposition of the support of f.. Note that the size of
45

---

the support of f is O 2nd 2ogg=" log=="), and therefore its E- E-WSPD has size Onn 3d log==1). By
e.g. (26), the e- G-WSPD can certainly be constructed in (n + log A log e ) time. Finally, the number of tree edges is bounded by the number of vertices, One 2d log E )) Combining
bounds on the number of edges of each type gives the desired result.
Construction of minimum-cost flow instance. As input. the minimum- flow problem is given directed graph G EV) with edge costs de: E R and demands 77: V - R satisfying Evey (7())= = 0. The solution to the minimum-cost flow problem a flow, i.e. a function f: E R2.0 satisfying Emw-wer f(u + v) - -EE f'w + u) = r7(4) for all u E V, minimizing the weighted cost Lu-ver f( + v) dc(u We have constructed our desired sparse graph G. to complete the instance of minimum-coss flow for which we call [16], we require edge costs dg and demands 1).
Now given this sparse graph G. we define edge costs so that the cost of the path from a to is bounded by (1 + O())) min -||", A}. Note that there does not necessarily exist path in G between every u and v in V. However, by construction of the well-separated pair decomposition W and edge set Eg, we guarantee that there exists unique path from every E Á to every b E B. We define a height parameter hy in a bottom-up manner. For each point p E Á AU B, we set the cost of the edge p > pa(p) in Ej and the edge pa(p) - p in Eg as de(p + pa(p)) = dc(pa(p) + p) =0
and additionally define hi (pa(p)) = hp(p) : O. Then for each other vertex v E VT in a bottom-up (x
order, we define (pa(v))) = min ( (x()) + p0)-))l)) Set the cost of the edge v - pa(p) in E1 and pa(v) + u in Eg as
do(v pa(t)) do (pa(v) + v) = hx(pa(u)) hx()).
Note this edge cost dg is always nonnegative since h (pa(v)) 2 (x(v) for all E Vr. Then for each edge u + v E Fg. we define the cost of the edge as
da(u - w) == min( - hn( - h2 (u), 1|u -l|-}-
For each w, 2 E VG such that w + z is not an edge in EG, we let do(w, 2) denote the shortest path distance between and z in G.
Lemma F4. The graph G and shortest path distance de satisfy
(1 22") min ||a b1I" S do(a, b) S (1+ (2 + 4db) e) min(A, 1|a 1laa-1l%% b11"}
for any constant p 2 1 and all E An V12 and b E B V2.
Proof. The structure of the proof will be as follows. We first prove the upper bound on IE' [da 6)]. then prove the lower bound on E [da 6)].
1. Upper bound. Define the canonical path from a to b the path aa b) == m cD, c02 ~TT; b, where u w*y v denotes the unique shortest path from u to using only edges from tree T and (em, <D,) is the edge added E for the WSPD pair (P, (P.Q) where a P and b Q. By construction of the set of cross edges, there exists a unique such edge (D,, 02) for each (P,Q) E W. Moreover, let c(T)= = Ex-ver do(u, v) denote the cost of the path I' in G.
We note (O) S 3/2 for all in the quadtree by construction of h1x Moreover, ((T ((Temna m9 b)) < Ca--vvTT...aa-- |/u vil, then e(T cm ~ b)) = X by construction of dG any path between two leaves with exactly one cross edge. Then
do b) < (a ~a b)) I min
X
Σ
||u u/IP
N N--OET (amob)
To conclude the upper bound on do(a, b). we prove En-veTan
I|u- S (1 + +)llaa - bl].
Consider the one cross edge (2) of path Pem ~ b). We note tion of the cross edges, and by construction of W, one has
([), E W by construc-
||c0, 10, (1 + 25)² ||e - bll < (1 + 23) 1lla- -l1..
46

---

Furthermore,
u a|| - 2d l 2a ||a b|l?
since the length of each vertical edge is exponentially increasing and W is well-separated pair decomposition. Analogously, Ew-rsErD web 1|w =|| I 2d ep ||a bllp. One can then conclude
Σ
1 U1l - (1 (2 + +4d3) 8) 1|a b||-
M. On
Lower bound. We prove the canonical path from a to is the unique path from for any E An V1 and Bn V2. Note to get from a to b, at least one cross edge required. Moreover, in Ti, the point a can only reach vertices en for cells O in the quadtree 'T' where E C. This is true since Ti only has edges going up the tree. Similarly, in T2, the point b can only be reached by vertices for cells D in the quadtree T where b E O. But by assumption, each cross edge corresponds to a pair of cells in the WSPD W. Since there exists a unique pair of cells (O,, in W where 1 and E 22 there must exist unique path in G' from to (going through the WSPD pair (2)).
With similar application as above of the fact that D1, Dl is well-separated pair, we argue that the cost de (a, b) comparable to ||a 6|1 Note da b) = e(T'cam(( b)) by the proof above that the path from a to b is unique. Then by definition of dG on the cross edges, observe
(a b) e(T (Ceaa(a ~ b)) = min X,
||u - u/.
N N--EET can (a-e))
Consider the cross edge (0,, (2) of path Peam (a b). We note tion of the cross edges, and by construction W, one has
(O1, [2) W by construc-
||o, <D2 1l (1 25) ||a bll 2 (1 - 22 .11a bl1..
Then the result follows after Σ Eu-veeTan(a-b) 1/u wIl 2 ||<D. - cD, 11 since the cross edge (O), [2) is an edge of the path I' (a my b).
We now describe the demands for the min-cost flow instance. For each a Å, we take its copy in V1 and assign it demand of na)) fa(a). For each b B, we take its copy in V/2 and assign it demand of (bb = -bb). Then we finally assign a demand of m7(1) for all other v e Vo. This
concludes the construction of the minimum-cos flow instance.
Shortcutting flows. Once minimum cost flow is computed G using the algorithm of (16). employ the procedure described in 2 to shortcut the minimum cost flow on G and obtain a feasible transport plan.
F.4 Algorithm for X-Capped 1- I-Wasserstein distance
The algorithm we use to compute approximate discrete X-capped optimal transport under Euclidean metrics follows similar high-level approach as Section F.3 We first construct a near-linear sized graph whose shortest path distance approximates Euclidean costs in expectation. We then capped the edge lengths so that shortest path distances between points in Á and B approximate the X-capped Euclidean metric. We solve the minimum-cost flow problem on this graph by designing an efficient approximate primal-dual oracle and running the multiplicative weight method using this oracle as 47). Once an approximate solution to minimum-cost flow is computed via multiplicative weights, shortcut paths to obtain transport plan.
47

---

Sparse undirected graph and minimum-cost flow instance. We now construct an undirected sparse graph whose shortest path distance well-appros -ell-approximated by a tree metrie. Our construction is slightly modified from the standard (e.g. 3 (2) to allow for approximate X-capped distances at the expense of a logarithmic increase the number of edges.
Suppose without loss of generality, a,b ||a- bll = 1 and A, B S [- If not, again one can rescale and translate Å and B so that this condition is satisfied. Let D = 1-22', 24/1 + o be randomly shifted grid cell containing Å, B, where is chosen uniformly random from [-4, AJe. We construct a quadtree I recursively with root cell D For a cell O E T, if to Tog N7 then denote D as a cell. Otherwise, we add the child cells D' E C[]] where O' n (Á U B) # 2 to the quadtree I and recurse on each L'. This completes the construction of the tree embedding.
cell.
Let V be the set of cells in quadtree T union with À U B. For each point p E uu B, we add an undirected edge (p, en) to E. where D is the unique leaf cell containing p. We call such edges leaf edges. For every non-leaf cell O and every child cell O' e C[O] of D we add an undirected edge (co. between the centers of the parent-child pair to E. We call these edges free edges or vertical edges.
Now construct an e-well separated pair decomposition W with quadtree cells from T as in Section For each pair of cells (L),, U2) in the WSPD W, let U denote their least common ancestor in T. Define DK as the minimum sidelength cell T containing D that satisfies Ea; Z Tog 20 Analogously
F.3
define D22 as the minimum sidelength cell of I containing D2 that satisfies lmg log & fo We add a single edge (c77, on2) between the centers of these ancestors for each (01, D2) E W, and call these edges cross edges or horizontal edges. Then define the sparse graph G (V, E). The following
properties on G hold.
Lemma F.S. The sparse graph G = (V, F.) satisfies [V] = Onne 2d log A loge and | =
O(ne 3d log A log==1)
Proof. The vertex set V of our graph G consists of centers of all quadtree cells in the quadtree T. We note that the quadtree T on m points with spread (or diameter if minimum pairwise distance is assumed to be 1) A has size at most (m log A). Then the quadtree built on the support of f (contains B3 has size One 2d log A log -1] since the size of the support of f is bounded above by
[P] + |B| Onn 1ine-221oog log 11) by Lemma F.1
The edge set consists of each tree edge plus one cross edge for each pair in the well-separated pair decomposition of the support of ja. Note that the size of the support of ji is O(ne 2d log 8-'). and therefore its g- WSPP has size O(ne 3d log==1). By e.g. (26). the e- E-WS can certainly be constructed in O (e d + log A) loge-1) time. Finally, the number of tree edges is bounded by the number of vertices, O(ne 22d log 1)). Combining bounds on the number of edges of each type gives the desired result.
We detine the cost de of each edge as follows. For each leaf edge (u, v), let du v) = ||u e|]- We will assume À 2 1; otherwise every transport plan has cost of since no pair a, b À x B has
cost less than X. We then construct height function, as in Section F3 For each leaf cell 0 in the quadtree T, define its height as h2 ([) = O. Then for each vertical edge (C0.c)) where D' E C[
define the height of D as hx (D) = min { { (x0)) + ||co - I ||} and the length of edge (aa )
as de crr) = (D) - hx Finally, for each horizontal edge (cT,,,T,)) define its cost as do (CT, ) = min (X- - (x,)
h2 (D2), Ilen, - en |}}. For any pair of vertices u, E V, we denote de [u, v) as the shortest path distance between and v. We additionally define dr v) as the shortest path distance in G' between u and using only leaf edges and vertical edges. That is, du, v) the shortest path distance on the subgraph T of G with respect to edge costs do. Then the following properties of dG and dr hold.
Lemma F. For every a, b E Á B,
min llla - bl1, S E (da b)] < (1 8dE) min { b11,
where expectation is over random choice of ū used to define O*.
48

---

Proof. The structure of the proof will be as follows. We first prove the upper bound on E [da (a, b)], then prove the lower bound on (d.(a,b)) b)].
1. Upper bound. Define the canonical path from a to b as the path c ~y b) = a wr co, - <02 why b, where war denotes the unique shortest path from to using only tree edges and (01, [2) is the edge added to for the WSPD pair (P, Q) where a E and b Q. By construction of the of cross edges, there exists unique such edge (Di, (2) for each (P,Q) E W. Moreover, let c(T) = Eu-ner ( v) denote the cost of the path I in G.
We note hx(1) S for all C in the quadtree T by construction of hA Moreover, if ( (a ~9
b)) Σ (e--ee ||u z/Il then c(T ~ b)) = X + log
2v Σ (1 + log ( À by construction of de on any path between two leaves with exactly one cross edge (leaf edges have cost most vdr and is assumed À 2 1). Then Tog
E [ddaaa, b) S E (e(TT (oTTaam (a ~ b])]
2vde log A
S E min
1
A
Σ
||u-
R M--eETan(a--))
2v log A
min
X,
E
1/u w11
U-T (am+b)
by Jensen's 's Inequality, where expectation is over the random shift of in the root of the quadtree
T. To conclude the upper bound on E [da(a, b)], we prove [u-vvETa(o--)) ||u (1 + +))lal b/1.
wll
Consider the one cross edge (O,, of path ca ~ b). If (O,, E W, then by construction of W we note ||e0, CD2 (1 2s) ||a bl]. Furthermore,
Σ
||u - u|l < 2Vden, I 22Vd= aaallll000 bl|.
and analogously Σ Ew-sseq web ||w - || 2vds ||a bl1. One can then conclude
Σ
||u w|| s 1 + (2 + 4VQ) 8) ||a bll-
M- (am-))))
Now suppose (2) & W. Then by construction of the cross edges, én, = Lr, I log A & where Das is the smallest cell in the quadtree T containing both a and b. In this case, it then follows that
E
Σ
| lu- 211
E
Σ
1|u
U|ll
+ E [[lon, cD,1l
+E
||w z||
0--$E0 C worb
2
E [40.a E [||e, c021l] -0,11]+22 2vd +vd log
E [f0,,
log A
4
E (8 + +E[[ - -c,,
log A
We additionally note that by triangle inequality,
|| .0.1
Σ
||u w/1
+ 11a bll
Σ
w 10-22|
w-
49

---

Therefore, E [llen, 002 I bound
||a
1v/d E [6 Substituting this inequality gives the bll
8vde log
E
| wll
S laa
bll+
E [e
(21)
a E-ttETan(aaa
Now to bound IE (f0,)] we observe that if as is the smallest cell containing both a and b, then the child cell D. of Lab containing a must contain b. Therefore, there must be face on the boundary of D. which intersects the line segment from a to 6. Furthermore, since Dab is convex and contains both a and b, we observe that the entire line segment from a must be contained
in Dob- We use these two properties to claim Pr (k... 2 ddj|a-bll 27 Then summing over all possible lengths of the cell ab in the quadtree T gives the bound on the expectation
log A
E [1
Pr Pr(. 20 -2'
j ==log Tog log A
vālla - bl
23
20+1
j =llog Tog 5
log
Vd
-1aa bll
j j= log
$ (Vdlog A) ||aa b/1.
where the last inequality assumes log
log S log A. Substituting this bound above into Equa-
tion 21 gives
Σ
||m
w|l
S (1 + 8dz). ||a bl1-
LI an
This concludes the proof of the upper bound.
2. Lower bound. Given a path r, define Vert(T)) and Hor(T') as the set of vertical and cross edges in H, respectively. Again define <(T) == Eu-rorr (u, v) to the cost of the path in G. For each vertex of the graph G, let D. denote the cell of the quadtree corresponding to vertex u and let lev(v) := log lo. denote the level of the vertex v in the quadtree.
Let be an arbitrary path from a to b in G. If ((') > Em-wer ||ue u/l, then since is a path from b can the triangle inequality conclude
<(T)
||u wlll
||a
b1l
min
[||a
b11. A].
Now suppose e(T) < > 1/ - v/|. Then there must be at least one edge w - in the path with do(w, 2) < ||l 2|I.
Suppose (w,, 2) is vertical edge. Then since do 2) < ||w- - =|| and (w, 2) is a vertical edge, hD or (D,)
Suppose, without loss of generality, hx(Dw) subpath of I starting from a and ending w. Additionally let I'2 the subpath of starting from w and ending at b. Then the composition of I2 and F2 is equal to I.. Therefore,
2. Let I2 be the
<(T) <(T)) + <TT2).
Additionally, note each subpath can be decomposed into vertical and horizontal edges:
<(T,)= = Σ
d)(u,u)+
Σ
do(u, w) for i E (1,2).
Vert(T)
(Tr(T))
50

---

By construction of W and the set of cross edges from W, lev(u) = lev(v) for every cross edge (u, v) E. Therefore, reach w from I, there must be at least one vertical edge u, - U in (er for each i < lev (m.,) such that le,( > lev(D..,)
i. Likewise, to reach b from w in [,, there must be at least one vertical edge un - rg in (T2 for each i < (ev( such that lev( < lev(D,:) = i. Then by the fixed depth and definition of child cells of the quadtree T, lxx (D) = hx (D") for all pairs of cells where lev(D) lev(D"). Using the definition of de((; v) for vertical edges (u, v), we deduce that
Σ
da(u, N
Σ hx( 1hg 1i
hx(Ow) - 0=
A
2
u- (T)
Note that hy (Da) = 0 by definition, where D. is the leaf cell in T containing a. We similarly observe Ve-soeVer((Ta) do(u, 2 2. Hence, c(T) 2 X.
Now suppose (w,2) is cross edge. Then by definition of dg(w,2) min -000-21p 211, hx( 11 (m])) it must be the case that (w, z) = À - (_(w)) - hx (m2). Let i denote the subpath of T from a to and let 12 denote the subpath of r from to b. Then the plan
is equal to the composition 11 o (w + 2) o T2- Then in the same manner as the previous case, we argue (e-re da(u,v) Z hx hx( and En-neVan(T.)) do (u, v) 2 hx())). It immediately follows that
<(T) 2 + d =) +
2 hD + -A hx(ww) (O.)) + hx( =
To conclude the lower bound, we observe that the path I from a to b was chosen arbitrarily among all possible paths in G.
Lemma E.7. For every u, E V, dr(u,e) So
loo ( c
do(u, v)-
Proof. We use the definitions of (er)(T), Hor(T), ((T) and lev as in the proof of Lemma F:6 With some slight abuse of notation, we will use lev(O,) and lev(v) interchangeably.
well h (O.) and hx(v) for each vertex v of the quadtree T and corresponding cell D.. We first claim for every edge (u, v) E E, the shortest path from u to in G is through the edge (u, v) with cost dg (u, v).
Suppose (u, v)) is a vertical edge. Then lev(u) o lev(u). Suppose I' is the shortest path from u to u in G. Then there exists some edge -> in the path I such that lev(w) lev(u) and lev(2) lev(u). We note, by the definition of the cost dG on vertical edges and equal distance between all parent and child cell centers at fixed level, that c(T) do (w, 2) do (u, v). Therefore, the edge - is also a shortest path from u to v.
Suppose (u, v) is horizontal edge. Then da(u, v) min(A-Axxuu hax u|l by definition. If there exists a path I' from u to in G such that e(T) = Ep dp(p. 9) < ||u =/I. then by triangle inequality there must exist an edge w + 2 on this path such that do(w, 2) < || - z||- Suppose this edge (w, 2) is a vertical edge. Then without loss of generality assume
pa(z). We
note that by definition of dg, if do 2) < ||w 2|| then hx(w) }. Consider the subpaths r, and 22 of from u to w and from w to e, respectively. Then as in the proof of Lemma F.6 note
Σ
do(p.g) +
Σ
do(p.9) 2 nx(w)
(2 --ha(t)) do(u,v).
(ergeVer((TT)
(-rgEV
If instead (w, 2) is cross edge, then do(w, 2) > - hx(w) (x).). Since (w, 2) and (u, v) are both cross edges, lev lev(w) lev(z) and lev(u) = lev(l). If lew(w) S lev(u), then
do 2) A hx(w)) hx(2) 2 A ha(u) hx(2) 2 dc(u,)).
51

---

Else suppose lev(w) > lev(u). Let 11 be the subpath of I from u to and 22 be the subpath of I from v. Then for every j E (ee(()),ev(u))) there exist vertical edges uj - uj E IT and v/ + v E 22 where lev(",) = lev(u)) = < lev(r)) = lev(r)). Then
e(T) - =(T + (w 2) + 2(T2)
Σ
do(p.9) q) + [ hx(w) - +
E
da(p.
D-gEVerr(Ta)
(w))
lev(w)
da(ug, 12) + - -x(ww) - hx(z)] +
da(ul,
j=)ev(w)
j=lov(u)
(ev(w)
lev(ar)
Σ (hx(ws) - lo (u²)) + [A- hx(w) - [xx[=)]+++++ +
Σ
(1a((1)) - -xx())))(((())
(=))v(u)
j=lev(u)
(hx(w) - (u)) + - h2 - hx (x)] (bx(2) hx
= [A - hx(u) - (x) 2 do v).
This concludes the claim that for every edge (u, v) E E, the shortest path from to in G is through the edge (u, v) with cost (u, v)..
After the above claim, it suffices to prove du 2) S 0 1) -do(u, v) for every cross edge (u, v). Let (u, v) be an arbitrary cross edge. Suppose do (u. = - hx(w)= - -xx Then observe
dr(u,v)
E
(hx(z) hx(w))
Σ
(hx(w) hx))))
w-- lca(u,))
= =(hx(caa(u 0)) hx(u)) + (Aa(Ica(u, v)) hx(w))
2 A hx(4)= hx(v) = do(u, v).
2
If instead dg(u, v) = ||u U|1. then observe
dr(u, z
Σ
1|w=x11
Σ
1| - =||
10 E- lea(u,,)))
10 tn--kEuwallaauu
2Vdlog A
||u -|| 2Vdlog A
do (u,
by construction that for each cross edge u, v, with corresponding cells D, and D, in the quadtree. fo,
to, 2 log - This concludes the desired statement of the Lemma.
We then define the demands y V by n(a) ja(a) for all E nbb) v(b) for all b E B, and n(() for all e E V (ÁUB). This completes the construction of the minimum-cost flow instance we solve.
Approximate minimum-cost flow oracle. We deseribe a simple O ) -approximatee primal- dual oracle for minimum-com cost flow on G. The oracle, given an arbitrary demand function 77. will compute a primal-dual pair (f.y) satisfying the following properties:
f is feasible primal solution to the min-cost flow problem, i.e. 0 and E Emtan))E flu
f(v + u) n7(1) for all u E V,
2. y is O
log A
approximately feasible, i.e. v(u) y(v) < O log do(u, v) for all
u, v E E, and
3. the dual objective value dominates the primal objective, i.e. Evev v(v)(()) 2
E( [uu v)dc(u,
Then it is known from existing works, e.g. above, one can apply the multiplicative weight update method obtain a (1 + e) -aapprooxeaaee minimum-cost flow routing the demand 17.
44 that given an oracle satisfying the properties
52

---

Lemma E8. [147] Lemma 11] Given an oracle which returns primal-duat pair satisfying conditions (1-3) above upon receiving a demand function n, V + R. there exists an algorithm which computes a (1+ -ap-pooipiiit mininum-cost flow after O log n log² A) calls the oracle plus O(EE) operations per oracle call.
The oracle builds a primal-dual solution on the quadtree T in greedy manner. For vertex v E V, let T(v) denote the subtree T rooted at v and define ((T(v)) = EueT(r) nu)). Then for every vertical edge (u, v) E E where u = pa(u), set f(u + v) = ([(()) if (T(()) 2 0 and f(v + u) (T(()) if (TT()) < O.
To construct a dual solution, let r denote the root vertex of the tree T and set u(r) = 0. Then in a top-down manner, for every vertical edge (u,v) E E where pa(u), we set y(u) v(r) + v) sgnn((((())) where sgn(z) denotes the sign function.
Lemma E.9. Given an arbitrary demand function y, the oracle described above computes a primal- dual pair (J,y) in O(ne -3d log A log loge=))
E-1) time satisfying
1. Em(aneer f(u - v) S(v uu) --](u))e for all EV,
2. u(u) - u(u) 0 (1oe2) da(u, v). and
3. Ever u(v)n(o) 2 E(ED)E (ʃ(u v) + s(v u)) dd
Proof. 1. By definition of f. have that for every leaf vertex of the tree T. s(u + pa(u)) (pa(u) + u) = (T(())) - v(u).
Now suppose is lonf vertex of (vertices of I and G are the same), and EeitunjEE f(u v) f(v + u) ru) for all u E T(w) and all w E Clu). Then, note that f(w - v) f(v + w) = nTT(w)) for all w E Clul by construction of f. Therefore, to guarantee the desired property, we must have f(v + pa())) - f (pa(v) + 2) = n7(1) + Ewec(e) ((T(w))) Observe that r7(2) + Eweclo) (T())) = (T(v)) by definition of (((())). Then by definition of f, it is indeed true that f(v + pa(v)) - f(a((v) + v) = (T(())
2. By construction of y, we observe that (b(u) - u(v)) S da(u, v) = dr(u, v) for every vertical edge (u, E E in the graph G. Then for every cross edge (w, z) E E, follows that
y(w) v(z)=
Σ
y(u) v(v)
Σ -ITDEE-T
dr(u, v) dr(w;
where w ~T z is the path the tree T from w to z. The result follows from Lemma where
shown that dw =) S O log A do (w, z).
3. We use complementary slackness and strong duality to prove the solution is an optimal primal- dual pair on the tree T. By definition of and y. we note that for every vertical edge (u, v) where is the child of y(u) y(v) do(u,,) v) if f( - u) > 0 (i.e. if (T())) (()) y(u)= = da(u, v) if f(v + u) > 0 (i.c. if ((T(u) < 0). Additionally if f(u + v) > 0 then f(u + u) O. Moreover, f(u + v) only strictly positive on vertical edges by construction of f. Therefore,
0) and
Σ (f(u + v) + f(v + u) do(u,o)= (u.F)E
Σ (u.))EB:((u-t)))
f(u + v)dc(u,
Σ
f(u + v) (u(u) u(u))
(u.))EE:(u--u))0
= E Σ
(f(u + v) - f(v - u)) v(u)
WEV
Σ
v(u)n(u).
wEV
The third equality above follows from rearranging terms and the last equality above follows from the first condition of the Lemma.
53

---

Shortcutting flows. Once a minimum cost flow computed on the geometric graph G, shortcutting the resulting flow on G to obtain an approximate transport plan in roughly linear time in EE] is standard. See e.g. 22 28
F55 Semi-discrete plan from discrete plan.
Now suppose one has computed an (er, appppoxxxiee X-capped discrete transport plan Êx between and v. We show how to transform ix into an (Er, -appprxxie X-capped semi-discrete transport
plan Ta between μl and v. For all O E P UsEB No. (b), define Tx(a, b) = μ((a) falem.b) over all
A(el)
D and b E B. Likewise, let (a, b) = p(a)
fa(m))) for all R \ P and b E B. What
A(r)
remains a procedure route the mass within each cell O of the neighborhoods NG (b).
We repeat the following local routing scheme for each b E B. If fb(,b)) < Abb), then let Ta b') = uua) (x(b.)) for all a E (b) and b' E B. Otherwise, define A())
P:(b) = { E Gea, ||e0 bll < {
and P.(b) = ()) D for all í Z. Let i" be such that [p. p(a)da I (x b) Sp, (1)) (la)dda Since (b, I A(6), we note (b-(b) S. (e()). Let (a = u(a) for all a E (b). Then for each O E (b) such that there exists a L' P.. where O' D. we replace O with its 2d child cells in Pi++1 (b). Initialize an undispersed mass counter (b, b) Fx b) - Sr.- u(a)do. For each D E Pi+4((b) \ (b) in increas- ing order of |le_ b11, set (a, u(a) min min(cgy (b.b). M(jdda) on all E O and subtract
J= Maadda
min ( (b, b), fn {(a)da} from 13 (b. (b b). After assigning (a b) for all b E B and a E Ne)().. we then set rx(a,b) (m(a) (x,a b)) FA()_6
(b.b,b') for all remaining E Ne(b) and 6' +
Lemma F.10. Suppose is continuous distribution with compact support A C I2 and v is discrete distribution with support C. R2d for some constant d 2 1. Let j be the discrete distribution formed from p as in Section F.1 and let d( b) = ||a - bll" for fixed constant p 2 1. If is an (c, 0)- approximate discrete X-capped transport plan between jú and v, then Tx is an (O(e), -yapprooximate semi-discred X-capped transport plan between μ and v.
Proof. Let T" be an arbitrary transport plan between js and v, where M(r*) = 1, and let f* be the transport plan from js to defined by f* (co, b) Sor" and 6 (b, b) = JN_(0) T" (a, d) for all 6 E B. By Lemma F.I we have that all cells U E P which are not contained in the neighborhood No. (b) of point are 2e-well separated from We use this property to find
T" (a. b)da for all b E B and O E P (b)
Σ BEB"
Ta b)d(a, b)da s (1+ de)))
Σ/
Ta b)d(a b)da.
BEB (b(
bEBacA
Now since ix is (Er, 0) approximitee discrete X-capped OT plan between A and v, we find
ΣΣ (a b)d(ab) s s1+
Σ Σ
b)dlaa
AEA
aE
Finally, we again use the property that all cells D E P which are not contained in the neighborhood (b) of point b are e-well separated from b to compare 4* with r". This property along with the fact that d(b, b) = 0 implies
Σ / 7" (a, b)d((a))))d
4* *"aa < (1 + 2vd=)"
BEB A
BEB
- -(1 -11 + 2v
(" b)da.
bE N.(b)
54

---

We combine the three above inequalities to deduce
Σ Σ/ BEB*
Tala,b)de I + (1+2Vde) 28 (1+ Σ
(a, (a,b)a(a b)d(a, b)da
BEB" A
Σ BEB (2)
b)dab)d((a, b)da
(1 + 2Vdsss
2p
Σ (a(ab)d((a b)da
N.(6)
I (1+ (2 dp =)
Σ
/ +" ()d(a, b)da
+
(a,ab)d(a, byda
DEB (b)
(1 + (1+)
Σ
T(a, b)da.
bEB N-(b)
now suffices show that
Σ/e (a (aaa b)d(a, b)da V (1 d8e) BEB (b)
E/
(a, b)d(a, b)da.
BEB N.
If this is true, then we substitute into the above inequality and conclude with the desired result.
By construction of Tx within (b) for each b E B, observe that if (x b) 0. then the cell DE Pia++1 (b) where a E satisfies |lq - bl| S ||c - bll for all E Pie+1 where there exists some a' E D' such that (a', b) p(a'), ). Therefore, for any a, a' Ne(b) where Tx (a, > 0 and
Tx(a', b) < p(a).
1 + vdr
||a - bl| S
2 vde
Ilaa - bll (1+ 2vde) - I|a' - bll
(1+2Vde))
if E S V/ For d(a, b) = ||a - bl1 it then immediately follows that
d(a, b) < (1 [ 2Vde) d(d', b) s (1 (1+29 (1+ 22) d(da,b)
for all such and a', and therefore one can argue
E/e0, (2 b)d(a, b)da s (1+ 2o
E/ ()d,)) b)da.
BEB N.
beB N.(b)
This concludes the statement of the Lemma, since we assume d and p are constants.
F.6 Putting Everything Together
We finally prove Theorem A.1 by simply combining all above lemmas in Appendix
Theorem A.1. Let be a continuous distribution with compact support A C Rd v be discrete distribution with support B C IRd >> be a parameter and d. p be constants. Suppose Q1 is the time complexity to compute Jn u(a)da for any hypercube O S A and Q2 is the time complexity compute, given point b IRd and constant 0, the radius
0 for which the Euclidean ball B (b, r) of radius r centered at satisfies JEO..)) u(a)da = c. Then a (E, appppoxximaat X- robust semiidiscrete transport plan can computed in ne- logg e- (1) log A + Q1 + Q2) time if d(a, 6) = ||a b1le and in O (ne -2d logg" -1 (s-2d-4 log" A log + Q + Q2)) time with probability at least / ifd = ||a bl1-
Proof. By Lemma E.1 we have that the set of hypercubes P satisfies [P] O(ne 2d ==ge=') and
can be constructed in ( logs-1 +log n)) time. Using P, construct the discrete distribution
55

---

ä in ((1P)-Q1) time as in Appendix F.1 Then by Lemmas and F.10 we know that if one computes an (O(E), appprxiaaee X-capped transport plan ix between f and v, an (O()), apprrxxmiaa transport plan тλ between μ and can then be computed in O (IP) + (E- log A Q2) time (we assume i" can be computed in Q2 time for each 6 E B). It therefore suffices to bound the time to construct the (O(=), 0) apprrximmtte X-capped transport plan fx.
For = 1, we note that the graph G Appendix F.4 can be computed in O((E)) = O(ne -3d log A loge time by Lemma F.5 Moreover, by Lemma F66 we note that the min- imum cost flow on G has expected cost within (1 ( factor of the cost of an opti- mal X-capped transport plan between f and v. Then the tree oracle described in Section [.4 takes O(V]) (ns 2d log A loge-)) time. By Lemmas F.8 and F. the tree oracle can be boosted into an (O(e), appproxxiie X-capped transport plan in ( + El)e= log log" A) O(n=-3d-
3d-4 log" A log loge-!) time.
For p 2 1, we note that the graph G in Appendix
F.3 can be computed in O(
O(ne ((ne-adpp log A 1og time Lemma F3 Then the minimum-cost flow on G can be computed in (( the minimum-cost flow on G to form an (O(e).
log A time by the algorithm of [16]. One can then shortcut
0)-approximate x-capped transport plan between jı X-capped transport plan
and v as in [2]. The resulting transport plan an ())=))))-apppooiii
between f and y/ by Lemma
We conclude that an (O(e), O)-approximate X-robust transport plan between μ and can be computed in O(n 00
log A) time for p 2 and in Onn ------d--ppp logg A log n log==! time for p = with probability at least 1/2. To get an (e, -a-aproxxie X-robust transport plan, one can simply choose appropriately smaller value of s.
56


Reviews:
[{"review_id": "xq85Rev33l", "reviewer": "Reviewer_3gQ4", "summary": "This paper presents two categories of semi-discrete OT which are the α-partial transport problem and TV-regularized optimal transport formulation with regularization parameter λ. The two categories can be demonstrated in one common formulation. The theoretical analysis reveals a fundamental connection between these problems, demonstrating that their optimal solutions can be characterized as restricted Laguerre diagrams - a geometric insight that enables us to establish their algorithmic equivalence. Building on this theoretical foundation, novel computational algorithm are developed. A rigorous analysis of the algorithm's precision and computational complexity is presented.", "questions": "1. Could numerical experiments be incorporated to validate the algorithm's efficacy?\n2. These two robust formulations of semi-discrete optimal transport provide mathematically grounded mechanisms. It will enhance noise resistance through distinct regularization paradigms. Could this theoretical claim be empirically validated through some experiments?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n1. Provides a novel characterization of optimal solutions as restricted Laguerre diagrams, deepening the geometric understanding of partial semi optimal transportation with TV regularization. The theory is mathematically rigorous.\n2. Establishes a strong equivalence between α-partial and λ-robust OT, enabling solvers for one problem to be adapted to the other.\n3. The proposed method achieves exact solutions with guaranteed precision, while also providing an efficient approximation scheme, with their respective computational complexities rigorously analyzed. \nWeaknesses:\nThe theoretical findings lack empirical validation and practical substantiation.", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your thoughtful review and comments. We address your concerns below.\n\n**Summary of Contributions.**\n Our paper presents the following key contributions:\n - A characterization of partial and robust variants of the optimal transport (OT) problem in the semi-discrete setting,\n - Theoretical guarantees on the approximation error when using an approximate $\\lambda$-ROT solver to compute an approximate $\\alpha$-OPT plan, and,\n - Two algorithmic solutions for the $\\lambda$-ROT problem.\n\nWhile the original submission emphasized our theoretical contributions, we appreciate your interest in practical implementation and experimental evaluation. In response, we include below preliminary experimental results. A more extensive evaluation will be added in the next version.\n\n**Implementation Details.**\nWe provide a Python implementation of our algorithm. However, it relies on two core procedures that lack efficient implementations in Python:\n1. Estimating continuous mass within regions defined by the arrangement of the restricted Voronoi cells and their expansions.\n2. Maintaining a dynamic tree data structure as per Sleator and Tarjan (reference [44] in the appendix).\n\nTo address these challenges in our prototype:\n- For (1), we approximate the mass within each region using a fine grid, summing the masses of grid points lying inside the region.\n- For (2), we simulate the dynamic tree by naïvely iterating over all edges.\n\nAn efficient implementation will require optimized versions of these components, which remains a technical challenge in Python.\n\n---\n---\n**Preliminary Results on Algorithm Performance.**\nWe evaluated our implementation using a continuous distribution sampled from a Gaussian distribution inside the unit square. The discrete support consists of $n$ samples drawn from the same distribution, each assigned a mass of $1/n$. We set $\\lambda = 0.2$ and $\\varepsilon = 0.02$. Results are averaged over 10 runs for each $n$.\n\n**Key Findings:**\n- The additive error stays within the target $\\varepsilon = 0.02$.\n- The number of iterations in the two steps of our algorithm (Step 1 and 2) remains well below our upper bounds proved in the paper ($6n$ and $12n$ as described in Lemma 4.3 and D.9).\n- The total path and cycle lengths computed in the two steps grow subquadratically in $n$, aligning with our $O(n^2)$ complexity analysis (Section 4.3).\n\nThe details of our experiments are presented in the following table.\n\n| n | Error | Step 1 iterations | Step 2 iterations | Step 1 paths and cycles length | Step 2 paths and cycles length | Regions |\n| :--: | :--: | :--: | :--: | :--: | :--: | :--: |\n| 10 | 0.009 | 1 | 4 | 41 | 361 | 198  |\n| 20 | 0.018 | 2 | 4 | 320 | 991  | 529  |\n| 30 | 0.021 | 3 | 5 | 688 | 1719  | 862  |\n| 40 | 0.015 | 3 | 5 | 979 | 2385  | 1212  |\n| 50 |  0.018 | 4 | 5 | 1603 | 3061  | 1544  |\n| 60 | 0.015 | 5 | 6 | 1973 | 3673  | 1855  |\n| 70 | 0.017 | 5 | 6 | 2512 | 4129  | 2130  |\n| 80 | 0.016 | 6 | 7 | 2845 | 4785  | 2421  |\n\n\n---\n---\n**Robustness of $\\alpha$-OPT.**\n\nWe examined the effect of 10% noise in the continuous distribution (added in the top-right corner) and measured how much noise was transported under various $\\alpha$ values. Results confirm the robustness of $\\alpha$-OPT: for $\\alpha$ < 0.9, almost no noise mass is transported.\n\n\n|Alpha| Cost| \\% inliers transported| \\% noise transported|\n| :--: | :--: | :--: | :--: | \n |0.700| 0.010| 0.777| 0.000|\n |0.725| 0.011| 0.804| 0.000|\n |0.750| 0.012| 0.830| 0.015|\n |0.775| 0.013| 0.859| 0.005|\n |0.800| 0.014| 0.887| 0.005|\n |0.825| 0.016| 0.914| 0.009|\n |0.850| 0.017| 0.942| 0.014|\n |0.875| 0.019| 0.970| 0.005|\n |0.900| 0.022| 0.996| 0.020|\n |0.925| 0.031| 0.999| 0.244|\n |0.950| 0.045| 0.999| 0.496|\n |0.975| 0.060| 0.999| 0.748|\n\n**Robustness of $\\lambda$-ROT.**\nWe conducted similar experiments for $\\lambda$-ROT. As $\\lambda$ increases, total mass transported grows, but noise transport remains limited when the transported mass is under 0.9.\n\n|Lambda| Cost| Cost ROT| \\% inliers transported| \\% noise transported| Mass transported|\n   | :--: | :--: | :--: | :--: | :--: | :--: | \n| 0.05| 0.014| 0.024| 0.896| 0.004| 0.807|\n| 0.10| 0.020| 0.031| 0.980| 0.009| 0.884|\n   | 0.15| 0.021| 0.037| 0.991| 0.012| 0.894|\n| 0.20| 0.023| 0.042| 0.996| 0.046| 0.902|\n| 0.25| 0.024| 0.047| 0.999| 0.073| 0.907|\n| 0.30| 0.023| 0.051| 0.999| 0.064| 0.907|\n| 0.35| 0.025| 0.055| 0.999| 0.141| 0.914|\n| 0.40| 0.027| 0.061| 1.0| 0.147| 0.915|\n| 0.45| 0.028| 0.063| 1.0| 0.212| 0.921|\n| 0.50| 0.034| 0.066| 1.0| 0.349| 0.935|\n| 0.55| 0.041| 0.068| 1.0| 0.506| 0.951|\n\nPlease let us know if any further clarification or experiments would strengthen the paper. We are happy to incorporate additional feedback into the next version of our paper."}, {"review_id": "1nfHgUHYS0", "reviewer": "Reviewer_umB2", "summary": "Paper studies two robust formulations of optimal transport (OT) : $\\alpha$-optimal partial transport ($\\alpha$-OPT), which transports only a fraction of the mass of the two input distributions, and $\\lambda$-robust optimal transport ($\\lambda$-ROT), which penalizes untransported mass via a total variation regularization. While both variants have been well-studied in fully discrete settings, e.g. by Caffarelli and MacCann (ref [9] of the paper), this work proposes to extend them to semi-discrete OT. The authors introduce novel characterizations of optimal solutions using restricted weighted Voronoi diagrams, which captures the points of the Voronoi cells that are within a given distance.\n\nAuthors show that any solver for $\\lambda$-ROT can be leveraged to solve $\\alpha$-OPT, thus unifying the algorithmic treatment of these problems. They propose two algorithmic contributions: (1) a high-precision combinatorial algorithm for semi-discrete $\\lambda$-ROT using a refined cost-scaling approach (2) a near-linear-time approximation algorithm. Theoretical analysis includes correctness proofs, runtime bounds, and approximation guarantees even when only a fraction of the mass is transported.\n\nThe paper does not include empirical experiments as it is primarily theoretical in nature. Authors clearly articulate limitations, such as the exponential dependence on dimension.", "questions": "- Could you provide numerical experiments (even in simple cases) that could validate experimentally your algorithm and show its practical utility?\n- how does this work relates to Bourne et al., 2018?\n- the connection between $\\alpha$-OPT and $\\lambda$-ROT has been thoroughly examined in the work of Caffarelli and MacCann (ref [9] of the paper). Could you clarify how your results relate to their duality findings?", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- The paper addresses an important problem: robust and partial OT in semi-discrete settings, for which no/few prior algorithms exists.\n- The characterizations via restricted Voronoi diagrams are very elegant and provide useful geometric insight.\n- The paper is exceptionally well-crafted, presenting the ideas clearly and accessibly despite the complexity of the subject matter.\n- Although I didn't check all the proofs, I think the algorithm and results are reasonable.\n\nWeaknesses:\n- While the contribution is theoretical in nature, numerical evaluations are not provided.\n- the connection between this work and Bourne et al., 2018, are not discussed\n\nBourne, D. P., Schmitzer, B., & Wirth, B. (2018). Semi-discrete unbalanced optimal transport and quantization. arXiv preprint arXiv:1808.01962.", "was_revised": true, "final_justification": "I keep a positive evaluation of the paper, leaning toward acceptance. I found the paper really clear and it provides new results about robust semi-discrete OT problems. My main concern is about the lack of experimental evaluation, beyond toy experimentation.", "rebuttal": "We appreciate your thoughtful review and comments. Below, we answer your main questions and concerns.\n\n---\n---\n>Could you provide numerical experiments (even in simple cases) that could validate experimentally your algorithm and show its practical utility?\n\n**Summary of Contributions.**\n Our paper presents the following key contributions:\n - A characterization of partial and robust variants of the optimal transport (OT) problem in the semi-discrete setting,\n - Theoretical guarantees on the approximation error when using an approximate $\\lambda$-ROT solver to compute an approximate $\\alpha$-OPT plan, and,\n - Two algorithmic solutions for the $\\lambda$-ROT problem.\n\nWhile the original submission emphasized our theoretical contributions, we appreciate your interest in practical implementation and experimental evaluation. In response, we include below preliminary experimental results. A more extensive evaluation will be added in the next version.\n\n**Implementation Details.**\nWe provide a Python implementation of our algorithm. However, it relies on two core procedures that lack efficient implementations in Python:\n1. Estimating continuous mass within regions defined by the arrangement of the restricted Voronoi cells and their expansions.\n2. Maintaining a dynamic tree data structure as per Sleator and Tarjan (reference [44] in the appendix).\n\nTo address these challenges in our prototype:\n- For (1), we approximate the mass within each region using a fine grid, summing the masses of grid points lying inside the region.\n- For (2), we simulate the dynamic tree by naïvely iterating over all edges.\n\nAn efficient implementation will require optimized versions of these components, which remains a technical challenge in Python.\n\n---\n**Preliminary Results on Algorithm Performance.**\nWe evaluated our implementation using a continuous distribution sampled from a Gaussian distribution inside the unit square. The discrete support consists of $n$ samples drawn from the same distribution, each assigned a mass of $1/n$. We set $\\lambda = 0.2$ and $\\varepsilon = 0.02$. Results are averaged over 10 runs for each $n$.\n\n**Key Findings:**\n- The additive error stays within the target $\\varepsilon = 0.02$.\n- The number of iterations in the two steps of our algorithm (Step 1 and 2) remains well below our upper bounds proved in the paper ($6n$ and $12n$ as described in Lemma 4.3 and D.9).\n- The total path and cycle lengths computed in the two steps grow subquadratically in $n$, aligning with our $O(n^2)$ complexity analysis (Section 4.3).\n\nThe details of our experiments are presented in the following table.\n\n| n | Error | Step 1 iterations | Step 2 iterations | Step 1 paths and cycles length | Step 2 paths and cycles length | Regions |\n| :--: | :--: | :--: | :--: | :--: | :--: | :--: |\n| 10 | 0.009 | 1 | 4 | 41 | 361 | 198  |\n| 20 | 0.018 | 2 | 4 | 320 | 991  | 529  |\n| 30 | 0.021 | 3 | 5 | 688 | 1719  | 862  |\n| 40 | 0.015 | 3 | 5 | 979 | 2385  | 1212  |\n| 50 |  0.018 | 4 | 5 | 1603 | 3061  | 1544  |\n| 60 | 0.015 | 5 | 6 | 1973 | 3673  | 1855  |\n| 70 | 0.017 | 5 | 6 | 2512 | 4129  | 2130  |\n| 80 | 0.016 | 6 | 7 | 2845 | 4785  | 2421  |\n\nPlease let us know if any further clarification or experiments would strengthen the paper. We are happy to incorporate additional feedback into the next version of our paper.\n\n---\n---\n> How does this work relate to Bourne et al., 2018? Could you clarify how your results relate to the duality findings of Caffarelli and McCann [9]?\n\nA: First, we appreciate the reviewer for bringing the work of Bourne et al. to our attention. We will cite their work in the final version. These papers also form connections among $\\alpha$-OPT, $\\lambda$-ROT, and weighted Voronoi diagrams. But we respectfully submit that we give a richer characterization that provides deeper insights into these problems. For example, [9] works with the weighted Voronoi diagram under the capped distance function and argues that the mass transported to $b$ in a $\\lambda$-ROT plan is contained within the weighted Voronoi cell centered at $b$. We note that the Voronoi diagram under the capped distance is counter-intuitive, e.g., the boundary of a Voronoi cell may not be connected, and thus hard to interpret.  Similarly, Bourne et al. prove an optimal $\\alpha$-OPT plan routes the mass of $b$ to a subset of the mass within its weighted Voronoi cell, and only find experimentally that the mass transported to each discrete point is contained within a ball (they do not prove the second claim). In contrast, we work with the weighted Voronoi diagram under the Euclidean (or squared Euclidean) distance, where the Voronoi diagram is more standard, and prove that the $\\alpha$-OPT plan routes $b$ to all mass of its restricted Voronoi cell, i.e. its weighted Voronoi cell intersected with the ball of radius equal to the dual weight. With our method, we additionally prove that every point whose mass is not fully routed in the optimal $\\alpha$-partial transport plan has a maximal dual weight, which is an even stronger guarantee. Thus we provide a more detailed characterization. For the sake of space, we did not go into as much detail on the differences between Section 2 and prior work. We appreciate the question and will make this difference more clear in a final version.\n\nFurthermore, [9] only studies the duality between exact solutions of $\\lambda$-ROT and $\\alpha$-OPT. In contrast, our main contribution  (see Section 3) is to show that an approximate solver for one could be used to compute an approximate solution for the other. This result is nontrivial, as shown by the pathological example in Section 3.\n\nWe finally note that our work provides two new combinatorial algorithms to approximate the $\\lambda$-ROT problem to varying degrees of precision."}, {"review_id": "T77eCGnG2j", "reviewer": "Reviewer_CmYB", "summary": "This paper studies two formulations of semi-discrete optimal transport (OT) problem: $\\alpha$-optimal partial transport and $\\lambda$-robust OT. The authors provide theoretical investigations on the optimal solutions in these formulations, connection between the algorithms for each of the settings and present their own algorithm for each of the formulations.", "questions": "- Could you conduct some experiments justifying your theoretical conclusions? \n- Could you please provide a comment regarding the originality of your approach which follows the ideas of the one suggested in (Agarwal et al., 2024)?", "limitations": "The limitations of the approach are stated only in the NeurIPS paper checklist. I kindly suggest the authors create a separate ‘Limitations’ section in the main text or Appendix of their paper and list the limitations there. One limitation is missing here - the paper establishes results only for the $p$-Wasserstein distance as the OT cost.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths.** The paper provides a lot of theoretical results on the connection and properties of partial and robust OT formulations, their solutions and corresponding algorithms.\n\n**Weaknesses.** My main concern corresponds to the **lack of experiments** in the paper. The authors present a lot of theoretical results related to partial/robust OT formulations, propose new algorithms to solve these problems, establish their convergence properties, but **do not test them in any experiments**. I am confused by this fact, since while the paper proposes the algorithms, it should provide at least some experimental validation of their performance. I am also interested in real-world applications of the proposed algorithms, i.e., experiments showing that semi-discrete partial/robust OT formulations and corresponding algorithms are useful in practice. However, at this step, even toy experiments are missing.\n\nThe related concern corresponds to the paper **clarity**. The paper is really hard to follow, most of the text is written in a technical, mathematical style starting from the ‘Introduction’ section. While I understand that it mainly focuses on theoretical investigations, the text should be written in an accessible and understandable form:\n-  I kindly suggest the authors exclude complex explanations of their results from the ‘Introduction’ section (lines 60-94) and rather include the paper contributions in a general form trying to avoid specific mathematical notations and definitions which are then introduced one more time in main sections. \n- I also have concerns regarding the clarity of OT formulations. In the classic OT problem, the marginals of the plan should coincide with input and target measures. However, in lines 98-104, you give the definition of the OT plan replacing the equality constraint on the marginals with the inequalities. While I understand that the purpose was to give a general definition applicable for the partial/robust OT cases, the explanation might be ambiguous for the reader and should be clarified. \n- I am also wondering why the paper lacks a ‘Conclusion’/’Discussion’ section with summarization of the achieved results in a general form. \n- Another minor concern corresponds to the lines 117-118 which operate with the notions of LP formulation/duality which were not given in the text before that.\n\nAnother restriction of the established theory corresponds to the usage of the $p$-Wasserstein distance in the partial/robust OT formulations. This point should be stated as a limitation of this work since general OT formulations might cover many other costs.\n\n**In summary,** I have doubts about the relevance of this paper for the current conference. The paper focuses on the theoretical results related to the  partial/robust OT formulations and does not include proper experimental validation of the suggested algorithms  (no experiments at all). Theoretical results might represent an interest for the ML community, but their motivation, ideas and originality should be clearly explained in the paper’s main text while now the exposition is very technical and hard to follow.  However, currently the exposition of the paper seems to be more suitable for a mathematical venue where the significance of the established theory (nearly 30 pages in Appendix) could be carefully checked.", "was_revised": true, "final_justification": "During the rebuttal, the authors have addressed my main concerns. That is, they conducted the experiments which I have suggested, and promised to improve the clarity of the paper, e.g., improve the paper text, add missing sections, etc. Thus, I increase my score to slightly positive (4). I do not assign higher score since *all* of the experiments were conducted during the rebuttal. I think that the question if this represent a major revision or not could be discussed among the Reviewers and Area Chairs.", "rebuttal": "We appreciate the reviewer’s comment regarding the density of technical details in the introduction. In view of this comment, we will fix the introduction to present the motivation, high-level ideas, and contributions in a clearer and more accessible manner, deferring technical details to later sections. We will also do our best to make the later sections more accessible.\n\n\n---\n---\n> Could you conduct some experiments justifying your theoretical conclusions?\n\n**Summary of Contributions.**\n Our paper presents the following key contributions:\n - A characterization of partial and robust variants of the optimal transport (OT) problem in the semi-discrete setting,\n - Theoretical guarantees on the approximation error when using an approximate $\\lambda$-ROT solver to compute an approximate $\\alpha$-OPT plan, and,\n - Two algorithmic solutions for the $\\lambda$-ROT problem.\n\nWhile the original submission emphasized our theoretical contributions, we appreciate your interest in practical implementation and experimental evaluation. In response, we include below preliminary experimental results. A more extensive evaluation will be added in the next version.\n\n**Implementation Details.**\nWe provide a Python implementation of our algorithm. However, it relies on two core procedures that lack efficient implementations in Python:\n1. Estimating continuous mass within regions defined by the arrangement of the restricted Voronoi cells and their expansions.\n2. Maintaining a dynamic tree data structure as per Sleator and Tarjan (reference [44] in the appendix).\n\nTo address these challenges in our prototype:\n- For (1), we approximate the mass within each region using a fine grid, summing the masses of grid points lying inside the region.\n- For (2), we simulate the dynamic tree by naïvely iterating over all edges.\n\nAn efficient implementation will require optimized versions of these components, which remains a technical challenge in Python.\n\n---\n**Preliminary Results on Algorithm Performance.**\nWe evaluated our implementation using a continuous distribution sampled from a Gaussian distribution inside the unit square. The discrete support consists of $n$ samples drawn from the same distribution, each assigned a mass of $1/n$. We set $\\lambda = 0.2$ and $\\varepsilon = 0.02$. Results are averaged over 10 runs for each $n$.\n\n**Key Findings:**\n- The additive error stays within the target $\\varepsilon = 0.02$.\n- The number of iterations in the two steps of our algorithm (Step 1 and 2) remains well below our upper bounds proved in the paper ($6n$ and $12n$ as described in Lemma 4.3 and D.9).\n- The total path and cycle lengths computed in the two steps grow subquadratically in $n$, aligning with our $O(n^2)$ complexity analysis (Section 4.3).\n\nThe details of our experiments are presented in the following table.\n\n| n | Error | Step 1 iterations | Step 2 iterations | Step 1 paths and cycles length | Step 2 paths and cycles length | Regions |\n| :--: | :--: | :--: | :--: | :--: | :--: | :--: |\n| 10 | 0.009 | 1 | 4 | 41 | 361 | 198  |\n| 20 | 0.018 | 2 | 4 | 320 | 991  | 529  |\n| 30 | 0.021 | 3 | 5 | 688 | 1719  | 862  |\n| 40 | 0.015 | 3 | 5 | 979 | 2385  | 1212  |\n| 50 |  0.018 | 4 | 5 | 1603 | 3061  | 1544  |\n| 60 | 0.015 | 5 | 6 | 1973 | 3673  | 1855  |\n| 70 | 0.017 | 5 | 6 | 2512 | 4129  | 2130  |\n| 80 | 0.016 | 6 | 7 | 2845 | 4785  | 2421  |\n\n\n---\n**Reduction from $\\alpha$-OPT to $\\lambda$-ROT.**\nWe evaluated our reduction from $\\alpha$-OPT to $\\lambda$-ROT by using the Sinkhorn algorithm to approximate $\\lambda$-ROT plans, followed by binary search over $\\lambda$ to compute an $\\alpha$-OPT plan. The continuous distribution in our experiments is an exponential distribution over the unit square, centered at the bottom-left corner, with 10% noise added to the top-right corner. The discrete set $B$ consists of $n=400$ samples from the same exponential distribution, each assigned mass $1/n$. We target a transported mass of $\\alpha$ = 0.9.\n\n*Additive Error:* For any fixed regularization parameter, we evaluate the additive error (with respect to the optimal $\\lambda$-ROT solution) incurred by Sinkhorn, and compare it to the additive error of our computed $\\alpha$-OPT solution (with respect to the optimal $\\alpha$-OPT). We observe that the error in $\\alpha$-OPT closely mirrors that of $\\lambda$-ROT across different choices of regularization parameter. Moreover, as the regularization parameter decreases, the additive errors in $\\alpha$-OPT also diminish. This behavior is consistent with the theoretical guarantee established in Theorem 3.5.\n\n|Regularization | Robust Cost | Partial Cost | Additive error of ROT | Additive error of OPT |\n| :--: | :--: | :--: | :--: | :--: |\n|0.01|0.176|0.084|0.000|0.000|\n|0.02|0.186|0.095|0.009|0.012|\n|0.03|0.182|0.100|0.012|0.017|\n|0.04|0.193|0.112|0.023|0.029|\n|0.05|0.205|0.125|0.034|0.041|\n|0.06|0.216|0.137|0.046|0.053|\n|0.07|0.227|0.148|0.057|0.065|\n|0.08|0.238|0.159|0.068|0.075|\n\n*Approximation factor:* Similarly, we evaluate the approximation factor for each method using the same setup. Once again, we find that the approximation factor in $\\alpha$-OPT closely mirrors that of $\\lambda$-ROT over a range of regularization parameters. As the regularization parameter decreases, the approximation factor in $\\alpha$-OPT also approaches one, further corroborating the theoretical prediction in Theorem 3.5.\n\n|Regularization | Robust Cost | Partial Cost | Relative error of ROT | Relative error of OPT |\n| :--: | :--: | :--: | :--: | :--: | \n|0.01|0.176|0.084|0.998|1.003|\n|0.02|0.186|0.095|1.050|1.142|\n|0.03|0.182|0.100|1.071|1.200|\n|0.04|0.193|0.112|1.136|1.346|\n|0.05|0.205|0.125|1.203|1.493|\n|0.06|0.216|0.137|1.269|1.638|\n|0.07|0.227|0.148|1.334|1.776|\n|0.08|0.238|0.159|1.397|1.904|\n\nPlease let us know if any further clarification or experiments would strengthen the paper. We are happy to incorporate additional feedback into the next version of our paper.\n\n---\n---\n> Could you please provide a comment regarding the originality of your approach which follows the ideas of the one suggested in Agarwal et al., 2024?\n\n**Challenges in Extending Cost-Scaling to Partial Transport:**\n Cost-scaling algorithms have proven effective for solving the complete optimal transport (OT) problem. However, extending these techniques to partial optimal transport introduces substantial challenges—even in fully discrete settings (see, e.g., Ramshaw and Tarjan, FOCS 2012 [37])—and these difficulties are further amplified in the semi-discrete setting.\nA central obstacle lies in the behavior of the dual weights: in complete OT, the optimal solution is invariant under translations of the dual weight vector for the discrete set $B$. In contrast, in partial OT, the absolute magnitudes of the dual weights are crucial, as they control the truncation of Voronoi cells and thus directly determine which regions of the continuous distribution are assigned to each discrete point.\n\n**Limitations of Existing Scaling Algorithms.**\n Recent cost-scaling algorithms (e.g., Agarwal et al., SODA 2024; NeurIPS 2024) increase the magnitudes of the dual weights for the discrete set $B$, and implicitly for the continuous domain $A$, over successive scales. However, these methods lack a mechanism to reduce or correct the dual weights once they are overestimated. This renders existing scaling algorithms unsuitable for partial OT, where accurate dual weight magnitudes are critical.\n\n**Our Solution.**\n To overcome this limitation, we introduce a novel step that uses consolidating paths and augmenting paths to reroute mass in the residual graph and correct overestimated dual weight magnitudes. The addition of this step makes it non-trivial to establish both the correctness of the modified algorithm and its efficiency matching that of Agarwal et al., NeurIPS 2024. These challenges necessitate several new ideas, as detailed in Lemmas 4.2 and 4.3 of the main text and Lemma D.9 in the Appendix.\n\n**Geometric Interpretation of this limitation and consolidating paths.**\n At the start of a scale $\\delta$, some mass from the continuous distribution may lie well inside the restricted Voronoi cells of the discrete points $B$, yet remain untransported—these are what we call violating deficit regions (see Figure 3). In the complete OT setting, such violations are naturally resolved as the plan evolves to transport all mass. In the partial OT setting, however, this issue is non-trivial and must be explicitly resolved to ensure correctness. Our consolidating paths are specifically designed to eliminate all violating deficit regions by restructuring the transport plan to ensure that the total transported mass lies strictly within the restricted Voronoi regions, thereby guaranteeing correctness of the partial OT plan."}, {"review_id": "6iekPP5XT0", "reviewer": "Reviewer_NHUU", "summary": "In this paper, the authors extend the deterministic semi-discrete optimal transport problem, which involves transporting a continuous distribution to a discrete one, to a stochastic setting. Noise and outliers are considered via two variants: $\\alpha$-partial optimal transport, where only an $\\alpha$-fraction of the probability mass is transported, and $\\lambda$-robust optimal transport, where mass can be discarded at a cost of $\\lambda$ per unit.\n\nThe authors propose a novel theoretical characterization of exact solutions using restricted Voronoi cells and weight functions on the discrete distribution domain, under specific conditions. For inexact solutions, they rigorously derive a connection between the two problem variants and introduce a reduction scheme along with corresponding error bounds. To solve the $\\lambda$-robust OT problem with high accuracy, the authors generalize the combinatorial cost-scaling paradigm, incorporating all necessary modifications. The resulting algorithm achieves a complexity bound of $n^{O(d)} \\log(\\Delta/\\varepsilon)$. All results hold for arbitrary cost functions, with additional acceleration techniques developed specifically for Wasserstein distances.", "questions": "How long does it take to compute the $\\Phi$ factor, both for known distributions (e.g., normal, exponential) and unknown distributions? What methods are available for this?\n\n\n Can your methods handle distributions for which only sample data is available?\n\n\nAre the obtained convergence rates optimal with respect to their dependencies on parameters such as n (sample size), d (dimensionality) and others?", "limitations": "The paper lacks practical validation of the proposed methods, particularly comparisons with deterministic counterparts in scenarios where robustness is essential. Even basic 2D experiments on synthetically corrupted data could have provided a proof of concept while demonstrating the methods' robustness and quality.\n\nWhile the theoretical contributions are significant, their novelty appears somewhat incremental, as they primarily combine (in a non-trivial way) existing ideas, problem formulations, and methodologies.\n\nMinor weakness: the notation $\\delta$ is ambiguous, it is used as a parameter in Section 3, and as units in Introduction.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Optimal Transport is an important topic, and new robust algorithms represent valuable contributions to both theory and practice.\n\nThe strongest side of the paper is its complex rigorous theoretical analysis of the proposed methods and properties of the solutions. These findings are interesting and cover all aspects of the considered setups and connections between them.\n\nDespite the amount of information, the paper’s writing is clear and well-structured making it easy to follow.\n\nFor weakness, see the Limitations section.", "was_revised": true, "final_justification": "I'm satisfied authors' answers.", "rebuttal": "We appreciate your thoughtful review and comments. We answer your main questions and concerns below.\n\n\n> How long does it take to compute the $\\Phi$ factor, both for known distributions (e.g. normal, exponential) and unknown distributions? \n\nThe factor of $\\Phi$ depends on how the distribution is represented. If the distribution is a parametric distribution with a constant number of parameters, e.g. normal or exponential (mixture), then a constant degree polygonal query region can be integrated over in constant time using numerical methods. On the other hand if the continuous distribution is represented as a histogram or piecewise-linear function, then the time depends on the complexity of the distribution inside the query region. In low dimensions, one could use some geometric data structures to expedite this computation, but in high dimensions, a simple method would be more practical. Finally, if the distribution is represented as a large representative sample, then the $\\Phi$ factor simply corresponds to assigning each sample point to the region of the partition it is contained in. In turn, our algorithm successfully computes a transport plan between a discrete distribution and a large sample of an unknown continuous distribution.\n\n---\n---\n> Can your methods handle distributions for which only sample data is available?\n\nOur algorithm can handle sample data from a distribution, as long as one global sample is used to approximate the continuous mass inside each region of the arrangement (rather than sampling from the distribution for each mass query).\n\n---\n---\n> Are the obtained convergence rates optimal with respect to their dependencies on parameters such as $n, d$?\n\nA: As noted in Section 2, as well as references [9] and Bourne et al., the optimal semi-discrete transport plan requires computing a (restricted) weighted Voronoi diagram. There are known lower bounds of $n^{\\Omega(d)}$ for the complexity of Voronoi diagrams in the worst case in $d$ dimensions (see e.g. Har-Peled \"Geometric Approximation Algorithms\"), and therefore our high-precision algorithm is near-optimal. While we do not have a rigorous lower bound for the optimal runtime of a low-precision algorithm in constant dimension, we suspect that optimal transport is at least as hard as nearest neighbor search (a naive plan would assign each point to its nearest neighbor). The runtime of our low-precision algorithm is near linear in $n$ and incurs comparable dependence in $\\varepsilon$ and $d$ to the best known algorithms for approximate nearest neighbor search.\n\n---\n---\n\n>The paper lacks practical validation of the proposed methods, particularly comparisons with deterministic counterparts in scenarios where robustness is essential. Even basic 2D experiments on synthetically corrupted data could have provided a proof of concept while demonstrating the methods' robustness and quality.\n\nWhile the original submission emphasized our theoretical contributions, we appreciate your interest in practical implementation and experimental evaluation. In response, we include below preliminary experimental results. A more extensive evaluation will be added in the next version.\n\n**Implementation Details.**\nWe provide a Python implementation of our algorithm. However, it relies on two core procedures that lack efficient implementations in Python:\n1. Estimating continuous mass within regions defined by the arrangement of the restricted Voronoi cells and their expansions.\n2. Maintaining a dynamic tree data structure as per Sleator and Tarjan (reference [44] in the appendix).\n\nTo address these challenges in our prototype:\n- For (1), we approximate the mass within each region using a fine grid, summing the masses of grid points lying inside the region.\n- For (2), we simulate the dynamic tree by naïvely iterating over all edges.\n\nAn efficient implementation will require optimized versions of these components, which remains a technical challenge in Python.\n\n---\n**Preliminary Results on Algorithm Performance.**\nWe evaluated our implementation using a continuous distribution sampled from a Gaussian distribution inside the unit square. The discrete support consists of $n$ samples drawn from the same distribution, each assigned a mass of $1/n$. We set $\\lambda = 0.2$ and $\\varepsilon = 0.02$. Results are averaged over 10 runs for each $n$.\n\n**Key Findings:**\n- The additive error stays within the target $\\varepsilon = 0.02$.\n- The number of iterations in the two steps of our algorithm (Step 1 and 2) remains well below our upper bounds proved in the paper ($6n$ and $12n$ as described in Lemma 4.3 and D.9).\n- The total path and cycle lengths computed in the two steps grow subquadratically in $n$, aligning with our $O(n^2)$ complexity analysis (Section 4.3).\n\nThe details of our experiments are presented in the following table.\n\n| n | Error | Step 1 iterations | Step 2 iterations | Step 1 paths and cycles length | Step 2 paths and cycles length | Regions |\n| :--: | :--: | :--: | :--: | :--: | :--: | :--: |\n| 10 | 0.009 | 1 | 4 | 41 | 361 | 198  |\n| 20 | 0.018 | 2 | 4 | 320 | 991  | 529  |\n| 30 | 0.021 | 3 | 5 | 688 | 1719  | 862  |\n| 40 | 0.015 | 3 | 5 | 979 | 2385  | 1212  |\n| 50 |  0.018 | 4 | 5 | 1603 | 3061  | 1544  |\n| 60 | 0.015 | 5 | 6 | 1973 | 3673  | 1855  |\n| 70 | 0.017 | 5 | 6 | 2512 | 4129  | 2130  |\n| 80 | 0.016 | 6 | 7 | 2845 | 4785  | 2421  |\n\n---\n**Robustness of $\\alpha$-OPT.**\n\nWe examined the effect of 10% noise in the continuous distribution (added in the top-right corner) and measured how much noise was transported under various $\\alpha$ values. Results confirm the robustness of $\\alpha$-OPT: for $\\alpha$ < 0.9, almost no noise mass is transported.\n\n\n|Alpha| Cost| \\% inliers transported| \\% noise transported|\n| :--: | :--: | :--: | :--: | \n |0.700| 0.010| 0.777| 0.000|\n |0.725| 0.011| 0.804| 0.000|\n |0.750| 0.012| 0.830| 0.015|\n |0.775| 0.013| 0.859| 0.005|\n |0.800| 0.014| 0.887| 0.005|\n |0.825| 0.016| 0.914| 0.009|\n |0.850| 0.017| 0.942| 0.014|\n |0.875| 0.019| 0.970| 0.005|\n |0.900| 0.022| 0.996| 0.020|\n |0.925| 0.031| 0.999| 0.244|\n |0.950| 0.045| 0.999| 0.496|\n |0.975| 0.060| 0.999| 0.748|\n\n**Robustness of $\\lambda$-ROT.**\nWe conducted similar experiments for $\\lambda$-ROT. As $\\lambda$ increases, total mass transported grows, but noise transport remains limited when the transported mass is under 0.9.\n\n|Lambda| Cost| Cost ROT| \\% inliers transported| \\% noise transported| Mass transported|\n   | :--: | :--: | :--: | :--: | :--: | :--: | \n| 0.05| 0.014| 0.024| 0.896| 0.004| 0.807|\n| 0.10| 0.020| 0.031| 0.980| 0.009| 0.884|\n   | 0.15| 0.021| 0.037| 0.991| 0.012| 0.894|\n| 0.20| 0.023| 0.042| 0.996| 0.046| 0.902|\n| 0.25| 0.024| 0.047| 0.999| 0.073| 0.907|\n| 0.30| 0.023| 0.051| 0.999| 0.064| 0.907|\n| 0.35| 0.025| 0.055| 0.999| 0.141| 0.914|\n| 0.40| 0.027| 0.061| 1.0| 0.147| 0.915|\n| 0.45| 0.028| 0.063| 1.0| 0.212| 0.921|\n| 0.50| 0.034| 0.066| 1.0| 0.349| 0.935|\n| 0.55| 0.041| 0.068| 1.0| 0.506| 0.951|\n\nPlease let us know if any further clarification or experiments would strengthen the paper. We are happy to incorporate additional feedback into the next version of our paper."}]

Author rebuttal:

Meta-review: The paper investigates the partial and robust variants of semi-discrete OT. All the reviewers agree that this is an important problem and the paper offers novel and interesting insights to this problem via the restricted Laguerre diagrams. The original submission did not contain any numerical verification results of the algorithm, however during the discussion period, such numerical results were provided to the satisfactory of the reviewers. For this reason, I think the paper has been improved and does meet the bar of acceptance. I strongly recommend that the authors incorporate these important information (experimental result, conclusion section) in the final paper if it is accepted.
Decision: Accept (poster)