Faster Algorithms for Reverse Shortest Path in Unit-Disk Graphs and Related Geometric Optimization Problems:
Improving the Shrink-And-Bifurcate Technique

As noted, the current fastest algorithms for all of the above problems are obtained via the shrink-and-bifurcate technique introduced by Avraham et al. [6]. The setting requires that critical values (values that the decision algorithm may compare the input value against) are distances defined by pairs of input objects for some distance function. The strategy consists of two stages, interval shrinking and bifurcation.

1. 1.

   The interval shrinking stage is a relaxed version of binary search – instead of insisting on finding the answer $r^{\ast }$ immediately, we seek an interval $(r^{-},r^{+}]$ that contains $r^{\ast }$ and that contains at most $L$ critical values for a given parameter $L$. In the case of Euclidean distances of points in the plane, Avraham et al. showed that this stage can be done in $O^{\ast }(n^{4/3}/L^{1/3}+D(n))$ time (using randomization), where $D(n)$ is the running time of the decision algorithm. This is potentially better than the $\tilde{O}(n^{4/3})$ time bound for distance selection.

2. 2.

   The bifurcation stage then uses this interval $(r^{-},r^{+}]$ to simulate the decision algorithm at the unknown $r^{\ast }$ more efficiently than standard parametric search when the decision algorithm is not parallelizable. Avraham et al.  showed that this simulation can be done in $\tilde{O}(L^{1/2}D(n))$ time.

When $D(n)=\tilde{O}(n)$, choosing $L=n^{2/5}$ to balance the costs of the two stages gives an overall running time of $O^{\ast }(n^{6/5})$.

Some of the applications use distance functions other than Euclidean distances of points in the plane, but the $O^{\ast }(n^{4/3}/L^{1/3}+D(n))$ bound for interval shrinking is still applicable. An exception is the reverse shortest path problem for general disk graphs. Here, the corresponding distance selection problem requires $O^{\ast }(n^{3/2})$ time, and the interval shrinking subproblem is solved in $O^{\ast }(n^{3/2}/L^{1/2}+D(n))$ time instead. When $D(n)=\tilde{O}(n)$, choosing $L=n^{1/2}$ to balance the costs of the two stages gives an overall running time of $O^{\ast }(n^{5/4})$ instead.

The new contributions in this paper are twofold:

1. 1.

   We improve all the previous $O^{\ast }(n^{6/5})$ time bounds – for reverse shortest path in unit-disk graphs, weighted unit-disk graphs, segment proximity graphs, visibility graphs over 1.5-dimensional terrains, and discrete Fréchet distance with one-sided shortcuts – to $\tilde{O}(n^{8/7})$. (In particular, for discrete Fréchet distance with one-sided shortcuts, this is the first improvement on the exponent in 10 years.)

   For reverse shortest path in arbitrary disk graphs (unweighted or weighted), we improve the previous $O^{\ast }(n^{5/4})$ time bound to $O^{\ast }(n^{6/5})$.

   Although the improvements are incremental, they are all obtained from one simple idea for improving interval shrinking (we keep the bifurcation part unchanged) – the simplicity and generality of our idea make it easy to apply to existing (and potentially future) applications of the shrink-and-bifurcate approach, and are the main takeaways of the paper. More precisely, we show how to lower the $O^{\ast }(n^{4/3}/L^{1/3}+D(n))$ time bound for interval shrinking subproblem to $\tilde{O}(n^{4/3}/L^{2/3}+D(n))$ (or in the arbitrary disk graph case, lower $O^{\ast }(n^{3/2}/L^{1/2}+D(n))$ to $O^{\ast }(n^{3/2}/L^{3/4}+D(n))$). The original interval shrinking method by Avraham et al. needs to modify range searching techniques for distance selection, constructing “partial biclique covers” and then applying random sampling. Our new method, in contrast, can be described in just a few lines, and simply uses a known distance selection algorithm as a black box on random subsets (see Section 2).

2. 2.

   Next, we focus more closely on the simplest of this class of problems, namely, reverse shortest path in (unweighted) unit-disk graphs. For this specific problem, we are able to further improve the running time to $\tilde{O}(n^{9/8})$. The improvement of the exponent from $6/5=1.2$ to $9/8=1.125$ is a bigger increment (for example, comparable in magnitude to Wang and Zhao’s original improvement [29] from $4/3\approx 1.333$ to $5/4=1.25$, or Kaplan et al.’s improvement [21] from $5/4=1.25$ to $6/5=1.2$).

   Interestingly, we obtain our result by combining the shrink-and-bifurcate technique with Wang and Zhao’s original grid-based approach! Although it may be natural to consider combining the two approaches, how to do so is not obvious at all. We need to put together several nontrivial technical ideas – this is the most intricate part of the paper (see Section 3).

Our first interval shrinking algorithm is simple, as presented in Algorithm 1.

Upon closer examination, one may notice that the basic algorithm has room for improvement. In particular, if $E$ has $\mathrm{\Theta }(\sqrt{L})$ non-isolated vertices in $P$ and $\mathrm{\Theta }(\sqrt{L})$ non-isolated vertices in $Q$, then clearly our samples $\widehat{P}\times Q$ and $P\times \widehat{Q}$ are wasteful. Indeed, we can eliminate such redundancy by sampling in $P$ and $Q$ simultaneously with different rates. We present the details in Algorithm 2.

After interval shrinking, Avraham et al. [6] devised a clever variant of parametric search via an idea called bifurcation (also redescribed in more generality in Kaplan et al.’s subsequent paper [21]) which accomplishes the following:

Typically, in applications of the technique, the algorithm $𝒜$ we want to simulate is the decision algorithm itself; here, the simulation produces an interval containing $r^{\ast }$ such that all input values lying in the interval have the same computation path and thus the same answer as $r^{\ast }$ (namely, “yes”), and so the left endpoint of this interval gives us precisely the value of $r^{\ast }$.

We will not redescribe the proof of the above lemma, as we will only use it as a black box. (The rough idea in bifurcation is to simulate multiple branches of $𝒜$ at the same time, until we have accumulated sufficiently many comparisons and will then resolve all of them as a batch using the decision oracle; comparisons with values outside the interval $(r^{-},r^{+}]$ are free.) Combining our improved interval shrinking method with bifurcation, we immediately obtain better results for all optimization problems stated in the introduction:

Our algorithms rely on the details of the known $\tilde{O}(n)$-time decision algorithm by Chan and Skrepetos [12], so we will first give a quick recap of their algorithm.

Given an input value $r$, we want to compute $d_r(s,t)$. We first create a uniform grid ${\mathrm{\Psi }}_r$ of side length $r/\sqrt{2}$. Chan and Skrepetos’ algorithm is a careful implementation of BFS using the grid. We start from $s$ and generate each level of the BFS tree (the “frontier”) one at a time. The key observation is that each grid cell only needs to be visited a constant number of times during the BFS, because (i) points in the same cell have the roughly same shortest-path distance $\pm 1$ from $s$ in $G_r(P)$ and so each cell “appears” in only $O(1)$ levels of the BFS tree, and (ii) each cell has only $O(1)$ neighboring cells – we say that two cells are neighboring if the minimum Euclidean distance between the two cells is at most $r$. Given one level of the BFS tree, we can generate the next level (i.e., “expand” the frontier) by solving the following subproblems between cells in the current level and their neighboring cells:

Chan and Skrepetos solved the above subproblem in $O(n_r+n_b)$ time, assuming that the input points are pre-sorted, by using envelopes of pseudo-lines. As a result, the total running time of their algorithm is linear after pre-sorting.

The above algorithm is not efficiently parallelizable (since the number of levels in the BFS may be large). Instead, we can apply the bifurcation technique to simulate the algorithm on $r^{\ast }$. However, one issue arises: the envelope computations requires critical values defined by triples of input elements, not pairs. As noted in other papers [21, 29], this can be fixed by switching to a different solution to Subproblem 5 which is slightly slower by a logarithmic factor but is simpler: namely, we compute the nearest red neighbor to each blue point by point location in the Voronoi diagram of the red points in $O((n_r+n_b)\log (n_r+n_b))$ time [15]; then for each blue point, we just compare $r$ with its nearest neighbor distance. This way, the critical values are all Euclidean distances of pairs of input points.

Another issue is the construction of the grid ${\mathrm{\Psi }}_r$, i.e., the assignment of input points to grid cells in ${\mathrm{\Psi }}_r$. This part is easily parallelizable, and so we can apply standard parametric search to construct the grid ${\mathrm{\Psi }}_{r^{\ast }}$ in $\tilde{O}(n)$ time, all done as preprocessing. (In fact, the time bound for constructing ${\mathrm{\Psi }}_{r^{\ast }}$ is $O(n\log n)$, as shown by Wang and Zhao [29].)

Applying the shrink-and-bifurcate technique to the preceding decision algorithm using our improved interval shrinking method yields an overall time bound of $\tilde{O}(n^{8/7})$, as we have already shown in Theorem 4(a). To obtain our best result, we will combine with another idea that was used in Wang and Zhao’s original paper [29]: heavy vs. light cells.

Recall that we have already constructed the grid ${\mathrm{\Psi }}_{r^{\ast }}$. Let $\mathrm{\Delta }$ be a parameter to be chosen later. Call a grid cell $\sigma$ heavy if $|P\cap \sigma |\ge \mathrm{\Delta }$, or light otherwise. The number of heavy cells is clearly at most $O(n/\mathrm{\Delta })$. Furthermore, a pair of neighboring cells $(\sigma ,{\sigma }^{\prime })$ is called a light pair if one of $\sigma$ and ${\sigma }^{\prime }$ is a light cell (the other may be light or heavy).

Our new strategy is as follows. First, we perform interval shrinking, but only with respect to distances of points from light pairs. Intuitively, there are fewer such distances and so interval shrinking should be doable faster. This allows us to simulate Chan and Skrepetos’ algorithm at $r^{\ast }$ on a modified graph where heavy cells are eliminated – or more precisely, points in a common heavy cell are contracted into a single vertex. When we “uncontract” the heavy cells, this causes an $O(n/\mathrm{\Delta })$ additive error to all the shortest-path distances in $G_{r^{\ast }}(P)$. As a final step, we show how to use a dynamic program to recover the value $r^{\ast }$. In the next three subsections, we provide the details of all these steps.

In the following, we show how to solve the interval shrinking problem faster for distances from light pairs. More precisely, define ${\delta }_{\text{light}}(p,q)=\delta (p,q)$ if $p\in \sigma$ and $q\in {\sigma }^{\prime }$ for some light pair $(\sigma ,{\sigma }^{\prime })$, and ${\delta }_{\text{light}}(p,q)=\mathrm{\infty }$ otherwise. We apply Algorithm 2 to do interval shrinking with respect to ${\delta }_{\text{light}}$.

To analyze the running time, we need an efficient implementation of the selection oracle for

Let $n_{\sigma }:=|P\cap \sigma |$ for each cell $\sigma$. Note that $|{\widehat{P}}_i\cap \sigma |=\tilde{O}(\frac{n_{\sigma }}{\lceil L/2^i\rceil })$ and $|{\widehat{Q}}_i\cap {\sigma }^{\prime }|=\tilde{O}(\frac{n_{{\sigma }^{\prime }}}{2^i})$ for every $(\sigma ,{\sigma }^{\prime })$ w.h.p. We already have a selection oracle for each $\delta ({\widehat{P}}_i\cap \sigma ,{\widehat{Q}}_i\cap {\sigma }^{\prime })$ that w.h.p. runs in time

To obtain a selection oracle for ${\delta }_{\text{light}}({\widehat{P}}_i,{\widehat{Q}}_i)$, we face a problem analogous to selection in a union of multiple sorted arrays, except that the arrays are implicitly given. Frederickson and Johnson [18] gave an algorithm for selection in $k$ sorted arrays of size $n$ running in $O(k\log n)$ time; the algorithm proceeds in $O(\log n)$ rounds, and in each round, it looks up only $O(1)$ elements in each array. In our scenario, each array lookup corresponds to a call to the selection oracle for $\delta ({\widehat{P}}_i\cap \sigma ,{\widehat{Q}}_i\cap {\sigma }^{\prime })$ for a light pair $(\sigma ,{\sigma }^{\prime })$. Thus, by Frederickson and Johnson’s algorithm, selection in the union ${\delta }_{\text{light}}({\widehat{P}}_i,{\widehat{Q}}_i)$ can be done in total time

since $\min \{n_{\sigma },n_{{\sigma }^{\prime }}\}\le \mathrm{\Delta }$ for each light pair $(\sigma ,{\sigma }^{\prime })$, and ${\sum }_{\text{neighboring pair }(\sigma ,{\sigma }^{\prime })}(n_{\sigma }+n_{{\sigma }^{\prime }})=O(n)$. Therefore, the total time for interval shrinking with respect to ${\delta }_{\text{light}}$ is $\tilde{O}({\mathrm{\Delta }}^{1/3}n/L^{2/3}+D(n))$ w.h.p.

Let $H_r$ be the graph obtained from $G_r(P)$ by contracting the points inside each heavy cell of ${\mathrm{\Psi }}_r$ into a single vertex. Since there are $O(n/\mathrm{\Delta })$ contracted vertices in $H_r$, we have $d_r(s,p)-k\le d_{H_r}(s,p)\le d_r(s,p)$ with $k=O(n/\mathrm{\Delta })$, for every $p\in P$.

It is straightforward to modify Chan and Skrepetos’ decision algorithm to do BFS on this modified graph $H_r$. Expanding the frontier again reduces to solving Subproblem 5 for pairs of neighboring cells. Note that for two neighboring cells that are both heavy, we can determine whether there is an edge between them by pre-computing the bichromatic closest pair between the points in the neighboring cells (via Voronoi diagrams) and testing $r$ against this closest pair distance – call this a special distance. This pre-computation takes $O(n\log n)$ total time.

After performing interval shrinking as described in previous subsection, we apply the bifurcation technique (Lemma 3) to simulate BFS on $H_{r^{\ast }}$. The critical values indeed are all Euclidean distances from light pairs only, except for the $O(n)$ special distances above, but initially before the simulation starts, we can resolve the comparisons with all $O(n)$ special distances by binary search using $O(\log n)$ calls to the decision algorithm in $\tilde{O}(n)$ time.

At the end of the simulation, we will then have computed $d_{H_{r^{\ast }}}(s,p)$ for all $p\in P$, which yield estimates to $d_{r^{\ast }}(s,p)$ with additive error $O(n/\mathrm{\Delta })$. To summarize:

Finally, we will compute $r^{\ast }$ using the estimates obtained in Lemma 6. This final step is done via a dynamic program, and does not require any knowledge of the details of the earlier steps.

To prove Lemma 7, it is helpful to view the problem as computing a minimum bottleneck path from $s$ to $t$ with at most $\lambda$ links (as mentioned in the introduction), i.e., a path from $s$ to $t$ with at most $\lambda$ links that minimizes the maximum Euclidean length of the edges along the path.

There is a straightforward dynamic programming solution to this problem: Let $D_i(p)$ denote the minimum bottleneck value over all paths from $s$ to $p$ with $i$ links. We want $r^{\ast }={\min }_{i\le \lambda }{D}_i(t)$. We have the following recursive formula: for each $p\in P$,

with $D_0(s)=0$ and $D_0(p)=\mathrm{\infty }$ for all $p\in P-\{s\}$ as base cases. Unfortunately, in this dynamic program, the number of table entries is $O(n\lambda )$, which is $O(n^2)$ in the worst case.

We use the given estimates ${\tilde{d}}_p$ to speed up the dynamic program. The main observation is that it suffices to generate $D_i(p)$ when $i-k\le {\tilde{d}}_p\le i$. The number of such table entries is thus reduced to $O(kn)$.

More precisely, define a subset $P_i=\{p\in P:i-k\le {\tilde{d}}_p\le i\}$ for each $i$. Since a point belongs to $O(k)$ of these subsets, ${\sum }_i|P_i|=O(kn)$. For each $p\in P_i$, define

with ${\widehat{D}}_0(s)=0$. (All unspecified ${\widehat{D}}_i(p)$ values are implicitly set to $\mathrm{\infty }$.)

We claim that $r^{\ast }={\min }_{i\le \lambda }{\widehat{D}}_i(t)$.