Efficient Algorithms for Geometric Partial Matching

TL;DR

This paper introduces efficient algorithms for the geometric partial matching problem in the plane, providing exact and approximate solutions with improved runtime using primal-dual flow augmentation and dynamic data structures.

Contribution

The paper presents novel algorithms for geometric partial matching that achieve faster runtimes for exact and approximate solutions, applicable to any $L_p$-norm matching objectives.

Findings

Exact algorithm runs in $O((n + k^2) ext{polylog} n)$ time.

Approximate $(1 + ext{varepsilon})$-algorithm runs in $O((n + k extsqrt{k}) ext{polylog} n imes ext{log} extvarepsilon^{-1})$ time.

Additional exact algorithm for planar transportation problem with $O( extminig ^2, r n^{3/2}ig) ext{polylog} n$ runtime.

Abstract

Let $A$ and $B$ be two point sets in the plane of sizes $r$ and $n$ respectively (assume $r \leq n$), and let $k$ be a parameter. A matching between $A$ and $B$ is a family of pairs in $A \times B$ so that any point of $A \cup B$ appears in at most one pair. Given two positive integers $p$ and $q$, we define the cost of matching $M$ to be $c(M) = \sum_{(a, b) \in M}\|{a-b}\|_p^q$ where $\|{\cdot}\|_p$ is the $L_p$-norm. The geometric partial matching problem asks to find the minimum-cost size-$k$ matching between $A$ and $B$. We present efficient algorithms for geometric partial matching problem that work for any powers of $L_p$-norm matching objective: An exact algorithm that runs in $O((n + k^2) {\mathop{\mathrm{polylog}}} n)$ time, and a $(1 + \varepsilon)$-approximation algorithm that runs in $O((n + k\sqrt{k}) {\mathop{\mathrm{polylog}}} n \cdot \log\varepsilon^{-1})$ time. Both algorithms are based on the primal-dual flow augmentation scheme; the main improvements involve using dynamic data structures to achieve efficient flow augmentations. With similar techniques, we give an exact algorithm for the planar transportation problem running in $O(\min\{n^2, rn^{3/2}\} {\mathop{\mathrm{polylog}}} n)$ time.

Equations

$${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$c(\tau)$}}}\coloneqq\sum_{(a,b)\in A\times B}c(a,b)\cdot\tau(a,b).$$

$${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$\phi_{f}(v)$}}}\coloneqq\phi(v)+\sum_{w\mathord{\shortrightarrow}v\in\vec{E}}{f(w\mathord{\shortrightarrow}v)}-\sum_{v\mathord{\shortrightarrow}w\in\vec{E}}{f(v\mathord{\shortrightarrow}w)}.$$

$${\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{$c(f)$}}}\coloneqq\sum_{v\mathord{\shortrightarrow}w\in\operatorname{supp}(f)}c(v\mathord{\shortrightarrow}w)\cdot f(v\mathord{\shortrightarrow}w).$$

$${\vec{E}}_f:=\{v\text{shortrightarrow}w\mid f(v\text{shortrightarrow}w)=0\}\cup \{w\text{shortrightarrow}v\mid f(v\text{shortrightarrow}w)=1\}.$$

$$f^{\prime }(v\text{shortrightarrow}w)=\{\begin{array}{ll}0& w\text{shortrightarrow}v\in {\vec{E}}_f\text{ and }g(w\text{shortrightarrow}v)=1,\\ 1& v\text{shortrightarrow}w\in {\vec{E}}_f\text{ and }g(v\text{shortrightarrow}w)=1,\\ f(v\text{shortrightarrow}w)& \text{otherwise.}\end{array}$$

$$c_{\pi }(v\text{shortrightarrow}w):=c(v\text{shortrightarrow}w)-\pi (v)+\pi (w).$$

Full text

We present efficient algorithms for geometric partial matching problem that work for any powers of $L_{p}$-norm matching objective: An exact algorithm that runs in $O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time, and a $(1+\varepsilon)$-approximation algorithm that runs in $O((n+k\sqrt{k})\mathop{\mathrm{polylog}}n\cdot\log\varepsilon^{-1})$ time. Both algorithms are based on the primal-dual flow augmentation scheme; the main improvements involve using dynamic data structures to achieve efficient flow augmentations. With similar techniques, we give an exact algorithm for the planar transportation problem running in $O(\min\{n^{2},rn^{3/2}\}\mathop{\mathrm{polylog}}n)$ time.

1 Introduction

Given two point sets $A$ and $B$ in the plane, we consider the problem of finding the minimum-cost partial matching between $A$ and $B$. Formally, suppose $A$ has size $r$ and $B$ has size $n$ where $r\leq n$. Let $G(A,B)$ be the undirected complete bipartite graph between $A$ and $B$, and let the cost of edge $(a,b)$ be {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{c(a,b)}}}=\mathopen{\|}a-b\mathclose{\|}_{p}^{q}, for some positive integers $p$ and $q$. A matching $M$ in $G(A,B)$ is a set of edges sharing no endpoints. The size of $M$ is the number of edges in $M$. The cost of matching $M$, denoted $c(M)$, is defined to be the sum of costs of edges in $M$. For a parameter $k$, the problem of finding the minimum-cost size-$k$ matching in $G(A,B)$ is called the geometric partial matching problem. We call the corresponding problem in general bipartite graphs (with arbitrary edge costs) the partial matching problem.111Partial matching is also called imperfect matching or imperfect assignment [17, 9].

We also consider the following generalization of bipartite matching. Let $\phi:A\cup B\to\mathbb{Z}$ be an integral supply-demand function with positive value on points of $A$ and negative value on points of $B$, satisfying $\sum_{a\in A}\phi(a)=-\sum_{b\in B}\phi(b)$. Let {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{U}}}\coloneqq\max_{p\in A\cup B}\mathopen{|}\phi(p)\mathclose{|}. A transportation map is a function $\tau:A\times B\to\mathbb{R}_{\geq 0}$ such that $\sum_{b\in B}\tau(a,b)=\phi(a)$ for all $a\in A$ and $\sum_{a\in A}\tau(a,b)=-\phi(b)$ for all $b\in B$. We define the cost of $\tau$ to be

[TABLE]

The transportation problem asks to compute a transportation map of minimum cost.

Related work.

Maximum-size bipartite matching is a classical problem in graph algorithms. Upper bounds include the $O(m\sqrt{n})$ time algorithm by Hopcroft and Karp [10] and the $O(m\min\{\sqrt{m},n^{2/3}\})$ time algorithm by Even and Tarjan [7], where $n$ is the number of nodes and $m$ is the number of edges. The first improvement in over thirty years was made by Mądry [15], which uses an interior-point algorithm, runs in $O(m^{10/7}\mathop{\mathrm{polylog}}n)$ time.

The Hungarian algorithm [13] computes a minimum-cost maximum matching in a bipartite graph in roughly $O(mn)$ time. Faster algorithms have been developed, such as the $O(m\sqrt{n}\log(nC))$ time algorithms by Gabow and Tarjan [8] and the improved $O(m\sqrt{n}\log C)$ time algorithm by Duan et al. [6] assuming the edge costs are integral; here $C$ is the maximum cost of an edge. Ramshaw and Tarjan [17] showed that the Hungarian algorithm can be extended to compute a minimum-cost partial matching of size $k$ in $O(km+k^{2}\log r)$ time, where $r$ is the size of the smaller side of the bipartite graph. They also proposed a cost-scaling algorithm for partial matching that runs in time $O(m\sqrt{k}\log(kC))$, again assuming that costs are integral. By reduction to unit-capacity min-cost flow, Goldberg et al. [9] developed a cost-scaling algorithm for partial matching with an identical running time $O(m\sqrt{k}\log(kC))$, again only for integral edge costs.

In geometric settings, the Hungarian algorithm can be implemented to compute an optimal perfect matching between $A$ and $B$ (assuming equal size) in time $O(n^{2}\mathop{\mathrm{polylog}}n)$ [11] (see also [21, 1]). This algorithm computes an optimal size-$k$ matching in time $O(kn\mathop{\mathrm{polylog}}n)$. Faster approximation algorithms have been developed for computing perfect matchings in geometric settings [21, 22, 4, 18]. Recall that the cost of the edges are the $q$th power of their $L_{p}$-distances. When $q=1$, the best algorithm to date by Sharathkumar and Agarwal [19] computes $(1+\varepsilon)$-approximation to the value of optimal perfect matching in $O(n\mathop{\mathrm{polylog}}n\cdot\mathop{\mathrm{poly}}\varepsilon^{-1})$ expected time with high probability. Their algorithm can also compute a $(1+\varepsilon)$-approximate partial matching within the same time bound. For $q>1$, the best known approximation algorithm to compute a perfect matching runs in $O(n^{3/2}\mathop{\mathrm{polylog}}n\log(1/\varepsilon))$ time [18]; it is not obvious how to extend this algorithm to the partial matching setting.

The transportation problem can also be formulated as an instance of the minimum-cost flow problem. The strongly polynomial uncapacitated min-cost flow algorithm by Orlin [16] solves the transportation problem in $O((m+n\log n)n\log n)$ time. Lee and Sidford [14] give a weakly polynomial algorithm that runs in $O(m\sqrt{n}\mathop{\mathrm{polylog}}(n,U))$ time, where $U$ is the maximum amount of node supply-demand. Agarwal et al. [2, 3] showed that Orlin’s algorithm can be implemented to solve 2D transportation in time $O(n^{2}\mathop{\mathrm{polylog}}n)$. In case of $O(1)$-dimension Euclidean space, by adapting the Lee-Sidford algorithm, they developed a $(1+\varepsilon)$-approximation algorithm that runs in $O(n^{3/2}\mathop{\mathrm{poly}}\varepsilon^{-1}\mathop{\mathrm{polylog}}(n,U))$ time. They also gave a Monte-Carlo algorithm that computes an $O(\log^{2}(1/\varepsilon))$-approximate solution in $O(n^{1+\varepsilon})$ time with high probability. Recently, Khesin, Niklov, and Paramonov [12] obtained a $(1+\varepsilon)$-approximation in low-dimensional Euclidean space that runs in randomized $O(n\mathop{\mathrm{poly}}\varepsilon^{-1}\mathop{\mathrm{polylog}}(n,U))$ time.

Our results.

There are three main results in this paper. First in Section 2 we present an efficient algorithm for computing an optimal partial matching in the plane.

Theorem 1.1**.**

Given two point sets $A$ and $B$ in the plane each of size at most $n$ and an integer $k\leq n$, a minimum-cost matching of size $k$ between $A$ and $B$ can be computed in $O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time.

We use bichromatic closest pair (BCP) data structures to implement the Hungarian algorithm efficiently, similar to Agarwal et al. [1] and Kaplan et al. [11]. But unlike their algorithms which take $\Omega(n)$ time to find an augmenting path, we show that after $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, an augmenting path can be found in $O(k\mathop{\mathrm{polylog}}n)$ time. The key is to recycle (rather than rebuild) our data structures from one augmentation to the next. We refer to this idea as the rewinding mechanism.

Next in Sections 3, we obtain a $(1+\varepsilon)$-approximation algorithm for the geometric partial matching problem in the plane by providing an efficient implementation of the unit-capacity min-cost flow algorithm by Goldberg et al. [9].

Theorem 1.2**.**

Given two point sets $A$ and $B$ in $\mathbb{R}^{2}$ each of size at most $n$, an integer $k\leq n$, and a parameter $\varepsilon>0$, a $(1+\varepsilon)$-approximate min-cost matching of size $k$ between $A$ and $B$ can be computed in $O((n+k\sqrt{k})\mathop{\mathrm{polylog}}n\cdot\log\varepsilon^{-1})$ time.

The main challenge here is how to deal with the dead nodes, which neither have excess/deficit nor have flow passing through them, but still contribute to the size of the graph. We show that the number of alive nodes is only $O(k)$, and then represent the dead nodes implicitly so that the Hungarian search and computation of a blocking flow can be implemented in $O(k\mathop{\mathrm{polylog}}n)$ time.

Finally in Section 4 we present a faster algorithm for the transportation problem in $\mathbb{R}^{2}$ when the two point sets are unbalanced.

Theorem 1.3**.**

Given two point sets $A$ and $B$ in $\mathbb{R}^{2}$ of sizes $r$ and $n$ respectively with $r\leq n$, along with supply-demand function $\phi:A\cup B\to\mathbb{Z}$, an optimal transportation map between $A$ and $B$ can be computed in $O(\min\{n^{2},rn^{3/2}\}\mathop{\mathrm{polylog}}n)$ time.

Our result improves over the $O(n^{2}\mathop{\mathrm{polylog}}n)$ time algorithm by Agarwal et al. [3] for $r=o(\sqrt{n})$. Similar to their algorithm, we also use the strongly polynomial uncapacitated minimum-cost flow algorithm by Orlin [16], but additional ideas are needed for efficient implementation. Unlike in the case of matchings, the support of the transportation problem may have size $\Omega(n)$ even when $r$ is a constant; so naïvely we can no longer spend time proportional to the size of support of the transportation map. However, with careful implementation we ensure that the support is acyclic, and one can find an augmenting path in $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ time with proper data structures, assuming $r\leq\sqrt{n}$.

2 Minimum-cost partial matchings using Hungarian algorithm

In this section, we solve the geometric partial matching problem and prove Theorem 1.1 by implementing the Hungarian algorithm for partial matching in $O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time.

A node $v$ is matched by matching $M$ if $v$ is the endpoint of some edge in $M$; otherwise $v$ is unmatched. Given a matching $M$, an augmenting path $\Pi=(a_{1},b_{1},\ldots,a_{\ell},b_{\ell})$ is an odd-length path with unmatched endpoints ($a_{1}$ and $b_{\ell}$) that alternates between edges outside and inside of $M$. The symmetric difference $M\oplus\Pi$ creates a new matching of size $\mathopen{|}M\mathclose{|}+1$, called the augmentation of $M$ by $\Pi$. The dual to the standard linear program for partial matching has one dual variable for each node $v$, called the potential $\pi(v)$ of $v$. Given potential $\pi$, we can define the reduced cost of the edges to be {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{c_{\pi}(v,w)}}}\coloneqq c(v,w)-\pi(v)+\pi(w). Potential $\pi$ is feasible on edge $(v,w)$ if $c_{\pi}(v,w)$ is nonnegative. Potential $\pi$ is feasible if $\pi$ are feasible on every edge in $G$. We say that an edge $(v,w)$ is admissible under potential $\pi$ if $c_{\pi}(v,w)=0$.

Fast implementation of Hungarian search.

The Hungarian algorithm is initialized with $M\leftarrow\emptyset$ and $\pi\leftarrow 0$. Each iteration of the Hungarian algorithm augments $M$ by an admissible augmenting path $\Pi$, discovered using a procedure called the Hungarian search. The algorithm terminates after $k$ augmentations, exactly when $\mathopen{|}M\mathclose{|}=k$; Ramshaw and Tarjan [17] showed that $M$ is guaranteed to be an optimal partial matching.

The Hungarian search grows a set of reachable nodes $X$ from all unmatched $v\in A$ using augmenting paths of admissible edges. Initially, $X$ is the set of unmatched nodes in $A$. Let the frontier of $X$ be the edges in $(A\cap X)\times(B\setminus X)$. $X$ is grown by relaxing an edge $(a,b)$ in the frontier: Add $b$ into $X$, modify potential to make $(a,b)$ admissible, preserve $c_{\pi}$ on other edges within $X$, and keep $\pi$ feasible on edges outside of $X$. Specifically, the algorithm relaxes the min-reduced-cost frontier edge $(a,b)$, and then raises $\pi(v)$ by $c_{\pi}(a,b)$ for all $v\in X$. If $b$ is already matched, then we also relax the matching edge $(a^{\prime},b)$ and add $a^{\prime}$ into $X$. The search finishes when $b$ is unmatched, and an admissible augmenting path now can be recovered.

In the geometric setting, we find the min-reduced-cost frontier edge using a dynamic bichromatic closest pair (BCP) data structure, similar to [3, 21]. Given two point sets $P$ and $Q$ in the plane and a weight function $\omega:P\cup Q\to\mathbb{R}$, the BCP is two points $a\in P$ and $b\in Q$ minimizing the additively weighted distance $c(a,b)-\omega(a)+\omega(b)$. Thus, a minimum reduced-cost frontier edge is precisely the BCP of point sets $P=A\cap X$ and $Q=B\setminus X$, with $\omega=\pi$. Note that the “state” of this BCP is parameterized by $X$ and $\pi$.

The dynamic BCP data structure by Kaplan et al. [11] supports point insertions and deletions in $O(\mathop{\mathrm{polylog}}n)$ time and answers queries in $O(\log^{2}n)$ time for our setting. Each relaxation in the Hungarian search requires one query, one deletion, and at most one insertion (aside from the potential updates). As $\mathopen{|}M\mathclose{|}\leq k$ throughout, there are at most $2k$ relaxations in each Hungarian search, and the BCP can be used to implement each Hungarian search in $O(k\mathop{\mathrm{polylog}}n)$ time.

Rewinding mechanism.

We observe that exactly one node of $A$ is newly matched after an augmentation. Thus (modulo potential changes), we can obtain the initial state of the BCP for the $(i+1)$-th Hungarian search from the $i$-th one with a single BCP deletion.

If we remember the sequence of points added to $X$ in the $i$-th Hungarian search, then at the start of the $(i+1)$-th Hungarian search we can rewind this sequence by applying the opposite insert/delete operation to each BCP update in reverse order to obtain the initial state of the $i$-th BCP. With one additional BCP deletion, we have the initial state of the $(i+1)$-th BCP. The number of insertions/deletions is bounded by the number of relaxations per Hungarian search which is $O(k)$. Therefore we can recover, in $O(k\mathop{\mathrm{polylog}}n)$ time, the initial BCP data structure for each Hungarian search beyond the first. We refer to this procedure as the rewinding mechanism.

Potential updates.

We modify a trick from Vaidya [21] to batch potential updates. Potential is tracked with a stored value $\gamma(v)$, while the true value of $\pi(v)$ may have changed since $\gamma(v)$ was last recorded. This is done by aggregating potential changes into a variable $\delta$, which is initially 0 at the very beginning of the algorithm. Whenever we would raise the potential of all nodes in $X$, we raise $\delta$ by that amount instead. We maintain the following invariant: $\pi(v)=\gamma(v)$ for $v\not\in X$, and $\pi(v)=\gamma(v)+\delta$ for $v\in X$.

At the beginning of the algorithm, $X$ is empty and stored values are equal to true values. When $a\in A$ is added to $X$, we update its stored value to $\pi(a)-\delta$ for the current value of $\delta$, and use that stored value as its BCP weight. Since the BCP weights are uniformly offset from $\pi(v)$ by $\delta$, the pair reported by the BCP is still minimum. When $b\in B$ is added to $X$, we update its stored value to $\pi(b)-\delta$ (although it won’t be added to a BCP set). When a node is removed from $X$ (e.g. by augmentation or rewinding), we update the stored potential $\gamma(v)\leftarrow\pi(v)+\delta$, again for the current value of $\delta$. Unlike Vaidya [21], we do not reset $\delta$ across Hungarian searches for the sake of rewinding.

There are $O(k)$ relaxations and thus $O(k)$ updates to $\delta$ per Hungarian search. $O(k)$ stored values are updated per rewinding, so the time spent on potential updates per Hungarian search is $O(k)$. Putting everything together, our implementation of the Hungarian algorithm runs in $O((n+k^{2})\mathop{\mathrm{polylog}}n)$ time. This proves Theorem 1.1.

3 Approximating min-cost partial matching through cost-scaling

In this section we describe an approximation algorithm for computing a min-cost partial matching. We reduce the problem to computing a min-cost circulation in a flow network (Section 3.1). We adapt the cost-scaling algorithm by Goldberg et al. [9] for computing min-cost flow of a unit-capacity network (Section 3.2). Finally, we show how their algorithm can be implemented in $O\left((n+k^{3/2})\mathop{\mathrm{polylog}}(n)\log(1/\varepsilon)\right)$ time in our setting (Section 3.3).

3.1 From matching to circulation

Given a bipartite graph $G$ with node sets $A$ and $B$, we construct a flow network $N=(V,\vec{E})$ in a standard way [17] so that a min-cost matching in $G$ corresponds to a min-cost integral circulation in $N$.

Flow network.

Each node in $G$ becomes a node in $N$ and each edge $(a,b)$ in $G$ becomes an arc $a\mathord{\shortrightarrow}b$ in $N$; we refer to these nodes and arcs as bipartite nodes and bipartite arcs. We also include a source node $s$ and sink node $t$ in $N$. For each $a\in A$, we add a left dummy arc $s\mathord{\shortrightarrow}a$ and for each $b\in B$ a right dummy arc $b\mathord{\shortrightarrow}t$. The cost $c(v\mathord{\shortrightarrow}w)$ is equal to $c(v,w)$ if $v\mathord{\shortrightarrow}w$ is a bipartite arc and [math] if $v\mathord{\shortrightarrow}w$ is a dummy arc. All arcs in $N$ have unit capacity.

Let $\phi:V\to\mathbb{Z}$ be an integral supply/demand function on nodes of $N$. The positive values of $\phi(v)$ are referred to as supply, and the negative values of $\phi(v)$ as demand. A pseudoflow $f:\vec{E}\to[0,1]$ is a function on arcs of $N$. The support of $f$ in $N$, denoted as $\operatorname{supp}(f)$, is the set of arcs with positive flow. Given a pseudoflow $f$, the imbalance of a node is

[TABLE]

We call positive imbalance excess and negative imbalance deficit. A node is balanced if it has zero imbalance. If all nodes are balanced, the pseudoflow is a circulation. The cost of a pseudoflow is defined to be

[TABLE]

The minimum-cost flow problem (MCF) asks to find a circulation of minimum cost.

If we set $\phi(s)=k$, $\phi(t)=k$, and $\phi(v)=0$ for all $v\in A\cup B$, then an integral circulation $f$ corresponds to a partial matching $M$ of size $k$ and vice versa. Moreover, $c(M)=c(f)$. Hence, the problem of computing a min-cost matching of size $k$ in $G$ transforms to computing an integral circulation in $N$. The following lemma will be useful for our algorithm.

Lemma 3.1**.**

Let $N$ be the network constructed from the bipartite graph $G$ above.

1. (i)

For any integral circulation $g$ in $N$, the size of $\operatorname{supp}(g)$ is at most $3k$. 2. (ii)

For any integral pseudoflow $f$ in $N$ with $O(k)$ excess, the size of $\operatorname{supp}(f)$ is $O(k)$.

3.2 A cost-scaling algorithm

Before describing the algorithm, we need to introduce a few more concepts.

Residual network and admissibility.

If $f$ is an integral pseudoflow on $N$ (that is, $f(v\mathord{\shortrightarrow}w)\in\{0,1\}$ for every arc in $\vec{E}$), then each arc $v\mathord{\shortrightarrow}w$ in $N$ is either idle with $f(v\mathord{\shortrightarrow}w)=0$ or saturated with $f(u\mathord{\shortrightarrow}v)=1$.

Given a pseudoflow $f$, the residual network $N_{f}=(V,\vec{E}_{f})$ is defined as follows. For each idle arc $v\mathord{\shortrightarrow}w$ in $\vec{E}$, we add a forward residual arc $v\mathord{\shortrightarrow}w$ in $N_{f}$. For each saturated arc $v\mathord{\shortrightarrow}w$ in $\vec{E}$, we add a backward residual arc $w\mathord{\shortrightarrow}v$ in $N_{f}$. The set of residual arcs in $N_{f}$ is therefore

[TABLE]

The cost of a forward residual arc $v\mathord{\shortrightarrow}w$ is $c(v\mathord{\shortrightarrow}w)$, while the cost of a backward residual arc $w\mathord{\shortrightarrow}v$ is $-c(v\mathord{\shortrightarrow}w)$. Each arc in $N_{f}$ also has unit capacity. By Lemma 3.1, $N_{f}$ has $O(k)$ backward arcs if $f$ has $O(k)$ excess.

A residual pseudoflow $g$ in $N_{f}$ can be used to change $f$ into a different pseudoflow on $N$ by augmentation. For simplicity, we only describe augmentation for the case where $f$ and $g$ are integral. Specifically, augmenting $f$ by $g$ produces a pseudoflow $f^{\prime}$ in $N$ where

[TABLE]

Using LP duality for min-cost flow, we assign potential $\pi(v)$ to each node $v$ in $N$. The reduced cost of an arc $v\mathord{\shortrightarrow}w$ in $N$ with respect to $\pi$ is defined as

[TABLE]

Similarly we define the reduced cost of arcs in $N_{f}$: the reduced cost of a forward residual arc $v\mathord{\shortrightarrow}w$ in $N_{f}$ is $c_{\pi}(v\mathord{\shortrightarrow}w)$, and the reduced cost of a backward residual arc $w\mathord{\shortrightarrow}v$ in $N_{f}$ is $-c_{\pi}(v\mathord{\shortrightarrow}w)$. Abusing the notation, we also use $c_{\pi}$ to denote the reduced cost of arcs in $N_{f}$.

The dual feasibility constraint asks that $c_{\pi}(v\mathord{\shortrightarrow}w)\geq 0$ holds for every arc $v\mathord{\shortrightarrow}w$ in $\vec{E}$; potential $\pi$ that satisfy this constraint is said to be feasible. Suppose we relax the dual feasibility constraint to allow some small violation in the value of $c_{\pi}(v\mathord{\shortrightarrow}w)$. We say that a pair of pseudoflow $f$ and potential $\pi$ is $\theta$-optimal [20, 5] if $c_{\pi}(v\mathord{\shortrightarrow}w)\geq-\theta$ for every residual arc $v\mathord{\shortrightarrow}w$ in $\vec{E}_{f}$. Pseudoflow $f$ is $\theta$-optimal if it is $\theta$-optimal with respect to some potential $\pi$; potential $\pi$ is $\theta$-optimal if it is $\theta$-optimal with respect to some pseudoflow $f$. Given a pseudoflow $f$ and potential $\pi$, a residual arc $v\mathord{\shortrightarrow}w$ in $\vec{E}_{f}$ is admissible if $c_{\pi}(v\mathord{\shortrightarrow}w)\leq 0$. We say that a pseudoflow $g$ in $N_{f}$ is admissible if $g(v\mathord{\shortrightarrow}w)>0$ only on admissible arcs $v\mathord{\shortrightarrow}w$, and $g(v\mathord{\shortrightarrow}w)=0$ otherwise.222The same admissibility/feasibility definitions will be used later in Section 4. However, the algorithm in Section 4 maintains a 0-optimal $f$ and therefore admissible residual arcs always have $c_{\pi}(v\mathord{\shortrightarrow}w)=0$. We will use the following well-known property of $\theta$-optimality.

Lemma 3.2**.**

Let $f$ be an $\theta$-optimal pseudoflow in $N$ and let $g$ be an admissible pseudoflow in $N_{f}$. Then $f$ augmented by $g$ is also $\theta$-optimal in $N$.

Using Lemma 3.1, the following lemma can be proved about $\theta$-optimality:

Lemma 3.3**.**

Let $f$ be a $\theta$-optimal integer circulation in $N$, and $f^{*}$ be an optimal integer circulation for $N$. Then, $c(f)\leq c(f^{*})+6k\theta$.

Estimating the value of $c(f^{*})$.

We now describe a procedure for estimating $c(f^{*})$ within a polynomial factor, which is used to initialize the cost-scaling algorithm.

Let $T$ be a minimum spanning tree of $A\cup B$ under the cost function $c$. Let $e_{1},e_{2},\ldots,e_{n-1}$ be the edges of $T$ sorted in nondecreasing order of length. Let $T_{i}$ be the forest consisting of the nodes of $A\cup B$ and edges $e_{1},\ldots,e_{i}$. We call a matching $M$ intra-cluster if both endpoints of each edge in $M$ lie in the same connected component of $T_{i}$. The following lemma will be used by our cost-scaling algorithm:

Lemma 3.4**.**

Let $A$ and $B$ be two point sets in the plane. Define $i^{*}$ to be the smallest index $i$ such that there is an intra-cluster matching of size $k$ in $T_{i^{*}}$. Set $\overline{\theta}$ $\coloneqq n^{q}\cdot c(e_{i^{*}})$. Then

1. (i)

The value of $i^{*}$ can be computed in $O(n\log n)$ time. 2. (ii)

$c(e_{i^{*}})\leq c(f^{*})\leq\overline{\theta}$. 3. (iii)

There is a $\overline{\theta}$-optimal circulation in the network $N$ with respect to the all-zero potential, assuming $\phi(s)=k$, $\phi(t)=-k$, and $\phi(v)=0$ for all $v\in A\cup B$.

As a consequence of Lemmas 3.4(ii) and 3.3, we have:

Corollary 3.5**.**

The cost of a $\underline{\theta}$-optimal integral circulation in $N$ is at most $(1+\varepsilon)c(f^{*})$, where {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{\underline{\theta}}}}\coloneqq\frac{\varepsilon}{6k}\cdot c(e_{i^{*}}).

Overview of the algorithm.

We are now ready to describe our algorithm, which closely follows Goldberg et al. [9]. The algorithm works in scales. In the beginning of each scale, we fix a cost scaling parameter $\theta$ and maintain potential $\pi$ with the following property:

($\ast$)

There exists a $2\theta$-optimal integral circulation in $N$ with respect to $\pi$.

For the initial scale, we set $\theta\leftarrow\overline{\theta}$ and $\pi\leftarrow 0$. By Lemma 3.4(iii), property (($\ast$) ‣ 3.2) is satisfied initially. Each scale of the algorithm consists of two stages. In the scale initialization stage, Scale-Init computes a $\theta$-optimal pseudoflow $f$. In the refinement stage, Refine converts $f$ into a $\theta$-optimal (integral) circulation $g$. In both stages, $\pi$ is updated as necessary. If $\theta\leq\underline{\theta}$, we return $g$. Otherwise, we set $\theta\leftarrow\theta/2$ and start the next scale. Note that property (($\ast$) ‣ 3.2) is satisfied in the beginning of each scale.

By Corollary 3.5, when the algorithm terminates, it returns an integral circulation $\tilde{f}$ in $N$ of cost at most $(1+\varepsilon)c(f^{*})$, which corresponds to a $(1+\varepsilon)$-approximate min-cost matching of size $k$ in $G$. The algorithm terminates in $\log_{2}(\overline{\theta}/\underline{\theta})=O(\log(n/\varepsilon))$ scales.

Scale initialization.

In the first scale, we compute a $\overline{\theta}$-optimal pseudoflow by simply setting $f(v\mathord{\shortrightarrow}w)\leftarrow 0$ for all arcs in $\vec{E}$. For subsequent scales, we begin with a $2\theta$-optimal circulation $f$ in $N$. First, we raise the potential of all nodes in $A$ by $\theta$, all nodes in $B$ by $2\theta$, and $t$ by $3\theta$. The potential of $s$ is unchanged. Such potential change increases the reduced cost of all forward arcs to at least $-\theta$.

Next, for each backward arc $w\mathord{\shortrightarrow}v$ in $N_{f}$ with $c_{\pi}(w\mathord{\shortrightarrow}v)<-\theta$, we set $f(v\mathord{\shortrightarrow}w)\leftarrow 0$ (that is, make arc $v\mathord{\shortrightarrow}w$ idle), which replaces the backward arc $w\mathord{\shortrightarrow}v$ in $N_{f}$ with forward arc $v\mathord{\shortrightarrow}w$ of positive reduced cost. After this step, the resulting pseudoflow must be $\theta$-optimal as all arcs of $N_{f}$ have reduced cost at least $-\theta$.

The desaturation of each backward arc creates one unit of excess. Since there are at most $3k$ backward arcs, the pseudoflow has at most $3k$ excess after Scale-Init. There are $O(n)$ potential updates and $O(k)$ arcs to desautrate, so the time required for Scale-Init is $O(n)$.

Refinement.

The procedure Refine converts a $\theta$-optimal pseudoflow with $O(k)$ excess into a $\theta$-optimal circulation, using a primal-dual augmentation algorithm. A path in $N_{f}$ is an augmenting path if it begins at an excess node and ends at a deficit node. We call an admissible pseudoflow $g$ in $N_{f}$ an admissible blocking flow if $g$ saturates at least one arc in every admissible augmenting path in $N_{g}$. In other words, there is no admissible excess-deficit path in the residual network after augmentation by $g$. Each iteration of Refine finds an admissible blocking flow to be added to the current pseudoflow in two steps:

Hungarian search: a Dijkstra-like search that begins at the set of excess nodes and raises potential until there is an excess-deficit path of admissible arcs in $N_{f}$. 2. 2.

Augmentation: construct an admissible blocking flow by performing depth-first search on the set of admissible arcs of $N_{f}$. It suffices to repeatedly extract admissible augmenting paths until no more admissible excess-deficit paths remain.

The algorithm repeats these steps until the total excess becomes zero. The following lemma bounds the number of iterations in the Refine procedure at each scale.

Lemma 3.6**.**

Let $\theta$ be the scaling parameter and $\pi_{0}$ the potential function at the beginning of a scale, such that there exists an integral $2\theta$-optimal circulation with respect to $\pi_{0}$. Let $f$ be a $\theta$-optimal pseudoflow with excess $O(k)$. Then Refine terminates within $O(\sqrt{k})$ iterations.

Proof.

We sketch the proof, which is adapted from Goldberg et al. [9]. Let $f_{0}$ be the assumed $2\theta$-optimal integral circulation with respect to $\pi_{0}$, and let $\pi$ be the potential maintained during Refine. Let {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{d(v)}}}\coloneqq(\pi(v)-\pi_{0}(v))/\theta, that is, the increase in potential at $v$ in units of $\theta$. We divide the iterations of Refine into two phases: before and after every (remaining) excess node has $d(v)\geq\sqrt{k}$. Each Hungarian search raises excess potential by at least $\theta$, since we use blocking flows. Thus, the first phase lasts at most $\sqrt{k}$ iterations.

At the start of the second phase, consider the set of arcs {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{E^{+}}}}\coloneqq\{v\mathord{\shortrightarrow}w\in\vec{E}\mid f(v\mathord{\shortrightarrow}w)<f_{0}(v\mathord{\shortrightarrow}w)\}. One can argue that the remaining excess with respect to $f$ is bounded above by the size of any cut separating the excess and deficit nodes [9, Lemma 4]. The proof examines cuts {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{Y_{i}}}}\coloneqq\{v\mid d(v)>i\} for $0\leq i\leq\sqrt{k}$. By $\theta$-optimality of $f$ and $2\theta$-optimality of $f_{0}$, one can show that each arc in $E^{+}$ crosses at most 3 cuts. Furthermore, the size of $E^{+}$ is $O(k)$, bounded by the support size of $f$ and $f_{0}$. Averaging, there is a cut among $Y_{i}$s of size at most $3k/\sqrt{k}$, so the total excess remaining is $O(\sqrt{k})$. Each iteration of Refine eliminates at least one unit of excess, so the number of second phase iterations is also at most $O(\sqrt{k})$. ∎

In the next subsection we show that after $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, an iteration of Refine can be performed in $O(k\mathop{\mathrm{polylog}}n)$ time (Lemma 3.8). By Lemma 3.6 and the fact the algorithm terminates in $O(\log(n/\varepsilon))$ scales, the overall running time of the algorithm is $O((n+k^{3/2})\mathop{\mathrm{polylog}}n\log(1/\varepsilon))$, as claimed in Theorem 1.2.

3.3 Fast implementation of refinement stage

We now describe a fast implementation of Refine. The Hungarian search and augmentation steps are similar: each traversing through the residual network using admissible arcs starting from the excess nodes. Due to lack of space, we only describe the former.

At a high level, let $X$ be the subset of nodes visited by the Hungarian search so far. Initially $X$ is the set of excess nodes. At each step, the algorithm finds a minimum-reduced-cost arc $v\mathord{\shortrightarrow}w$ in $N_{f}$ from $X$ to $V\setminus X$. If $v\mathord{\shortrightarrow}w$ is not admissible, the potential of all nodes in $X$ is increased by $\left\lceil c_{\pi}(v\mathord{\shortrightarrow}w)/\theta\right\rceil$ to make $v\mathord{\shortrightarrow}w$ admissible. If $w$ is a deficit node, the search terminates. Otherwise, $w$ is added to $X$ and the search continues.

Implementing the Hungarian search efficiently is more difficult than in Section 2 because (a) excess nodes may show up in $A$ as well as in $B$, (b) a balanced node may become imbalanced later in the scales, and (c) the potential of excess nodes may be non-uniform. We therefore need a more complex data structure.

We call a node $v$ of $N$ dead if $\phi_{f}(v)=0$ and no arc of $\operatorname{supp}(f)$ is incident to $v$; otherwise $v$ is alive. Note that $s$ and $t$ are always alive. Let ${A}^{\text{\char 96}}$ denote the set of alive nodes in $A$; define ${B}^{\text{\char 96}}$ similarly. There are only $O(k)$ alive nodes, as each can be charged to its adjoining $\operatorname{supp}(f)$ arcs or its imbalance. We treat alive and dead nodes separately to implement the Hungarian search efficiently. By definition, dead nodes only adjoin forward arcs in $N_{f}$. Thus, the in-degree (resp. out-degree) of a node in ${A}\setminus{A}^{\text{\char 96}}$ (resp. ${B}\setminus{B}^{\text{\char 96}}$) is 1, and any path passing through a dead node has a subpath of the form $s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}b$ for some $b\in B$ or $a\mathord{\shortrightarrow}v\mathord{\shortrightarrow}t$ for some $a\in A$. Consequently, a path in $N_{f}$ may have at most two consecutive dead nodes, and in the case of two consecutive dead nodes there is a subpath of the form $s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}w\mathord{\shortrightarrow}t$ where $v\in{A}\setminus{A}^{\text{\char 96}}$ and $w\in{B}\setminus{B}^{\text{\char 96}}$. We call such paths, from an alive node to an alive node with only dead interior nodes, alive paths. Let the reduced cost $c_{\pi}(\Pi)$ of an alive path $\Pi$ be the sum of $c_{\pi}$ over its arcs. We say $\Pi$ is weakly admissible if $c_{\pi}(\Pi)\leq 0$.

We find the min-reduced-cost alive path of lengths $1$, $2$, and $3$ leaving $X$, then relax the cheapest among them (raise potential of $X$ by $c(\Pi)$ and add every node of $\Pi$ into $X$). Essentially, relaxing alive paths “skips over” dead nodes. Since reduced costs telescope on paths, weak admissibility of an alive path depends only on the potential of its alive endpoints. Thus, we can query the minimum alive path using a partial assignment of $\pi$ on only the alive nodes, leaving $\pi$ over the dead nodes untracked. We now describe a data structure for each path length. Note that our “time budget” per Hungarian search is $O(k\mathop{\mathrm{polylog}}n)$.

Finding length-1 paths.

This data structure finds a min-reduced-cost arc from an alive node of $X$ to an alive node of $V\setminus X$. There are $O(k)$ backward arcs, so the minimum among backward arcs can be maintained explicitly in a priority queue and retrieved in $O(1)$ time.

There are three types of forward arcs: $s\mathord{\shortrightarrow}a$ for some $a\in{A}^{\text{\char 96}}$, $b\mathord{\shortrightarrow}t$ for some $b\in{B}^{\text{\char 96}}$, and bipartite arc $a\mathord{\shortrightarrow}b$ with two alive endpoints. Arcs of the first (resp. second) type can be found by maintaining ${A}^{\text{\char 96}}\setminus X$ (resp. ${B}^{\text{\char 96}}\cap X$) in a priority queue, but should only be queried if $s\in X$ (resp. $t\not\in X$). The cheapest arc of the third type can be maintained using a dynamic (additively weighted) bichromatic closest pair (BCP) data structure between ${A}^{\text{\char 96}}\cap X$ and ${B}^{\text{\char 96}}\setminus X$, with reduced cost as the weighted pair distance. Such BCP data structure can be implemented so that insertions/deletions can be performed in $O(\mathop{\mathrm{polylog}}k)$ time [11].

Finding length-2 paths.

We describe how to find a cheapest path of the form $s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}b$ where $v$ is dead and $b\in{B}^{\text{\char 96}}$. A cheapest path $a\mathord{\shortrightarrow}v\mathord{\shortrightarrow}t$ can be found similarly. Similar to length-1 paths, we only query paths starting at $s$ if $s\in X$ and paths ending at $t$ if $t\not\in X$.

Note that $c_{\pi}(s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}b)=c(v,b)+\pi(b)-\pi(s)$. Since $\pi(s)$ is common in all such paths, it suffices to find a pair $(v,w)$ between ${A}\setminus{A}^{\text{\char 96}}$ and ${B}^{\text{\char 96}}\setminus X$ minimizing $c(v,w)+\pi(w)$. This is done by maintaining a dynamic BCP data structure between ${A}\setminus{A}^{\text{\char 96}}$ and ${B}^{\text{\char 96}}\setminus X$ with the cost of a pair $(v,w)$ being $c(v,w)+\pi(w)$. We may require an update operation for each alive node added to $X$ during the Hungarian search, of which there are $O(k)$, so the time spent during a search is $O(k\mathop{\mathrm{polylog}}n)$.

Since the size of ${A}\setminus{A}^{\text{\char 96}}$ is at least $r-k$, we cannot construct this BCP from scratch at the beginning of each iteration. To resolve this, we use the idea of rewinding from Section 2, with a slight twist. There are now two ways that the initial BCP may change across consecutive Hungarian searches: (1) the initial set $X$ may change as nodes lose excess through augmentation, and (2) the set of alive/dead nodes in $A$ may change. The first is identical to the situation in Section 2; the number of excess depletions is $O(k)$ over the course of Refine. For the second, the alive/dead status of a node can change only if the blocking flow found passes through the node. By Lemma 3.7 below, there are $O(k)$ such changes per Hungarian search, which can be done in $O(k\mathop{\mathrm{polylog}}n)$ time.

Finding length-3 paths.

We now describe how to find the cheapest path of the form $s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}w\mathord{\shortrightarrow}t$ where $v\in{A}\setminus{A}^{\text{\char 96}}$ and $w\in{B}\setminus{B}^{\text{\char 96}}$. Note that $c_{\pi}(s\mathord{\shortrightarrow}v\mathord{\shortrightarrow}w\mathord{\shortrightarrow}t)=c(v\mathord{\shortrightarrow}w)-\pi(s)+\pi(t)$. A pair $(v,w)$ between ${A}\setminus{A}^{\text{\char 96}}$ and ${B}\setminus{B}^{\text{\char 96}}$ minimizing $c(v,w)$ can be found by maintaining a dynamic BCP data structure similar to the case of length-2 paths.

This BCP data structure has no dependency on $X$—the only update required comes from membership changes to ${A}^{\text{\char 96}}$ or ${B}^{\text{\char 96}}$ after an augmentation. Applying Lemma 3.7 again, there are $O(k)$ alive/dead updates caused by an augmentation, so the time for these updates per Hungarian search is $O(k\mathop{\mathrm{polylog}}n)$.

Updating potential.

Potential updates for alive nodes can be handled in a batched fashion as in Section 2. The three data structures above have no dependency on the dead node potential; we leave them untracked as described before. The Hungarian search remains intact since alive nodes are visited in the same order as when using arc-by-arc relaxations. However, we need values of $\pi$ on all nodes at the end of a scale (for the next Scale-Init) and for individual dead nodes whenever they become alive (after augmentation).

We can reconstruct a “valid” potential in these situations. To recover potential for $v\in{A}\setminus{A}^{\text{\char 96}}$ we set $\pi(v)\leftarrow\pi(s)$, and for $v\in{B}\setminus{B}^{\text{\char 96}}$ we set $\pi(v)\leftarrow\pi(t)$. Straightforward calculation shows that such potential (1) preserves $\theta$-optimality, and (2) makes $\Pi$ (arc-wise) admissible for any weakly admissible alive path $\Pi$. Hence, a blocking flow composed of weakly admissible alive paths is admissible under the recovered potential.

The following lemma is crucial to the analysis of running time for the Hungarian search, bounding both the number of relaxations and potential update/recovery operations.

Lemma 3.7**.**

Both Hungarian search and augmentation stages explore $O(k)$ nodes, and the blocking flow found in augmentation stage is incident to $O(k)$ nodes.

Augmentation can also be implemented in $O(k\mathop{\mathrm{polylog}}n)$ time, after $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, using similar data structures. We thus obtain the following:

Lemma 3.8**.**

After $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, each iteration of Refine can be implemented in $O(k\mathop{\mathrm{polylog}}n)$ time.

4 Transportation algorithm

Given two point sets $A$ and $B$ in $\mathbb{R}^{2}$ of sizes $r$ and $n$ respectively and a supply-demand function $\phi:A\cup B\to\mathbb{Z}$ as defined in the introduction, we present an $O(rn^{3/2}\mathop{\mathrm{polylog}}n)$ time algorithm for computing an optimal transport map between $A$ and $B$. By applying this algorithm in the case of $r\leq\sqrt{n}$ and the one by Agarwal et al. [3] when $r>\sqrt{n}$, we prove Theorem 1.3. We use a standard reduction to the uncapacitated min-cost flow problem and use Orlin’s algorithm [16] as well as some of the ideas from Agarwal et al. [3] for efficient implementation under the geometric settings. We first present an overview of the algorithm and then describe its fast implementation that achieves the desired running time.

4.1 Overview of the algorithm

Orlin’s algorithm follows an excess-scaling paradigm and the primal-dual framework. It maintains a scale parameter $\Delta$, a flow function $f$, and potential $\pi$ on the nodes. Initially $\Delta$ is equal to the total supply, $f=0$, and $\pi=0$. We fix a constant parameter $\alpha\in(0.5,1)$. A node $v$ is called active if the magnitude of imbalance of $v$ is at least $\alpha\Delta$. At each step, using the Hungarian search, the algorithm finds an admissible excess-to-deficit path between active nodes in the residual network and pushes a flow of amount $\Delta$ along this path.333Note that this augmentation may convert an excess node into a deficit node. Repeat the process until either active excess or deficit nodes are gone; when this happens, $\Delta$ is halved. The sequence of augmentations with a fixed value of $\Delta$ is called an excess scale.

The algorithm also performs two preprocessing steps at the beginning of each excess scale. First, if $f(v\mathord{\shortrightarrow}w)\geq 3n\Delta$, $v\mathord{\shortrightarrow}w$ is contracted to a single node $z$ with $\phi(z)=\phi(v)+\phi(w)$.444Intuitively, $f(v\mathord{\shortrightarrow}w)$ is so high that future scales cannot deplete the flow on $v\mathord{\shortrightarrow}w$. Second, if there are no active excess nodes and $f(v\mathord{\shortrightarrow}w)=0$ for every arc $v\mathord{\shortrightarrow}w$, then $\Delta$ is aggresively lowered to $\max_{v}\phi(v)$.

When the algorithm terminates, an optimal circulation in the contracted network is found. We use the algorithm described in Agarwal et al. [3] to recover an optimal circulation for the original network in $O(n\mathop{\mathrm{polylog}}n)$ time. Orlin showed that the algorithm terminates within $O(n\log n)$ scales and performs a total of $O(n\log n)$ augmentations. In the next subsection, we describe an algorithm that, after $O(n\mathop{\mathrm{polylog}}n)$ time preprocessing, finds an admissible excess-to-deficit path in $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ amortized time. Summing this cost over all augmentations, we obtain the desired running time.

4.2 An efficient implementation

Recall in the previous sections that we could bound the running time of the Hungarian search by the size of $\operatorname{supp}(f)$. Here, the number of active imbalanced nodes at any scale is $O(r)$, and the length of an augmenting path is also $O(r)$. Therefore one might hope to find an augmenting path in $O(r\mathop{\mathrm{polylog}}n)$ time, by adapting the algorithms described in Sections 2 and 3. The challenge is that $\operatorname{supp}(f)$ may have $\Omega(n)$ size, therefore an algorithm which runs in time proportional to the support size is no longer sufficient. Still, we manage to implement Hungarian search in time $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$, by exploiting a few properties of $\operatorname{supp}(f)$ as described below.

We note that each arc of $\operatorname{supp}(f)$ is admissible with reduced cost [math], so we prioritize the relaxation of support arcs as soon as they arrive in $X\times(V\setminus X)$, over the non-support arcs. This strategy ensures the following crucial property.

Lemma 4.1**.**

If the support arcs $\operatorname{supp}(f)$ are relaxed as soon as possible, $\operatorname{supp}(f)$ is acyclic.

Next, similar to Section 3, we call node $u$ alive if (a) $u$ is an active imbalanced node or (b) if $u$ is incident to an arc of $\operatorname{supp}(f)$; $u$ is dead otherwise. Unlike in Section 3, once a node becomes alive it cannot be dead again. Furthermore, a dead node may become alive only at the beginning of a scale (after the value of $\Delta$ is reduced). Also, an augmenting path cannot pass through a dead node. Therefore, we can ignore all dead nodes during Hungarian search, and update the set of alive/dead nodes at the beginning of a scale.

Let {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{B^{\text{\char 84}}}}}\subseteq{B}^{\text{\char 96}} be the set of nodes that are either (a) active imbalanced nodes or (b) incident to exactly one arc of $\operatorname{supp}(f)$. Lemma 4.1 implies that ${B}^{\text{\char 96}}\setminus B^{\text{\char 84}}$ has size $O(r)$. We can therefore find the min-reduced-cost arc between $X\cap{A}^{\text{\char 96}}$ and ${B}^{\text{\char 96}}\setminus(B^{\text{\char 84}}\cup X)$ using a BCP data structure as in Section 2, along with lazy potential updates and the rewinding mechanism. The total time spent by Hungarian search on the nodes of ${B}^{\text{\char 96}}\setminus B^{\text{\char 84}}$ will be $O(r\mathop{\mathrm{polylog}}n)$. We subsequently focus on handling $B^{\text{\char 84}}$.

Handling $B^{\text{\char 84}}$.

We now describe how we query a min-reduced-cost arc between $X\cap{A}^{\text{\char 96}}$ and $B^{\text{\char 84}}\setminus X$. Each node $b\in B^{\text{\char 84}}$ is incident to exactly one arc in $\operatorname{supp}(f)$. We partition these nodes into clusters depending on their unique neighbor in $N_{f}$. That is, for a node $a\in{A}^{\text{\char 96}}$, let {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{B_{a}^{\text{\char 84}}}}}\coloneqq\{b\in B^{\text{\char 84}}\mid{a\mathord{\shortrightarrow}b}\in\operatorname{supp}(f)\}. We refer to $B_{a}^{\text{\char 84}}$ as the star of $a$.

The crucial observation is that $a$ is the only node in $N_{f}$ reachable from each $b\in B_{a}^{\text{\char 84}}$, so once the Hungarian search reaches a node of $B_{a}^{\text{\char 84}}$ and thus $a$ (recall we prioritize relaxing support arcs), the Hungarian search need not visit any other nodes of $B_{a}^{\text{\char 84}}$, as they will only lead to $a$. Hence, as soon as one node of $B_{a}^{\text{\char 84}}$ is reached, all other nodes of $B_{a}^{\text{\char 84}}$ can be discarded from further consideration. Using this observation, we handle $B^{\text{\char 84}}$ as follows.

We classify each $a\in{A}^{\text{\char 96}}$ as light or heavy: heavy if $\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}\geq\sqrt{n}$, and light if $\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}\leq 2\sqrt{n}$. Note that if $\sqrt{n}\leq\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}\leq 2\sqrt{n}$ then $a$ may be classified as light or heavy. We allow this flexibility to implement reclassification in a lazy manner. Namely, a light node is reclassified as heavy once $\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}>2\sqrt{n}$, and a heavy node is reclassified as light once $\mathopen{|}B_{a}^{\text{\char 84}}\mathclose{|}<\sqrt{n}$. This scheme ensures that the star of $a$ has gone through at least $\sqrt{n}$ updates between two successive reclassifications, and these updates will pay for the time spent in updating the data structure when $a$ is re-classified.

For each heavy node $a\in{A}^{\text{\char 96}}\setminus X$, we maintain a BCP data structure between $B_{a}^{\text{\char 84}}$ and $X\cap{A}^{\text{\char 96}}$. Next, for all light nodes in ${A}^{\text{\char 96}}\setminus X$, we collect their stars into a single set {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}{\emph{B_{<}^{\text{\char 84}}}}}\coloneqq\bigcup_{\text{a light}}B_{a}^{\text{\char 84}}. We maintain one single BCP data structure between $B_{<}^{\text{\char 84}}$ and ${A}^{\text{\char 96}}\cap X$. Thus, at most $r$ different BCP data structures are maintained for stars.

Using these data structures, we can compute and relax a min-reduced-cost arc $v\mathord{\shortrightarrow}w$ between ${A}^{\text{\char 96}}\cap X$ and $B^{\text{\char 84}}\setminus X$. If $w$ lies in some star $B_{a}^{\text{\char 84}}$, then we also add $a$ into $X$. If $a$ is light, then we delete $B_{a}^{\text{\char 84}}$ from $B_{<}^{\text{\char 84}}$ and update the BCP data structure of $B_{<}^{\text{\char 84}}$. If $a$ is heavy, then we stop querying the BCP data structure of $B_{a}^{\text{\char 84}}$ for the remainder of the search. Finally, since $a$ becomes part of $X$, $a$ is added to all $O(r)$ BCP data structures. Recall that $r\leq\sqrt{n}$ by assumption. Adding arc $v\mathord{\shortrightarrow}w$ thus involves performing $O(\sqrt{n})$ insertion/deletion operations in various BCP data structures, thereby taking $O(\sqrt{n}\mathop{\mathrm{polylog}}n)$ time.

Putting it together.

While proof is omitted, the following lemma bounds the running time of the Hungarian search.

Lemma 4.2**.**

Assuming all BCP data structures are initialized correctly, the Hungarian search terminates within $O(r)$ steps, and takes $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ time.

Once an augmenting path is found and the augmentation is performed, the set of imbalanced nodes and the support arcs change. We thus need to update the sets $B^{\text{\char 84}}$, $B_{a}^{\text{\char 84}}$s, and $B_{<}^{\text{\char 84}}$. This can be accomplished in $O(r\mathop{\mathrm{polylog}}n)$ amortized time. When we begin a new Hungarian search, we use the rewinding mechanism to set various BCP data structures in the right initial state. Finally, when we move from one scale to another, we also update the sets ${A}^{\text{\char 96}}$ and ${B}^{\text{\char 96}}$. Omitting all the details, we conclude the following.

Lemma 4.3**.**

Each Hungarian search can be performed in $O(r\sqrt{n}\mathop{\mathrm{polylog}}n)$ time.

Since there are $O(n\log n)$ augmentations and the flow in the original network can be recovered from that in the contracted network in $O(n\mathop{\mathrm{polylog}}n)$ time [3], the total running time of the algorithm is $O(rn^{3/2}\mathop{\mathrm{polylog}}n)$, as claimed in Theorem 1.3.

Acknowledgment.

We thank Haim Kaplan for discussion and suggestion to use Goldberg et al. [9] algorithm. We thank Debmalya Panigrahi for helpful discussions.