Preconditioning for the Geometric Transportation Problem

TL;DR

This paper introduces a nearly-linear time approximation scheme for the geometric transportation problem in fixed dimensions, improving efficiency by combining preconditioning techniques with geometric reductions.

Contribution

It presents the first approximation algorithm with nearly-linear time complexity in the number of points, leveraging generalized preconditioning and geometric insights.

Findings

Achieves a (1+ε)-approximate transportation map in nearly-linear time for fixed dimensions.

Reduces the problem to a minimum cost flow on a sparse graph using geometric arguments.

Develops a new preconditioner for the flow problem, enhancing computational efficiency.

Abstract

In the geometric transportation problem, we are given a collection of points $P$ in $d$-dimensional Euclidean space, and each point is given a supply of $\mu(p)$ units of mass, where $\mu(p)$ could be a positive or a negative integer, and the total sum of the supplies is $0$. The goal is to find a flow (called a transportation map) that transports $\mu(p)$ units from any point $p$ with $\mu(p) > 0$, and transports $-\mu(p)$ units into any point $p$ with $\mu(p) < 0$. Moreover, the flow should minimize the total distance traveled by the transported mass. The optimal value is known as the transportation cost, or the Earth Mover's Distance (from the points with positive supply to those with negative supply). This problem has been widely studied in many fields of computer science: from theoretical work in computational geometry, to applications in computer vision, graphics, and machine learning. In this work we study approximation algorithms for the geometric transportation problem. We give an algorithm which, for any fixed dimension $d$, finds a $(1+\varepsilon)$-approximate transportation map in time nearly-linear in $n$, and polynomial in $\varepsilon^{-1}$ and in the logarithm of the total supply. This is the first approximation scheme for the problem whose running time depends on $n$ as $n\cdot \mathrm{polylog}(n)$. Our techniques combine the generalized preconditioning framework of Sherman, which is grounded in continuous optimization, with simple geometric arguments to first reduce the problem to a minimum cost flow problem on a sparse graph, and then to design a good preconditioner for this latter problem.

Equations

$$\sum _{p\in P_{+},q\in P_{-}}{f}_{pq}\Vert p-q{\Vert }_2$$

$$\forall p\in P_{+}:$$

$$\forall q\in P_{-}:$$

$$\forall p\in P_{+},q\in P_{-}:$$

$$E_2=\{(u,v):\ell \in \{0,1,\dots ,L\},C\in {\mathcal{C}}_{\ell },u,v\in N_{\ell }\cap C\text{s.t.}u\ne v\}.$$

$$E_3=\{(N_{\ell -1}(u),u):\ell \in \{1,2,\dots ,L\},u\in N_{\ell }\}.$$

$${\mathrm{dist}}_G(u,v)=\min \sum _{i=1}^k\Vert u_i-u_{i-1}{\Vert }_2,$$

$$\begin{gathered} \{(N_{\ell }(p),N_{\ell -1}(p):\ell (p,q)+1\le \ell \le L\}\cup \{(N_{\ell (p,q)}(p),N_{\ell (p,q)}(q))\} \\ \cup \{(N_{\ell }(q),N_{\ell -1}(q):\ell (p,q)+1\le \ell \le L\}\cup \{(N_L(p),p),(N_L(q),q)\}. \end{gathered}$$

$$\mathbb{P}(C_{\ell }(p)\ne C_{\ell }(q))\le \frac{\Vert p-q{\Vert }_1}{{\Delta }_{\ell }}\le \frac{\sqrt{d}\Vert p-q{\Vert }_2}{{\Delta }_{\ell }}.$$

$$\sum _{(u,v)\in E}{f}_{uv}\Vert u-v{\Vert }_2$$

$$\forall u,v\in V:$$

$$\forall p\in P:$$

$$\forall u\in V\setminus P:$$

$$\mathbb{E}[\text{textsc}cost(G,\mu )]\le (1+O({\epsilon }_{0}L))\text{textsc}cost(P,\mu ).$$

$$\sum _{i=1}^k{a}_i\Vert p_i-q_i{\Vert }_2\le \sum _{i=1}^k{a}_i{\mathrm{dist}}_G(p_i,q_i),$$

$$\mathbb{E}[\text{textsc}cost(G,\mu )]$$

$$\le (1+O({\epsilon }_{0}L))\sum _{i=1}^k{f}_{p_i,q_i}^{\prime }\Vert p_i-q_i{\Vert }_2$$

$$=(1+O({\epsilon }_{0}L))\text{textsc}cost(P,\mu ),$$

$$\Vert f\Vert$$

$$Af=b,.$$

$$\Vert T\Vert =\underset{x\ne 0}{\sup }\frac{\Vert T(x){\Vert }_Y}{\Vert x{\Vert }_X}.$$

$$\kappa (T)=\underset{S}{\inf }\Vert T\Vert \Vert S\Vert ,$$

$$\Vert Af-b{\Vert }_Y\le \beta \Vert A{\Vert }_{\Vert \cdot \Vert \to \Vert \cdot {\Vert }_Y}\Vert f^{\ast }\Vert .$$

$$O\left({\kappa }^2(m+n+M)\log (m)({\epsilon }^{-2}+\log (1/\beta ))\right).$$

$$\Vert g{\Vert }_1$$

$$A{C}^{-1}g=b,$$

$$(Af{)}_u=\sum _{v:(u,v)\in \vec{E}}{f}_{uv}-\sum _{v:(v,u)\in \vec{E}}{f}_{vu}=\sum _{v:(u,v)\in E}{f}_{uv}.$$

$$\Vert h{\Vert }_Y\le \min \{\Vert f{\Vert }_X:f\in {\mathbb{R}}^n,Hf=h\}\le \gamma \Vert h{\Vert }_Y.$$

$$\Vert Hf{\Vert }_Y$$

$$\Vert S(h){\Vert }_X$$

$$\Vert B\tilde{b}{\Vert }_1\le \min \{\Vert f{\Vert }_c:f\in {\mathbb{R}}^{\vec{E}},Af=\tilde{b}\}\le \gamma \Vert B\tilde{b}{\Vert }_1,$$

$$\Vert B\tilde{b}{\Vert }_1=\sum _{p\in P}|{\tilde{b}}_p|\Vert p-N_L(p){\Vert }_2+\frac{1}{L+1}\sum _{\ell =0}^L\frac{{\epsilon }_0{\Delta }_{\ell }}{4}\sum _{\tilde{C}\in {\tilde{\mathcal{C}}}_{\ell }}\left|\sum _{u\in \tilde{C}\cap V}{\tilde{b}}_u\right|.$$

$$\text{textsc}cost(G,\tilde{b})=\min \{\Vert f{\Vert }_c:f\in {\mathbb{R}}^{\vec{E}},Af=\tilde{b}\}.$$

$$\frac{{\epsilon }_0{\Delta }_{\ell }}{4}\sum _{\tilde{C}\in {\tilde{\mathcal{C}}}_{\ell }}\left|\sum _{u\in \tilde{C}\cap V}{\tilde{b}}_u\right|\le \text{textsc}cost(G,\tilde{b})-\sum _{p\in P}|{\tilde{b}}_p|\Vert p-N_L(p){\Vert }_2.$$

$$\delta (u,f)=\sum _{v:(u,v)\in \vec{E}}{f}_{uv}-\sum _{v:(v,u)\in \vec{E}}{f}_{vu}-{\tilde{b}}_u=\sum _{v:(u,v)\in E}{f}_{uv}-{\tilde{b}}_u.$$

$$\delta (u,f)=-\sum _{v\in \tilde{C}\cap V}{\tilde{b}}_v=-\sum _{v:N_{\ell -1}(v)=u}{\tilde{b}}_v.$$

Full text

In this work we study approximation algorithms for the geometric transportation problem. We give an algorithm which, for any fixed dimension $d$, finds a $(1+\varepsilon)$-approximate transportation map in time nearly-linear in $n$, and polynomial in $\varepsilon^{-1}$ and in the logarithm of the total supply. This is the first approximation scheme for the problem whose running time depends on $n$ as $n\cdot\mathrm{polylog}(n)$. Our techniques combine the generalized preconditioning framework of Sherman, which is grounded in continuous optimization, with simple geometric arguments to first reduce the problem to a minimum cost flow problem on a sparse graph, and then to design a good preconditioner for this latter problem.

1 Introduction

In the Geometric Transportation problem, we are given a set $P$ of $n$ points in $d$-dimensional Euclidean space, and a function $\mu:P\to\mathbb{Z}$ that assigns a (positive or negative integer) weight to each point, so that $\sum_{p\in P}{\mu(p)}=0$. We call $\mu$ the supply function. We can think of each point $p\in P$ as either a pile of earth, or a hole, and $\mu(p)$ gives, respectively, the amount of earth (if positive) or the size of the hole (if negative). The constraint on $\mu$ means that the total space in the holes equals the total amount of earth in the piles. Our goal is to find the most efficient way to transport the earth to the holes, where the cost of transporting a unit of mass from $p$ to $q$ equals the distance $\|p-q\|_{2}$ it must travel, measured in the Euclidean norm. More formally, let $P_{+}=\{p\in P:\mu(p)\geq 0\}$ be the “piles” and $P_{-}=\{q\in P:\mu(q)<0\}$ be the “holes”. We want to solve the minimum cost flow problem given by the following linear program

[TABLE]

Above, the constraints enforce that all demands are satisfied, i.e. all the earth goes into holes, and no hole is overfilled. A vector $f\in\mathbb{R}^{P_{+}\times P_{-}}$ that satisfies these constraints is called a transportation map for $(P,\mu)$. The cost of an optimal transportation map is called the transportation cost of $(P,\mu)$, or also the Earth Mover’s Distance, or $1$-Wasserstein distance between the measures $\phi=\mu|_{P_{+}}$ and $\psi=-\mu|_{P_{-}}$, and is denoted $\textsc{cost}(P,\mu)$.

The geometric transportation problem is a discretized version of the continuous optimal mass transportation problem of Monge and Kantorovich. This is a classical problem in mathematics, with many beautiful connections to partial differential equations, geometry, and probability theory, among others [Vil03]. The discrete problem above has also found a large number of applications, for example to shape matching [GV02, GD04] and image retrieval [RTG00, LCL04] in computer vision, and interpolation between distributions in graphics [BvdPPH11]. It is also of intrinsic interest as a natural problem in computational geometry.

The geometric transportation problem can be modeled as a minimum cost flow problem in the complete bipartite graph with bi-partition $P_{+}\cup P_{-}$. This graph is dense, and solving the minimum cost flow problem exactly using the algorithm of Lee and Sidford would take time $O(n^{2.5}\ \mathrm{polylog}(U))$,111We allow the implicit constants in the asymptotic notation to depend on the dimension $d$, except when the asymptotic notation is used in an exponent. where $U=\sum_{p\in P_{+}}\mu(p)$. If instead we settle for an $(1+\varepsilon)$-approximation, then the recent algorithm of Sherman [She17] gives running time $O(n^{2+o(1)}\varepsilon^{-2})$. However, there exist faster (in certain parameter regimes) approximation algorithms which exploit the geometric structure of the problem: after a long line of work, Sharathkumar and Agarwal [SA12a] gave an algorithm that computes a $(1+\varepsilon)$-approximation in time $O(n\sqrt{U}\ \mathrm{polylog}(U,\varepsilon,n))$, and, recently, Agarwal et al. [AFP*+*17] gave several algorithms with different trade-offs, among them an $(1+\varepsilon)$-approximation algorithm running in time $O(n^{3/2}\varepsilon^{-d}\ \mathrm{polylog}(U,n))$. A nearly linear time $(1+\varepsilon)$-approximation algorithm was given by Sharathkumar and Agarwal [SA12b] in the special case of unit supplies; their algorithm runs in time $O(n\varepsilon^{-O(d)}\ \mathrm{polylog}(n))$. Until now, no such algorithm was known for general supply functions.

A related line of work focuses on estimating the transportation cost, without necessarily computing a transportation map. An influential paper in this direction was that of Indyk [Ind07], whose algorithm gives a constant factor approximation to the transportation cost in the case of unit supplies, in time $O(n\ \mathrm{polylog}(n))$. This result was extended to arbitrary supplies and to approximation factor $1+\varepsilon$ by Andoni et al. [ANOY14], whose algorithm runs in time $O(n^{1+o(1)})$, with the $o(1)$ factor in the exponent hiding dependence on $\varepsilon$. The result of Sherman [She17], mentioned above, together with the existence of sparse Euclidean spanners, implies an $O(n^{1+o(1)}\varepsilon^{-O(d)})$ time algorithm to estimate the transportation cost.222We are indebted to an anonymous reviewer for this observation. Indeed, using the well-separated pair decompositions of Callahan and Kosaraju [CK93] (see also Chapter 3 of [HP11]), one can construct in time $O(n\log n+n\varepsilon^{-d})$ a graph $G$ on $P$ with $O(n\varepsilon^{-d})$ edges, so that the shortest path in $G$ between any $p,q\in P$ has total length at most $(1+\varepsilon)\|p-q\|_{2}$. It is then simple to construct a minimum cost flow problem on $G$ whose optimal value equals the transportation cost of $(P,\mu)$ up to a factor of $1+\varepsilon$. This minimum cost flow problem in $G$ can then be solved with Sherman’s algorithm to achieve the claimed running time. It is not immediately clear, however, how to use the flow in $G$ to construct a transportation map of comparable cost in nearly linear time.

There is a related line of work [PA06, AS14, ABRW18] that studies the transportation problem when the cost of transporting mass from $p$ to $q$ is $\|p-q\|_{2}^{r}$, giving the $r$-Wasserstein distance. This appears to be a more challenging problem, and we do not address it further.

1.1 Our Results and Methods

Let us recall that the aspect ratio (or spread) of a pointset $P$ is defined as the ratio of its diameter to the smallest distance between two distinct points. Our main result is a nearly linear time $(1+\varepsilon)$-approximation algorithm for the geometric transportation problem, as captured by the following theorem.

Theorem 1**.**

There exists a randomized algorithm that on input an $n$-point set $P\subset\mathbb{Q}^{d}$ with aspect ratio $\Delta$, and supply function $\mu:P\to\mathbb{Z}$, runs in time $O(n\varepsilon^{-O(d)}\log(\Delta)^{O(d)}\log n)$, and with probability at least $1/2$ finds a transportation map with cost at most $(1+\varepsilon)~{}\cdot~{}\textsc{cost}(P,\mu)$.

There also exists a randomized algorithm that on input an $n$-point set $P\subset\mathbb{Q}^{d}$ and supply function $\mu:P\to\mathbb{Z}$ such that $U=\sum_{p\in P_{+}}{\mu(p)}$, runs in time $O(n\varepsilon^{-O(d)}\log(U)^{O(d)}\log(n)^{2})$, and with probability at least $1/2$ finds a transportation map with cost at most $(1+\varepsilon)~{}\cdot~{}\textsc{cost}(P,\mu)$.

In constant dimension, the dependence of the running time on the aspect ratio $\Delta$ or total supply $U$ is polylogarithmic, and the dependence on the approximation $\varepsilon$ is polynomial. The dependence on $n$ is just $O(n\log n)$ (respectively $O(n\log(n)^{2})$). This is in contrast with prior work which either had a much larger dependence on $n$, or a polynomial dependence on $U$.

In the proof of this result, we employ a combination of geometric tools, and tools from continuous optimization developed for the general minimum cost flow problem. As a first step, we reduce the transportation problem to an (uncapacitated) minimum cost flow problem on a random sparse graph. This construction is closely related to the prior work of Sharathkumar and Agarwal in the case of unit capacities [SA12b], and also to the estimation algorithm for the transportation cost in [ANOY14]. In particular, the sparse graph and minimum cost flow instance we construct are a simplification of the minimum cost flow instance in [ANOY14]. Together with the recent work of Sherman [She17], this reduction is enough to get a $O(n^{1+o(1)}\varepsilon^{-O(d)})$ time algorithm to estimate the transportation cost. As mentioned above, this running time can also be achieved using a Euclidean spanner, and the random sparse graph we use can, in fact, be seen as a randomized spanner. Our graph, however, has a nice hierarchical structure not shared by other spanner constructions. In particular, there is a quadtree of the input point set such that any edge leaving a cell of the quadtree goes to either a sibling cell or a parent cell. This property is useful both for improving the running time further, and for computing a transportation map.

To further improve the running time, we open up Sherman’s elegant framework, and combine it with classical geometric constructions. Sherman’s framework is based on repeatedly finding a flow which is approximately optimal, and approximately feasible. This can be achieved quickly, for example using the multiplicative weights update method [AHK12]. Once an approximately feasible flow is found, the supplies and demands are updated to equal the still-unrouted part of the original supplies and demands, and the problem is solved again. Finally, all the flows computed so far are added up; if the resulting flow is very close to being feasible, it can be “rounded off” to a truly feasible one with a small loss in the cost. Sherman showed that, if the problem being solved is sufficiently well conditioned in the appropriate way, then this process of repeatedly finding approximately feasible flows converges very fast to a feasible one.

Unfortunately, most minimum cost flow problems are not well-conditioned in their natural formulations. For this reason, we need to find a preconditioner, i.e. a linear transformation, to be applied to the constraints of the minimum cost flow problem, which makes the problem well-conditioned. This is where we can further exploit the structure of the geometric transportation problem. Constructing a preconditioner in our random sparse graph is closely related to estimating the transportation cost by an embedding into $\ell_{1}$. We exploit ideas from the known embeddings of transportation cost into $\ell_{1}$ [Cha02, IT03] to construct a preconditioner with condition number that depends polynomially on the approximation factor and the logarithm of the aspect ratio. The insight that these simple and well-known techniques can be repurposed to give high-quality preconditioners is one of our main conceptual contributions.

The final algorithm does not require any specialized data structures and is relatively simple: once the sparse graph is constructed, and the preconditioner is applied, we iteratively run the multiplicative weights method several times, and, once we are close enough to a feasible flow, greedily round off the errors. This is sufficient for estimating the transportation cost, but does not yet give us a transportation map, because we transformed the original problem into a different minimum cost flow problem on a sparse graph. In the final section of our paper we show a simple way to efficiently transform a flow on this sparse graph into a transportation map, without increasing its cost.

2 Notation and Basic Notions

We use the notation $\|x\|_{p}$ for the $\ell_{p}$ norm of a vector $x\in\mathbb{R}^{d}$: $\|x\|_{p}=\left(\sum_{i=1}^{d}{|x_{i}|^{p}}\right)^{1/p}$. We will assume that the input $P\subset\mathbb{Q}^{d}$ to our problem lies in $[0,\Delta]^{d}$, and that the smallest Euclidean distance between any two points is at least $1$. Since any shifts, rotations, or dilations of the points do not change the problem, this assumption is equivalent to assuming that the aspect ratio of $P$, i.e. the ratio between the diameter of $P$ and the smallest distance between two distinct points, is bounded between $\Delta$ and $\sqrt{d}\Delta$.

In fact, at the cost of a small increase in the approximation factor, we can reduce to this case of bounded aspect ratio. The reduction produces a point set $P$ for which $\Delta$ is bounded by a factor that depends on $U=\sum_{p\in P_{+}}{\mu(p)}$.

Lemma 1** ([Ind07, ANOY14]).**

Suppose that there exists a function $T$, increasing in all its parameters, and an algorithm which, for any $\varepsilon,\delta<1$, on input of size $n$ in $[0,\Delta]^{d}\cap\mathbb{Z}^{d}$, runs in time $O(nT(\Delta,\varepsilon,\delta))$, and with probability $1-\delta$ computes a $1+\varepsilon$ approximation to the geometric transportation problem. Then there exists an algorithm that takes any input $P\subset\mathbb{R}^{d}$ of size $n$ and $\mu:P\to\mathbb{Z}$ such that $U=\sum_{p\in P_{+}}{\mu(p)}$, runs in time $O(nT(cdU/\varepsilon^{2},\varepsilon,\delta n))$ for an absolute constant $c$, and, with probability $1-\delta$, achieves an approximation of $1+O(\varepsilon)$ for the geometric transportation problem on $P$ and $\mu$.

With this lemma, the second algorithm in Theorem 1 follows from the first one. For this reason, we can focus our presentation on the case of bounded aspect ratio.

For a set $V$, we call $f\in\mathbb{R}^{V\times V}$ a flow if it is anti-symmetric, i.e. $f_{uv}=-f_{vu}$ for every $u,v\in V$. Intuitively, we think of $f_{uv}>0$ as flow going in the direction from $u$ to $v$. For a graph $G=(V,E)$, we define a flow on $G$ to be a flow $f$ supported on $E$, i.e. one such that $f_{uv}=0$ for any $(u,v)\not\in E$. The divergence of a flow $f$ at $u$ is the quantity $\sum_{v\in V}{f_{uv}}$, i.e. the excess of the flow leaving $u$ over the flow entering $u$. When the divergence at $u$ is [math], we say that flow conservation is satisfied at $u$.

In general, we assume that our graphs are stored in the adjacency list representation, and that flow values are stored with the adjacency list.

3 Reduction to Minimum Cost Flow on a Sparse Graph

The first step in our algorithm is to reduce the geometric transportation problem, which is naturally modeled as a minimum cost flow problem on a complete bipartite graph, to another minimum cost flow on a sparse random graph. The following construction is a simplified version of one used in [ANOY14].

3.1 Graph Construction

Recall that $P\subset[0,\Delta^{d}]$. We start by constructing a grid containing all the points in $P$ as follows. We sample a uniformly random point $x\in[0,\Delta]^{d}$, and define the cell $C_{0}=[-\Delta,\Delta]^{d}+x$. Note that all points in $P$ are in $C_{0}$. We say that $C_{0}$ is on level 0. The set of cells on level $\ell$ is labelled $\mathcal{C}_{\ell}$ and is constructed by taking each cell $C$ in $\mathcal{C}_{\ell-1}$ and dividing it into $2^{d}$ equally-sized cubes of side length half the side length of $C$, and retaining those that contain points in $P$. We say that $C$ is the parent of the cells into which it was divided. Conversely, those $2^{d}$ cells are called the children of $C$. This division is continued until level $L$, where all points in $P$ lie in different cells in $\mathcal{C}_{L}$. With probability $1$, no point of $P$ lands on the boundary of any cell, so we will ignore this event.

Since the side length of any cell on level $\ell$ is half of the side length of a cell on level $\ell-1$, any cell in $\mathcal{C}_{\ell}$ has side length $2^{1-\ell}\Delta$, and we refer to this length as $\Delta_{l}$. Also note that, since the closest two points in $P$ are at least distance $1$ apart, we have $L\leq\log_{2}(2\sqrt{d}\Delta)$, so $L=O(\log\Delta)$. Moreover, any point $p\in P$ lies in at most $L+1$ cells, one per level, and, since in $\mathcal{C}_{\ell}$ we only retain cells that contain points in $P$, we have that $|\mathcal{C}_{0}\cup\ldots\cup\mathcal{C}_{L}|\leq n(L+1)=O(n\log\Delta)$.

Let $\varepsilon_{0}$ be a small number such that $\varepsilon_{0}^{-1}$ is an even integer. This $\varepsilon_{0}$ is not the same as the value used when we speak of a $(1+\varepsilon)$-approximation, although the two are related, as we will see. Next, for each $\ell\in\{0,\ldots,L$, we take each cell $C\in\mathcal{C}_{\ell}$, and we divide it into $\varepsilon_{0}^{-d}$ subcells, each of side length $\varepsilon_{0}\Delta_{l}$. The set of all such subcells of all $C\in\mathcal{C}_{\ell}$ is denoted $\mathcal{\widetilde{C}}_{\ell}$. To each subcell we associate a net point at its centre and we denote the set of all net points on level $\ell$ as $N_{\ell}$. See Figure 1 for an illustration. Note that $|N_{\ell}|={|\mathcal{C}_{\ell}|}{\varepsilon_{0}^{-d}}$, so $|N_{0}\cup\ldots\cup N_{L}|\leq{n(L+1)}{\varepsilon_{0}^{-d}}=O({n}{\varepsilon_{0}^{-d}}\log\Delta)$.

If $p$ is a point in $C_{0}$, let $C_{\ell}(p)\in\mathcal{\widetilde{C}}_{\ell}$ be the unique subcell on level $\ell$ that contains $p$. Similarly, let $N_{\ell}(p)$ be the net point on level $\ell$ closest to $p$. Equivalently, $N_{\ell}(p)$ is the center of $C_{\ell}(p)$. We say that $N_{\ell}(p)$ and $N_{\ell-1}(p)$ are each other’s child and parent, respectively.

We are now ready to construct our graph $G=(V,E)$. We set $V=P\cup N_{0}\cup N_{1}\cup\dots\cup N_{L}$, the set of all of the points in the original problem and all of the netpoints we have constructed. By the previous discussion, $|V|=O({n}{\varepsilon_{0}^{-d}}\log\Delta)$. We build our graph using three types of edges. The first set, $E_{1}$, connects points $p$ in $P$ to their closest net point on level $L$, i.e. $E_{1}=\{(p,N_{L}(p)):p\in P\}$. The size of this set is $|E_{1}|=n$. The second set of edges, $E_{2}$, connects all net points of subcells of a given cell pairwise. Thus,

[TABLE]

Since $|\mathcal{C}_{\ell}|\leq n$, and, for any $C\in\mathcal{C}_{\ell}$, $|N_{\ell}\cap C|=\varepsilon_{0}^{-d}$, we have that $|E_{2}|\leq(L+1)n\varepsilon_{0}^{-2d}=O(n\varepsilon_{0}^{-2d}\log\Delta)$. The last set of edges, $E_{3}$, connects net points to their parent net points, i.e.

[TABLE]

The size of $E_{3}$ is less than the total number of net points, so $|E_{3}|=O({n}{\varepsilon_{0}^{-d}}\log\Delta)$. Then, the set of edges is defined as $E=E_{1}\cup E_{2}\cup E_{3}$, and this completes the description of $G=(V,E)$. The number of edges is dominated by the size of $E_{2}$, and is in $O(n\varepsilon_{0}^{-2d}\log\Delta)$. We will be calculating distances between points along paths in this graph, and we define the distance/cost of any edge to equal the Euclidean distance between the two endpoints of the edge.

For the rest of the paper we will use $G$ to refer to the random graph constructed above. Having stored the random shift of $C_{0}$, we can enumerate the vertices in time $O({n}{\varepsilon_{0}^{-d}}\log\Delta)$ by going through each point in $P$ and checking which cells it lies in. We can then store the vertices in a static dictionary [FKS84]. Similarly, we can construct the adjacency lists of the vertices in time $O(n\varepsilon_{0}^{-2d}\log\Delta)$.

3.2 Embedding Euclidean Distance into a Sparse Graph

Let $\mathrm{dist}_{G}(u,v)$ be the shortest path distance between two nodes $u$ and $v$ of $G$, i.e.

[TABLE]

with the minimum taken over paths $u=u_{0},\ldots,u_{k}=v$ connecting $u$ and $v$ in $G$. Our goal for this section is to show that, for any two $p$ and $q$ in $P$, the expected value of $\mathrm{dist}_{G}(u,v)$ is close to the Euclidean distance between them. In other words, we want to show that our random graph construction provides a randomized embedding of Euclidean distance into the shortest path metric of a sparse graph. This is similar to the embeddings used in [ANOY14] and [SA12a, SA12b], and can be seen as a randomized spanner construction, in which the graph $G$ only needs to have low stretch in expectation.

For two points $p$ and $q$ in $P$, we define $\ell(p,q)$ to be the level of the smallest cell containing both $p$ and $q$, which can depend on the random shift of the grid. This means that $C_{\ell(p,q)}(p)=C_{\ell(p,q)}(q)$, but $C_{\ell(p,q)+1}(p)\neq C_{\ell(p,q)+1}(q)$. The following definition will be useful.

Definition 1**.**

Given two points $p,q\in P$, the canonical path between $p$ and $q$ in the random graph $G$ consists of the edges

[TABLE]

The canonical path could be much longer than $\|p-q\|_{2}$ because, if $p$ and $q$ are two nearby points that are separated by a grid line, the canonical path will be much larger than the Euclidean distance. However, when $p$ and $q$ are very close together, the likelihood that such a grid line will fall between them is low. To formalize this intuition, we need the following well-known bound on the probability that $p$ and $q$ will be split at a certain level of the tree, which is elementary, and goes back at least to Arora’s work on Euclidean TSP [Aro98].

Lemma 2**.**

Over the random shift of $C_{0}$, $\forall p,q\in P$,

[TABLE]

Using this standard lemma, we can prove the bound on the distortion between Euclidean distance and the expected shortest path distance in $G$. The proof, which is standard, is deferred to the appendix.

Lemma 3**.**

Let $p,q\in P$. Then, for any shift of $C_{0}$, $\|p-q\|_{2}\leq\mathrm{dist}_{G}(p,q)$. Moreover, over the random shift of $C_{0}$, $\mathbb{E}[\mathrm{dist}_{G}(p,q)]\leq(1+O(\varepsilon_{0}L))\|p-q\|_{2}$.

3.3 Approximating the Transportation Cost with Minimum Cost

Flow

We are now ready to translate the transportation problem to a minimum cost flow problem in $G$, captured by the following linear program.

[TABLE]

Let us denote by $\textsc{cost}(G,\mu)$ the value of the optimal solution to the problem (1)–(4). Note that, since $G$ is a random graph, $\textsc{cost}(G,\mu)$ is also a random variable.

The next theorem shows that that $\textsc{cost}(G,\mu)$ approximates the transportation cost $\textsc{cost}(P,\mu)$.

Theorem 2**.**

For any shift of $C_{0}$, $\textsc{cost}(P,\mu)\leq\textsc{cost}(G,\mu)$. Moreover, over the random shift of $C_{0}$,

[TABLE]

Proof of Theorem 2.

To prove the first claim, we show that a flow ${f}$ satisfying (3)–(4) can be transformed into a transportation map ${f}^{\prime}$ for $P$ and $\mu$ without increasing its cost. Towards this goal, let us decompose $f$ into a sum of paths as $f=\sum_{i=1}^{k}{a_{i}f_{\pi_{i}}}$, where $\pi_{i}$ is a path from some point $p_{i}$ in $P_{+}$ to some point $q_{i}$ in $P_{-}$, $f_{\pi_{i}}$ is the unit flow traveling on edges of $\pi_{i}$ in the direction from $p_{i}$ to $q_{i}$, and $a_{1},\ldots,a_{k}$ are non-negative integers. The existence of such a decomposition is a standard fact, see e.g. Chapter 10 of [Sch03]. The cost of $f$ is at least $\sum_{i=1}^{k}{a_{i}\mathrm{dist}_{G}(p_{i},q_{i})}$. We can construct $f^{\prime}$ by defining $f^{\prime}_{pq}=\sum_{i:p_{i}=p,q_{i}=q}{a_{i}}$. The constraints on $f$ imply that $f^{\prime}$ is a feasible transportation map; its cost is

[TABLE]

where the inequality follows from Lemma 3, and the right hand side is, as already observed, at most the cost of $f$. The first claim then follows since we can take $f$ to be an optimal solution to (1)–(4).

For the second claim, we need to go in reverse. Let $f^{\prime}$ be an optimal transportation map, i.e. one achieving $\textsc{cost}(P,\mu)$. Let $(p_{1},q_{1}),\ldots,(p_{k},q_{k})$ be the edges in the support of $f^{\prime}$, where, for every $i$, $p_{i}\in P_{+}$ and $q_{i}\in P_{-}$. Let $\pi_{i}$ be a shortest path in $G$ between $p_{i}$ and $q_{i}$. We define a flow $f$ in $G$ by $f=\sum_{i=1}^{k}{f^{\prime}_{p_{i},q_{i}}f_{\pi_{i}}}$, where $f_{\pi_{i}}$ is, as before, a unit flow from $p_{i}$ to $q_{i}$ supported on the edges of $\pi_{i}$. The feasibility of $f$ follows from the fact that $f^{\prime}$ is a feasible transportation map. Moreover, the cost of $f$ is $\sum_{i=1}^{k}{f^{\prime}_{p_{i},q_{i}}\mathrm{dist}_{G}(p_{i},q_{i})}$. We then have

[TABLE]

where the final inequality follows from Lemma 3. ∎

4 Solving the Minimum Cost Flow Problem

In this section we describe the generalized preconditioning framework of Sherman [She17], and describe how to use it to solve the minimum cost flow problem (1)–(4) in nearly linear time.

4.1 Generalized Preconditioning Framework

We take a short detour to describe the generalized preconditioning framework of Sherman. Suppose that $\|\cdot\|$ is a norm on $\mathbb{R}^{m}$, and we are given inputs $A\in\mathbb{R}^{n\times m}$ and $b\in\mathbb{R}^{n}$. The framework is concerned with solving the following optimization problem over the variables $f\in\mathbb{R}^{m}$:

[TABLE]

This formulation can capture many important problems, including, crucially for us, the uncapacitated minimum cost flow problem,333Uncapacitated minimum cost flow is commonly known as the transportation problem. We do not use this terminology, in order to avoid confusion with the geometric transportation problem. as we explain later.

A central definition in the framework is that of the nonlinear condition number. Before we state it, let us recall that the norm of a linear map $T:X\to Y$, where $(X,\|\cdot\|_{X})$ and $(Y,\|\cdot\|_{Y})$ are two finite dimensional normed vector spaces, is defined as

[TABLE]

We adopt the same definition for non-linear maps, as well. For an $n\times m$ matrix $A$, and norms $\|\cdot\|_{X}$ on $X=\mathbb{R}^{n}$ and $\|\cdot\|_{Y}$ on the column span $Y$ of $A$, we use the notation $\|A\|_{\|\cdot\|_{X}\to\|\cdot\|_{Y}}$ to refer to the norm of the linear map from $X$ to $Y$ represented by $A$.

Having recalled these notions, we proceed to define the condition number.

Definition 2**.**

Let $(X,\|\cdot\|_{X})$ and $(Y,\|\cdot\|_{Y})$ be two finite dimensional normed vector spaces. The non-linear condition number $\kappa$ of a linear map $T:X\to Y$ is defined by

[TABLE]

where $S:Y\to X$ ranges over all (not necessarily linear) maps such that for all $x\in X$ we have $T(S(Tx))=Tx$.

For an $n\times m$ matrix $A$, and norms $\|\cdot\|_{X}$ on $X=\mathbb{R}^{m}$ and $\|\cdot\|_{Y}$ on the column span $Y$ of $A$, we define $\kappa_{\|\cdot\|_{X}\to\|\cdot\|_{Y}}(A)$ as $\kappa(T)$, where $T:X\to Y$ is the linear map represented by the matrix $A$.

This definition generalizes the standard condition number, which is the special case in which both norms are taken to be Euclidean.

The generalized preconditioning framework is based on composing rough approximation algorithms to get a high-quality approximation. The rough approximation algorithms are further allowed to violate the constraints of the problem slightly. Thus, they achieve a bi-criteria approximation, captured by the following definition.

Definition 3**.**

Let $\|\cdot\|_{Y}$ be a norm defined on the range of the matrix $A$ in (6). Let $f^{\ast}$ be an optimal solution to the problem (5)–(6). Then $f\in\mathbb{R}^{m}$ is called an $(\alpha,\beta)$-solution (with respect to $\|\cdot\|_{Y}$) to (5)–(6) if $\|f\|\leq\alpha\|f^{\ast}\|$ and, moreover,

[TABLE]

An algorithm that, when given inputs $A$ and $b$, outputs an $(\alpha,\beta)$-solution $f$ is called an $(\alpha,\beta)$-solver (with respect to $\|\cdot\|_{Y}$) for (5)–(6).

We use the next result of Sherman, which gives a solver for (5)–(6) with running time controlled by condition number of the constraint matrix.

Theorem 3** ([She17]).**

Let $\varepsilon,\beta>0$, and suppose that the norm in (5) is the $\ell_{1}^{m}$ norm. Let $\kappa=\kappa_{\|\cdot\|_{1}\to\|\cdot\|_{1}}(A)$, and let $M$ be an upper bound on the time necessary to compute matrix-vector products with $A$ and $A^{\top}$. Then there exists a $(1+\varepsilon,\beta)$-solver with respect to $\|\cdot\|_{1}$ for the problem (5)–(6) with running time bounded by

[TABLE]

This algorithm is based on repeatedly applying an $(\alpha,\beta)$ solver with much worse dependence on $\beta$; in fact, at the cost of a slightly worse dependence on $\kappa$, for the latter one can use a simple solver based on the multiplicative weights update method. To make our paper more self-contained, we sketch this latter algorithm in the appendix.

Rescaling $f$ coordinatewise gives us the following easy corollary.

Corollary 1**.**

Let $\varepsilon,\beta>0$, and let $c\in\mathbb{R}_{>0}^{m}$. Suppose that the norm in (5) is given by $\|f\|_{c}=\sum_{i=1}^{m}{c_{i}f_{i}}$. Let $\kappa=\kappa_{\|\cdot\|_{c}\to\|\cdot\|_{1}}(A)$, and let $M$ be the time necessary to compute matrix-vector products with $A$ and $A^{\top}$. Then there exists a $(1+\varepsilon,\beta)$-solver with respect to $\|\cdot\|_{1}$ for the problem (5)–(6) with running time bounded by $O\left(\kappa^{2}(m+n+M)\log(m)(\varepsilon^{-2}+\log(1/\beta))\right).$

Proof.

Define the diagonal matrix $C$ by $c_{ii}=c_{i}$. The problem (5)–(6) with norm $\|\cdot\|=\|\cdot\|_{c}$ can be reformulated, after the change of variables $g=Cf$, as

[TABLE]

Then the corollary follows from Theorem 3 since $\kappa_{\|\cdot\|_{c}\to\|\cdot\|_{1}}(A)=\kappa_{\|\cdot\|_{1}\to\|\cdot\|_{1}}(AC^{-1})$. The latter follows from two observations: first, $\|AC^{-1}\|_{\|\cdot\|_{1}\to\|\cdot\|_{1}}=\|A\|_{\|\cdot\|_{c}\to\|\cdot\|_{1}}$; second, the map $G$ from the column span $Y$ of $A$ to $\mathbb{R}^{m}$ satisfies $AGAf=Af$ for all $f$ if and only if $AC^{-1}(CG(AC^{-1}g))$ for all $g$, and $\|CG\|_{\|\cdot\|_{1}\to\|\cdot\|_{1}}=\|G\|_{\|\cdot\|_{1}\to\|\cdot\|_{c}}$. ∎

While Theorem 3 and Corollary 1 allow us to find solutions which are very close to being feasible, they do not give an exactly feasible solution. The final result we use from Sherman allows us to use a solver which exactly satisfies the constraints, but only achieves a large approximation ratio, to round an approximately feasible solution to an exactly feasible one.

To state the result, we need to define the composition of two solvers $\mathcal{F}$ and $\mathcal{F}^{\prime}$. The composed algorithm $\mathcal{F}^{\prime}\circ\mathcal{F}$ first calls $\mathcal{F}$ on $A$ and $b$ to get $f$; then it calls $\mathcal{F}^{\prime}$ on $A$ and the residual vector $b-Af$ to get a solution $f^{\prime}$; finally it outputs $f+f^{\prime}$.

Theorem 4** ([She17]).**

Let $C\geq 1$, and $\varepsilon,\beta>0$. Let $\kappa=\kappa_{\|\cdot\|\to\|\cdot\|_{Y}}(A)$, where $\|\cdot\|$ is the norm in (5), and $\|\cdot\|_{Y}$ is some norm on the column span of $A$. Then, if $\mathcal{F}$ is a $(1+\varepsilon,\varepsilon\beta/\kappa)$-solver, and $\mathcal{F}^{\prime}$ is a $(C,0)$-solver, the composition $\mathcal{F}^{\prime}\circ\mathcal{F}$ is a $(1+\varepsilon+C\varepsilon\beta,0)$-solver.

4.2 Preconditioner Construction

The formulation (5)–(6) captures the (uncapacitated) minimum cost flow problem with arbitrary demands. We will explain this for the graph $G=(V,E)$ defined in Section 3.2 and the minimum cost flow problem (1)–(4). Let us pick an arbitrary orientation on the edges $E$, and call the directed edges $\vec{E}$. Then take $A$ to be the directed vertex by edge incidence matrix of $G$: it is indexed by $V\times\vec{E}$, and for any node $u$, and any directed edge $e=(v,w)$, set $A_{u,e}=1$ if $u=v$, $A_{u,e}=-1$ if $u=w$, and $A_{u,e}=0$ otherwise. We represent a flow $f$ by its restriction to $\vec{E}$, seen as a vector in $\mathbb{R}^{\vec{E}}$. Slightly abusing notation, we will use the letter $f$ both for this vector, and for the flow, which is defined both for $(u,v)\in\vec{E}$, and for $(v,u)$, with $f_{vu}=-f_{uv}$. For a flow vector $f$, the product $Af$ gives us the vector of divergences, i.e. for any $u\in V$

[TABLE]

We define the vector $b\in\mathbb{R}^{V}$ to encode the supplies, i.e. for $p\in P$ we set $b_{p}=\mu(p)$, and for $u\in V\setminus P$ we set $b_{u}=0$. It follows that the constraint $Af=b$ encodes (3)–(4). Finally, let us denote the cost of an edge $e=(u,v)\in\vec{E}$ by $c(e)=\|u-v\|_{2}$. Then, for the norm in the objective (5), we choose $\|f\|_{c}=\sum_{e\in E}{c(e)|f_{e}|}$. With these choices of $A$, $b$, and $\|\cdot\|=\|\cdot\|_{c}$, an optimal solution to (5)–(6) gives an optimal solution to (1)–(4). For the rest of this section we will fix $A$, $b$, $c$, and $\|\cdot\|_{c}$ to be as just defined.

Unfortunately, we cannot directly use Corollary 1 to get the running time we are aiming for, since the condition number of the matrix $A$ could be large. We address this by designing a preconditioner: another matrix $B$, of full column rank, which can be applied quickly, and has the property that $BA$ has small condition number. This allows us to apply Corollary 1 to the problem (5)–(6) with the modified, but equivalent, constraint $BAf=Bb$, and get a fast algorithm for the minimum cost flow problem in $G$.

In our construction of the preconditioner $B$ we will use the following lemma, which was implicit in [She17].

Lemma 4**.**

Let $\|\cdot\|_{X}$ be a norm on $\mathbb{R}^{n}$, and let $H$ be an $n\times m$ matrix. Suppose that there exists a $\gamma>0$ such that for any $h$ in the column span $Y$ of $H$, the norm $\|\cdot\|_{Y}$ satisfies

[TABLE]

Then $\kappa_{\|\cdot\|_{X}\to\|\cdot\|_{Y}}(H)\leq\gamma$.

Proof.

Let us identify $H$ with the linear map from $X=\mathbb{R}^{n}$ to $Y$ represented by it. Define $S(h)=\arg\min\{\|f\|_{X}:Hf=h\}$, with ties between minimizers broken arbitrarily. It is clear that for any $f$, $HS(Hf)=Hf$. By assumption, we have the inequalities

[TABLE]

which are equivalent to $\|H\|\leq 1$ and $\|S\|\leq\gamma$. ∎

We will design a preconditioner matrix $B$ such that $\|\tilde{b}\|_{1}$ approximates the cost of the minimum cost flow with supply vector $\tilde{b}$. Then the fact that $BA$ has small condition number will follow from Lemma 4. The following lemma captures the construction of $B$, which is inspired by embedding of the Earth Mover Distance in $\ell_{1}$ [Cha02, IT03].

Lemma 5**.**

There exists a matrix $B\in\mathbb{R}^{V\times V}$ of full column rank with at most $O(|V|\log\Delta)$ nonzero entries, such that the following holds. For any $\tilde{b}\in\mathbb{R}^{V}$ such that $\sum_{u\in V}{\tilde{b}_{u}}=0$, we have

[TABLE]

for $\gamma=O(\log(\Delta)/\varepsilon_{0})$. Moreover, a flow $f$ satisfying $Af=\tilde{b}$ of cost at most $\gamma\|B\tilde{b}\|_{1}$ can be constructed in time $O({n}{\varepsilon_{0}^{-d}}\log\Delta)$ given $B\tilde{b}$.

Proof.

The matrix $B$ has one row associated with each vertex $u\in V$. For $p\in P$, we let $B_{p,p}=\|p-N_{L}(p)\|_{2}$ and $B_{p,u}=0$ for any $u\neq p$. Every other vertex $u\in V\setminus P$ is a net point. Suppose that $u$ is in $N_{\ell}$; then for any $v\in V$ such that $N_{\ell}(v)=u$, $B_{u,v}=\frac{\varepsilon_{0}\Delta_{\ell}}{4(L+1)}$, and for all other $v$, $B_{u,v}=0$. This defines the matrix $B$. Notice that each vertex $u\in V$ contributes to at most $L+2$ nonzero entries of $B$: one for each $\ell\in\{0,\ldots,L\}$ corresponding to $N_{\ell}(u)$, and one more when $u\in P$. It is easy to verify directly that $B$ has full column rank.

The definition of $B$ guarantees that

[TABLE]

The remainder of the proof is broken down into several claims. For the first claim, let us extend the definition of $\textsc{cost}(G,\cdot)$ to supplies which can be non-zero on net points as well by

[TABLE]

Claim 1**.**

$\|B\tilde{b}\|_{1}\leq\textsc{cost}(G,\tilde{b})$.

Proof.

In any feasible flow $f$, the total cost of the flow on edges incident to points $p\in P$ must equal the first term in (7). Then, we just need to show that the remaining terms are a lower bound on the total cost of the flow on the remaining edges. In fact we will show that, for each $\ell$,

[TABLE]

It follows that the average of the term on the left over all $\ell$ is also a lower bound on the right hand side, which proves the claim. To establish inequality (8) notice that, in any feasible $f$, flow of value at least $\left|\sum_{u\in\widetilde{C}\cap V}\tilde{b}_{u}\right|$ must enter or leave each $\widetilde{C}\in\mathcal{\widetilde{C}}_{\ell}$, and every edge $e$ from a vertex in $\widetilde{C}$ to a vertex outside of $\widetilde{C}$ has cost at least $c(e)\geq\varepsilon_{0}\Delta_{\ell}/2$ by the definition of $G$. Therefore, the total cost of flow leaving or entering $\widetilde{C}$ is at least $\frac{\varepsilon_{0}\Delta_{\ell}}{2}\left|\sum_{u\in\widetilde{C}\cap V}\tilde{b}_{u}\right|$. Adding up these terms over all $\widetilde{C}\in\mathcal{\widetilde{C}}_{\ell}$ accounts for the cost of any edge of $G$ at most twice, and (8) follows. ∎

Claim 2**.**

$\textsc{cost}(G,\tilde{b})\leq\gamma\|B\tilde{b}\|_{1}$* for $\gamma=O(L/\varepsilon_{0})=O(\log(\Delta)/\varepsilon_{0})$.*

Proof.

We prove the claim by constructing an appropriate feasible flow $f$. We proceed to construct $f$ by levels, from the bottom up. The construction will be efficient, also proving the claim after “moreover” in the statement of the lemma. On level $L$, for any $p\in P$, we set $f_{pu}=\tilde{b}_{p}$, where $u=N_{L}(p)$. Then, for levels $\ell=L-1,L-2,\ldots,1$, we execute the following procedure in every subcell $\widetilde{C}\in\mathcal{\widetilde{C}}_{\ell-1}$. Let us define, for any node $u$ and the current flow $f$, the surplus

[TABLE]

I.e. this is how much more flow leaves $u$ than should, as prescribed by $\tilde{b}$. Pick any two $u$ and $v$ in $N_{\ell}\cap\widetilde{C}$ such that $\delta(u,f)>0>\delta(v,f)$, and add $\min\{|\delta(u,f)|,|\delta(v,f)|\}$ units of flow to $f_{vu}$. Continue until we have that for all $u\in N_{\ell}\cap\widetilde{C}$ the surpluses $\delta(u,f)$ have the same sign. Since at every step of this procedure we make the surplus of at least one node [math], and we always decrease the absolute value of the surplus at any node, we will stop after at most $|N_{\ell}\cap\widetilde{C}|$ steps. Finally, for the node $u\in N_{\ell-1}$ which is the center of $\widetilde{C}$, and for all nodes $v\in N_{\ell}\cap\widetilde{C}$ for which $\delta(v,f)\neq 0$, we set $f_{uv}=\delta(v,f)$. After this final step, every node $u\in N_{\ell}\cap\widetilde{C}$ has surplus [math].

For $\ell=0$, we perform essentially the same procedure, but in the entire cell $C_{0}$. I.e. we pick any two $u$ and $v$ in $N_{\ell}\cap C_{0}$ such that $\delta(u,f)>0>\delta(v,f)$, and add $\min\{|\delta(u,f)|,|\delta(v,f)|\}$ units of flow to $f_{vu}$. Once again, after at most $|N_{0}|$ steps all surpluses will have the same sign, and, as we show below, will in fact be [math], so $f$ will be feasible.

An easy induction argument shows that, once we have processed levels $L,L-1,\ldots,\ell$, for any $\ell>0$, the surplus at any node $u\in N_{L}\cup\ldots\cup N_{\ell}$ is [math], and the same is true for the surplus at any $p\in P$. To show that the flow $f$ is feasible, it is enough to show that after processing the cell $C_{0}$ on level [math], all nodes have surplus [math]. Notice that, while processing any subcell $\widetilde{C}$, we do not change the total surplus $\sum_{u\in V\cap\widetilde{C}}{\delta(u,f)}$. This means that, for any $\ell>0$, after having processed a subcell $\widetilde{C}\in\mathcal{\widetilde{C}}_{\ell-1}$ with center $u$, we have

[TABLE]

In particular, we have that before we start processing $C_{0}$, $\sum_{u\in N_{0}}{\delta(u,f)}=-\sum_{u\in V}{\tilde{b}_{u}}=0$. Since, once again, every step we make during the processing of $C_{0}$ preserves the total surplus, and we stop when all surpluses have the same sign, then at the time we are done it must be the case that every node has surplus [math], and, therefore, $f$ is a feasible flow.

Next we bound the cost of the flow $f$ constructed above. The first term on the right hand side of (7) accounts exactly for the cost of the flow on edges incident on points $p\in P$, as we already observed. Let $0<\ell\leq L$. During each step executed while processing the subcell $\widetilde{C}\in\mathcal{\widetilde{C}}_{\ell-1}$, we pick some $u\in\widetilde{C}\cap N_{\ell}$, and add flow of value $\delta(u,f)$ along an edge of cost at most the diameter of $\widetilde{C}$, which is $\sqrt{d}\varepsilon_{0}\Delta_{\ell}$. Moreover, the absolute value of $\delta(u,f)$ never increases while processing $\widetilde{C}$, and by (9) it is bounded from above by $\left|\sum_{v:N_{\ell}(v)=u}{\tilde{b}_{u}}\right|$ before $\widetilde{C}$ is processed. Therefore, the cost added to the flow while processing $\widetilde{C}$ is bounded by

[TABLE]

Adding up these upper bounds on the cost over all subcells in $\mathcal{\widetilde{C}}_{\ell-1}$ gives us that the total cost incurred while processing subcells on level $\ell-1$ is at most

[TABLE]

An analogous argument shows that the total cost incurred while processing $C_{0}$ is at most

[TABLE]

Therefore, the cost of $f$ is bounded by

[TABLE]

where the inequality follows by a term by term comparison with (7). ∎

The two claims finish the proof of the lemma, together with the observation that the construction of the flow $f$ can be implemented recursively in time $O({n}{\varepsilon_{0}^{-d}}\log\Delta)$. ∎

Our main result for solving the minimum cost flow problem in $G$ follows.

Theorem 5**.**

A flow $f$ feasible for the problem (1)–(4), with cost at most $1+O(\varepsilon)$ factor larger than the optimal cost, can be computed in time

[TABLE]

Proof.

We reformulate the problem as minimizing $\|f\|_{c}$ over $f\in\mathbb{R}^{\vec{E}}$ satisfying $BAf=Bb$. Since $B$ has full column rank, the constraint is equivalent to the original constraint $Af=b$. Moreover, the time for applying $BA$ is bounded by the total number of nonzeros of $B$ and $A$, which is $O(|E|+|V|\log\Delta)$. By Lemma 5 the assumption of Lemma 4 is satisfied with $H=BA$, $h=Bb$, $\|\cdot\|_{X}=\|\cdot\|_{c}$ and $\|\cdot\|_{Y}=\|\cdot\|_{1}$. Therefore, $\kappa_{\|\cdot\|_{c}\to\|\cdot\|_{1}}(H)\leq\gamma=O(\log(\Delta)/\varepsilon_{0})$. This allows us to apply Corollary 1 and get a $(1+\varepsilon,\varepsilon\gamma^{-2})$-solver with respect to $\|\cdot\|_{1}$ with the running time claimed in the theorem. To get a truly feasible flow $f$, we compose this solver with the algorithm in Lemma 5 which, for a supplies vector $\tilde{b}$, computes a feasible flow with cost bounded by $\gamma\|B\tilde{b}\|_{1}\leq\gamma\cdot\textsc{cost}(G,\tilde{b})$, and is, therefore, a $(\gamma,0)$-solver. By Theorem 4, this composed algorithm is a $(1+2\varepsilon,0)$-solver. The total running time is dominated by the solver guaranteed by Corollary 1. ∎

5 Generating a Transportation Map

Theorem 5 guarantees we can obtain an approximately optimal flow in our graph $G$ in nearly linear time. We wish to turn this flow into a transportation map on $P$. To accomplish this, we will repeatedly transform our flow into other flows, without increasing the cost, and ending at a flow which is also a transportation map. We will no longer keep the flow supported on the edges of $E$, and will instead allow positive flow between arbitrary pairs of points. We will gradually make sure that there is positive flow only between points in $P$ (as opposed to net points).

We first begin by defining a notion we call uniform flow parity.

Definition 4**.**

Given a flow $f\in\mathbb{R}^{V\times V}$, and a vertex $u\in V$, $f$ is said to satisfy, or have uniform flow parity at $u$ if for all remaining $v\in V$, $f_{uv}$ has the same sign. In other words, either there is flow going out of $v$ or there is flow going into $v$, but there is no flow passing through $v$. The flow $f$ is said to satisfy uniform flow parity if $f$ has uniform flow parity on every $v\in V$.

The next easy lemma shows that a flow that satisfies uniform flow parity is supported on edges between vertices with non-zero divergence.

Lemma 6**.**

Suppose that a flow $f\in\mathbb{R}^{V\times V}$ satisfies uniform flow parity at a vertex $u\in V$, and that $\sum_{v\in V}{f_{uv}}=0$. Then $f_{uv}=0$ for all $v\in V$.

Proof.

Assume, towards contradiction, that $f$ has [math] divergence at $u$ but $f_{uv}\neq 0$ for some $v\in V$. Then there must be some $w\in V$ such that $f_{uw}$ has the opposite sign to $f_{uv}$, or $\sum_{v\in V}{f_{uv}}=0$ would not hold. This contradicts uniform flow parity. ∎

We apply Lemma 6 through the following corollary of it.

Corollary 2**.**

Let $f$ be a flow on $V\times V$. Suppose that for any $u\in V\setminus P$, $\sum_{v\in V}{f_{uv}}=0$ and for any $p\in P$, $\sum_{v\in V}{f_{pv}}=\mu(p)$. Then, if $f$ satisfies uniform flow parity, it is a transportation map for $P$ and $\mu$.

Proof.

By Lemma 6, $f$ is supported on $P\times P$. Moreover, since for any $p\in P_{+}$ we have $\sum_{q\in P}{f_{pq}}=\mu(p)>0$, and $f$ satisfies uniform flow parity, it must be the case that $f_{pq}\geq 0$ for all $q\in P$. Similarly, for any $p\in P_{-}$ it must be the case that $f_{pq}\leq 0$ for all $q\in P$, or, equivalently, $f_{qp}\geq 0$ for all $q\in P$. It follows that $f$ is in fact supported on $P_{+}\times P_{-}$, and is, therefore, a transportation map. ∎

Let us take a flow $f$ which is an approximately optimal solution to (1)–(4). We can extend $f$ to $V\times V$ by setting $f_{uv}=0$ for $(u,v)\not\in E$. If we can transform $f$ into another flow $f^{\prime}$ without increasing its cost, so that $f^{\prime}$ satisfies uniform flow parity, then we can use $f^{\prime}$ as an approximately optimal transportation map. Towards this goal, we define the Cancellation Procedure Cancel-Vertex($f,u$), given in Algorithm 1, and analyzed in the following lemma.

Lemma 7**.**

Let $f$ be a flow on $V\times V$, and $u\in V$. The while loop in Cancel-Vertex($f,u$) makes at most $|\{v:f_{uv}\neq 0\}|$ many iterations, and, letting $f^{\prime}$ be the flow after the procedure is called, we have the following properties:

All divergences are preserved, i.e.

[TABLE]

for all $v\in V$. 2. 2.

The flow $f^{\prime}$ satisfies uniform flow parity at $u$. 3. 3.

If $f$ satisfies uniform flow parity at some vertex $v\in V$, then so does $f^{\prime}$. 4. 4.

The size of the support of $f$ is at most that of $f^{\prime}$. 5. 5.

The cost of $f^{\prime}$ with respect to the cost function $c(u,v)=\|u-v\|_{2}$ is at most the cost of $f$.

Proof.

To bound the number of iterations of the while loop, observe that every iteration sets the flow of at least one edge incident on $u$ to [math], and that the flow on an edge incident on $u$ is modified only if it nonzero. We proceed to prove the other claims.

The first claim follows because every iteration of the while loop of the procedure preserves the divergence $\sum_{w\in V}{f_{zw}}$ at every node $z\in V$: if $z=u$, then we remove $x$ units of flow going into $u$, and $x$ units of flow leaving $u$; otherwise, if $z=v$ or $z=w$ we just redirect $x$ units of flow entering or leaving $v$ to another vertex; finally, if $z$ is any other vertex, then the flow through it is unaffected.

The fact that $f^{\prime}$ satisfies uniform flow parity on $u$ follows from the termination condition of the loop.

To prove that uniform flow parity is preserved, again we show that this is the case for any iteration of the while loop. For $u$, we observe that the the flow is unchanged if $u$ already satisfies uniform flow parity. The only other vertices for which the flow on incident edges can be modified during an iteration are $v$ and $w$. For $v$, we have that $f_{vu}>0$, so, if $v$ has uniform flow parity, then the flow on every edge leaving $v$ must be non-negative. We decrease the value of $f_{vu}$ without making it negative (by the definition of $x$), and we increase the value of $f_{vw}$, which was non-negative to begin with. Therefore, uniform flow parity is preserved for $v$. The argument for $w$ is analogous.

To see that that the number of edges in the support of $f$ does not increase, notice that every iteration of the while loop removes at least one edge from the support of $f$, and adds flow to at most one edge which previously had flow of [math], namely $(v,w)$.

Finally, to show that the cost of the flow does not increase, observe that every iteration of the while loop replaces $x$ units of flow of total cost $x(\|v-u\|_{2}+\|w-u\|_{2})$ by $x$ units of flow of total cost $x\|w-v\|_{2}$, which is no larger than the original cost by the triangle inequality. ∎

Suppose that we maintain $f$ in a sparse representation; namely, we keep an adjacency list of the edges on which $f$ is not zero, with the corresponding flow values. Then, after preprocessing the adjacency list of $u$ in linear time to find the edges with positive and negative flow, every iteration of the while loop in the Cancellation Procedure can be executed in constant time, and the total running time, by Lemma 7, is bounded by $O(\deg_{f}(u))$, where $\deg_{f}(u)=|\{v:f_{uv}\neq 0\}|$.

Corollary 2 and Lemma 7 imply that if we apply the Cancellation Procedure to the flow $f$ and every vertex of the graph $G$, then the resulting flow has cost no greater than that of $f$, is supported on a set of edges of size bounded by $|E|$, and is a transportation map. In the next theorem, which is our main result for constructing a transportation map from a flow $f$ in $G$, we show that if we apply the procedure first to netpoints in $N_{L}$, then to $N_{L-1}$, etc., then the total running time is nearly linear.

Theorem 6**.**

There exists an algorithm running in time $O\left({n\varepsilon^{-d}\log(\Delta)}+{\varepsilon^{-2d}}{\log(\Delta)}\right)$ that, given as input a flow $f$ which is feasible for the minimum cost flow problem defined in Section 3.3, outputs a transportation map for $P$ and $\mu$ of cost no larger than that of $f$.

Proof.

As discussed above, we (mentally) extend $f$ to all of $V\times V$ by defining $f_{uv}=0$ if $(u,v)\not\in E$. We call Cancel-Vertex($f,u$) on $f$ and every net point $u\in N_{L}$ in arbitrary order; then we call it on every net point $u\in N_{L-1}$, and so on, finally calling it on every net point in $N_{0}$. By Lemma 7, after these calls uniform flow parity must hold at all net points. Moreover, it also holds at points in $P$, since in $G$ they only have a single edge incident on them, so initially uniform flow parity holds trivially, and by Lemma 7 it is not destroyed by later calls to the Cancellation Procedure. Thus, the flow at termination has uniform flow parity at all vertices, and, by Corollary 2, is a transportation map. Since calls to the Cancellation Procedure do not increase the cost of $f$, the cost of the transportation map is at most the cost of the initial flow.

Let us then consider a time step when the Cancellation Procedure has been applied to all net points in $N_{\ell+1}\cup\ldots\cup N_{L}$ and some points of $N_{\ell}$, but no points in $N_{0}\cup\ldots\cup N_{\ell-1}$. By Lemma 6, the flow $f$ at such a step is supported on edges incident on $N_{0}\cup\ldots\cup N_{\ell}\cup P$. Recall that in $G$ net points on level $\ell$ only connect to other net points on the same level and in the same cell, and to parent net points on level $\ell-1$. Then we have that for any $C\in\mathcal{C}_{\ell}$ and the cell $C^{\prime}\in\mathcal{C}_{\ell-1}$ containing it, and for any point $u\in N_{\ell}\cap C$, $f_{uv}$ can be nonzero only for $v\in(N_{\ell}\cap C)\cup(P\cap C)\cup(N_{\ell-1}\cap C^{\prime})\}$. Let $n_{\text{below}}(C)=\lvert P\cap C\rvert$, and recall the notation $\deg_{f}(u)=|\{v\in V:f_{uv}\neq 0\}|$. We then have that $\deg_{f}(u)\leq n_{\text{below}}(C)+{2}{\varepsilon^{-d}}$, which, by Lemma 7, also dominates the running time of the Cancellation Procedure on $u$. The total complexity of running the Cancellation Procedure on all points in $N_{\ell}$ is then dominated by

[TABLE]

Summing the running times over the $L+1=O(\log\Delta)$ levels finishes the proof. ∎

Combining Theorem 2, used with $\varepsilon_{0}$ set to a sufficiently small multiple of $\varepsilon/L$, and Theorems 5 and 6 gives the first claim of Theorem 1, but with the approximation holding in expectation. At the cost of increasing $\varepsilon$ by a factor of $2$, Markov’s inequality shows that the approximation also holds with probability at least $1/2$. Then the second statement in Theorem 1 follows from the first one, and Lemma 1.

Appendix A Proof of Lemma 3

The first claim follows from the triangle inequality, since $\mathrm{dist}_{G}(p,q)$ is the length of a polygonal path between $p$ and $q$, and $\|p-q\|_{2}$ is the length of the straightline segment between them.

For the second claim, it is enough to show that the canonical path between $p$ and $q$ has the required length. Observe first that, if $u=N_{\ell-1}(v)$ is the parent net point of some $v\in N_{\ell}$, then $\|u-v\|_{2}=\frac{\sqrt{d}\varepsilon_{0}\Delta_{\ell}}{2}$. Similarly, if $r\in P$, then $\|r-N_{\ell}(r)\|_{2}\leq\frac{\sqrt{d}\varepsilon_{0}\Delta_{\ell}}{2}$. This shows that

[TABLE]

Using the canonical path to bound the shortest path distance between $p$ and $q$ in $G$, we see that

[TABLE]

This shows that, conditional on ${\ell(p,q)}$, $\mathrm{dist}_{G}(p,q)=\|p-q\|_{2}+O(\varepsilon_{0}\Delta_{\ell(p,q)})$. All that remains is to find the expected value of $\mathrm{dist}_{G}(p,q)$ by taking expectations over the distribution of ${\ell(p,q)}$. We get that the expected value of $\mathrm{dist}_{G}(p,q)$ is

[TABLE]

The final inequality above follows from Lemma 2. This completes the proof of the lemma.

Appendix B Sherman’s Algorithms

To make our paper more self-contained, in this section we briefly sketch a version of the algorithm from Theorem 3 which uses the multiplicative weight update method [AHK12] in a simple way, but achieves a slightly worse running time than the one stated in the theorem. We do not claim any novelty, and everything in this section is due to Sherman.

One of Sherman’s main results is an analysis of how composing rough solvers affects the error in the linear constraints. In particular, suppose that $\mathcal{F}$ is a $(1+\varepsilon,\frac{\varepsilon}{2\kappa})$-solver, $\mathcal{G}$ is $(2,\frac{1}{2\kappa})$-solver, and $\mathcal{G}^{t}$ is $\mathcal{G}$ composed with itself $t$ times. (The notion of composition here is the same one used in Theorem 4.) Then, the composition $\mathcal{G}^{t}\circ\mathcal{F}$ is a $(1+O(\varepsilon),\frac{\varepsilon 2^{-t-1}}{\kappa})$-solver. This reduces designing the solver from Theorem 3 to designing a $(1+\varepsilon,\varepsilon)$-solver that performs $O(\varepsilon^{-2}\log(m))$ multiplications with $A$ or $A^{\top}$ (and possibly lower-order operations). We will sketch how to get a solver whose complexity is instead $O(\varepsilon^{-2}\log(m)\log_{1+\varepsilon}(\kappa))$ multiplications with $A$ or $A^{\top}$.

By scaling, we can assume that $\|A\|_{\|\cdot\|_{1}\to\|\cdot\|_{1}}\leq 1$, i.e. every column of the matrix $A$ has $\ell_{1}$-norm bounded by $1$. It then follows that for any $f$ that satisfies the constraint $Af=b$, we have $\|f\|_{1}\geq\|b\|_{1}$, i.e. the optimal value of the optimization problem (5)–(6) is at least $\|b\|_{1}$. Moreover, if the condition number of $A$ is $\kappa$, then there exists a map $S$ from the span of $A$, equipped with $\ell_{1}$, to $\mathbb{R}^{m}$, equipped with $\ell_{1}$, such that $\|S\|\leq\kappa$ and $AS(b)=b$. It follows that $\|S(b)\|_{1}\leq\kappa\|b\|_{1}$, so the optimal value lies in $[\|b\|_{1},\kappa\|b\|_{1}]$. We can then guess a value in this range (up to a factor of $1+\varepsilon$), and reduce the problem to making $O(\log_{1+\varepsilon}\kappa)$ calls to an algorithm solving the feasibility problem $\{f:\|f\|_{1}\leq t,Ax=b\}\neq\emptyset$. This is the feasibility problem that we will approximately solve using the multiplicative weights framework.

In fact, it will be convenient to reformulate the above feasibility problem as $\{f:\|f\|_{1}\leq 1,Ax=b/t\}\neq\emptyset$. In order to get a $(1+\varepsilon,O(\varepsilon))$-solver, it is enough to output an $f$ such that $\|f\|_{1}\leq 1$ and $\|Af-\frac{1}{t}b\|_{1}\leq\varepsilon$ whenever the problem is feasible. We will write $f=f^{+}-f^{-}$, where both $f^{+}$ and $f^{-}$ have non-negative coordinates and $\sum_{i=1}^{m}{f^{+}_{i}}=\sum_{i=1}^{m}{f^{-}_{i}}=1$. We initialize $f^{+}$ and $f^{-}$ to $f^{+}_{i}=f^{-}_{i}=\frac{1}{m}$ for all $i$, and we update them iteratively using the multiplicative weights update method. In each step of the algorithm, we either have $\|Af-\frac{1}{t}b\|_{1}\leq\varepsilon$ and we are done, or we can find a vector $y\in\{-1,1\}^{n}$ such that $y^{\top}(Af-b)>\varepsilon$. We update

[TABLE]

where $\eta$ is a parameter. If we choose $\eta=\frac{\varepsilon}{2}$, a standard analysis of the multiplicative weights method [AHK12] shows that, after $O(\varepsilon^{-2}\log m)$ iterations, we would either have $\|A(f^{+}-f^{-})-\frac{1}{t}b\|_{1}\leq\varepsilon$, as desired, or the average $\bar{y}$ of all vectors $y$ computed so far satisfies $\pm\bar{y}^{\top}a_{i}>\bar{y}^{\top}b$ for every column $a_{i}$ of $A$, certifying that the problem is infeasible.