arXiv:1601.02867v1 [cs.DS] 12 Jan 2016

Department of Computer Science, TU Dortmund, Germany {denis.kurz,petra.mutzel}@tu-dortmund.de

Abstract. We present an algorithm for the k shortest simple path problem on weighted directed graphs (kSSP) that is based on Eppstein’s algorithm for a similar problem in which paths are allowed to contain cycles. In contrast to most other algorithms for kSSP, ours is not based on Yen’s algorithm [14] and does not solve replacement path problems. Its worst-case running time is on par with state-of-the-art algorithms for kSSP. Using our algorithm, one may find O(m) simple paths with a single shortest path tree computation and O(n + m) additional time per path in well-behaved cases, where n is the number of nodes and m is the number of edges. Our computational results show that on random graphs and large road networks, these well-behaved cases are quite common and our algorithm is faster than existing algorithms by an order of magnitude. Further, the running time is far better predictable due to very small dispersion. Keywords: directed graph, k-best, shortest path, simple path, weighted graph

1

Introduction

The k shortest path problem in weighted, directed graphs (kSP) asks for a set of k paths from a source s to a target target t in a graph with n nodes and m edges. Every path that is not output by an algorithm should be at least as long as any path in the output. Algorithms for this problem can be useful tools when it is hard to specify constraints that a solution should satisfy. Instead of computing only one shortest path, kSP algorithms generate k paths, and the user can then pick the one that suits their needs best. The best known algorithm for this problem runs in time O(m + n log n + k log k) and is due to Eppstein [3]. In the initialization phase, the algorithm builds a data structure that contains information about all s-t paths and how they interrelate with each other, in time O(m + n log n). This can even be reduced to O(m + n) if the shortest path tree (SP tree) is given in the input or if the SP tree can be computed in time O(m+n). In the enumeration phase, a path graph is constructed. The path graph is a quaternary min-heap where every path starting in the root correlates to an s-t path in the original graph. We require O(k log k) time for the enumeration phase if we want the output paths to be ordered by length. If the order in which

the paths are output does not matter, Frederickson’s heap selection algorithm [6] can be used to enumerate the paths after the initialization phase in time O(k). The k shortest simple path problem (kSSP), introduced by Clarke, Krikorian and Schwartz [2], seems to be more expensive, computationally. In contrast to kSP, the computed paths are required to be simple, i.e., they must not contain a cycle. The extra effort may be well-invested if many of the k shortest paths are non-simple and we are only interested in simple paths. The algorithm by Yen [14] used to have the best theoretical worst-case running time of O(kn(m + n log n)) for quite some time. Gotthilf and Lewenstein [7] improved upon this bound recently. They observed that kSSP can be solved by solving O(k) all pairs shortest path (APSP) instances. Using the APSP algorithm by Pettie [11], they obtain a new upper bound of O(kn(m + n log log n)). Vassilevska Williams and Williams [13] showed that, for constant k, an algorithm for kSSP with running time O(n3−ε ) for some positive ε (truly subcubic) would also yield algorithms with truly subcubic running times for some other problems, including APSP. A recent survey of the field is due to Eppstein [4]. The kSSP on undirected graphs seems to be significantly easier. Katoh et al. [9] proposed an algorithm that solves kSSP on undirected graphs in time O(k(m + n log n)). A subproblem occurring in Yen’s algorithm is the (restricted) replacement path problem. Given a shortest s-t path p in a graph, it asks for a set of paths as follows. For each i < |p|, the set has to include a shortest simple path that uses the first i−1 edges of p, but not the ith. This problem has to be solved O(k) times to find the k shortest simple paths using Yen’s algorithm. In the original version of Yen’s algorithm, the replacement paths are found using O(|p|) shortest path computations, resulting in time O(n(m+n log n)). Hershberger et al. [8] compute one SP tree rooted in s and one reversed SP tree rooted in t, respectively. They use these two trees to find a replacement path in constant time per edge on p, cutting down the time required to find all replacement paths to O(m + n log n) when Dijkstra’s algorithm is used. However, the paths generated this way are not guaranteed to be simple. Such non-simple paths can be detected in constant time and repaired by falling back to Yen’s replacement path computation for the path edge in question. Since they do not provide an upper bound for the number of non-simple paths that may occur using this method, the worst-case running time is again O(n(m + n log n)). Some approaches reuse one fixed reversed SP tree T0 rooted in t and computed during the initialization of their kSSP algorithm, in contrast to O(1) SP trees per replacement path instance. Pascoal [10] noticed that the replacement path that deviates from p at node v might be one that uses an edge (v, w) to an unused successor w of v and then follows the path from w to t in T0 . Therefore, they test whether the shortest such path is simple, and fall back to a full shortest path computation if it is not. Although they do not describe in detail how this check is done, it can be done in time O(m + n) per replacement path instance by partitioning the nodes into blocks as described by Hershberger et al. [8]. Feng [5] uses the reversed SP tree to partition V into three classes. For each edge (u, v) on p for which we want to compute a replacement path, red nodes have already 2

been used to reach v via p. A yellow node v is a non-red upstream node of some red node in T0 , i.e., the path from v to t in T0 contains a red node. All other nodes are green. They then do shortest path computations from v using Dijkstra’s algorithm like Yen. However, they are able to restrict the search to yellow nodes, resulting in a significantly smaller search space. Feng does not provide upper bounds on the size of this search space, resulting again in a worst-case running time of O(n(m + n log n)) for each replacement path instance. Our contribution. We propose an algorithm that was derived from Eppstein’s notion of a path graph [3]. Our algorithm achieves the same worst-case running time as Yen’s algorithm. Like Yen, we rely on shortest path (tree) computations. In contrast to Yen-based algorithms, however, our algorithm may draw O(m) simple paths from one shortest path tree computation. If the underlying graph is acyclic, the revised algorithm at the end of this paper requires O(n log n + k(m + n)) without further modifications. Alternatively, one could test whether the graph is acyclic and then use Eppstein’s algorithm. However, this method fails if the graph has just a single cycle, in which case our algorithm appears to be a good choice. Our algorithm works on multigraphs without modification. This is also true for every other kSSP algorithm we know of. After some definitions in Section 2, we propose a simplified version of our algorithm with running time O(km(m + n log n)) in Section 3. In Section 4, we show how this running time can be reduced to O(kn(m + n log n)), and how to reduce the number of shortest path tree computations in practice. Finally, we present the results of our computational studies in Section 5 to prove the efficiency of our algorithm.

2

Definitions

Let G = (V, E) be a directed graph with node set V and edge set E. Let s, t ∈ V be source and target nodes, respectively. We assume an implicit edge weight function c : E → R+ 0 throughout this paper. We denote the number of nodes |V | by n and the number of edges |E| by m. A path connecting v to w in G, or v-w path, is an edge sequence p = (e1 , e2 , . . . , er ), ei = (vi , wi ), with v = v1 , w = wr and wi = vi+1 for 1 ≤ i < r. For the sake of simplicity, we only consider combinations of G, s and t such that there exists an s-v path and a v-t path in G for every v ∈ V . A node u is said to be on the path p, denoted by u ∈ p, if u = w or u = vi for some i. If vi 6= vj 6= w for 1 ≤ i, j ≤ r, p is a simple path. The prefix (e1 , . . . , ei ) is a v-wi path and denoted by pi . The length c(p) of the path p is the sum of edge weights of its edges. If every v-w path is at least as long as p, it is called a shortest v-w path. We write G − p to denote the induced subgraph G[{v ∈ V | v ∈ / p}]. The k shortest simple path problem (kSSP) is an enumeration problem. Given a directed graph G = (V, E) with source node s ∈ V , target node t ∈ V , edge weights c, and some k ∈ N, we want to compute a set P comprising k simple paths from s to t in G such that c(p) ≤ c(p0 ) for every pair p ∈ P , p0 ∈ / P of simple paths. We obtain the k shortest path problem (kSP) if we do not require the computed paths to be simple. 3

A shortest path tree (SP tree) T of G is a subtree of G with node set V such that each v ∈ V has exactly one outgoing edge, which lies on a shortest v-t path, or no outgoing edges if no v-t path in G exists. We denote the latter case by v ∈ / T even if v is in the node set of T . Our algorithm will compute several SP trees, the first of which we call initial SP tree T0 . An edge e ∈ / T is called sidetrack w.r.t. T ; we will omit T in most cases. For a sidetrack e = (v, w), the sidetrack cost δT (e) is defined as (c(e) + d(w)) − d(v), where d(u) is the length of the unique u-t path in T . The sidetrack cost is therefore the difference between the length of a shortest v-t path and the length of a shortest v-t path that starts with e. The sidetrack set DT (v) of a node v ∈ V is the set of all sidetracks w.r.t. T with tails on the unique v-t path in T . When sidetracks are organized in heaps, we use sidetrack costs for comparison. Let p = (e1 , . . . , ek ), p0 = (f1 , . . . , fl ) be two s-t paths, and i∗ = max{i | ej = fj for 1 ≤ j < i}. Then, with respect to p, i∗ is the deviation index, the tail of ei∗ is the deviation node dev(p0 ), and ei∗ is the deviation edge of p0 . As is quite usual for kSSP algorithms, we will discover paths in a hierarchical fashion: a path p0 is added to the candidate set after p was extracted from the candidate set. In such cases, p is called the parent path of p0 . When p is omitted, the terms deviation node and edge are w.r.t. the parent path of p0 . By removing the deviation edge ∗ of p from p, p is split into its prefix path pref(p) := pi starting in s, and its suffix path suff(p) ending in t. The initial s-t path p0 in T0 has no parent path and thus no deviation edge. We define its suffix path to be p0 itself. We introduce a generalization of Eppstein’s representation [3] for paths. Eppstein represented paths as sequences of sidetracks, which were all sidetracks w.r.t. the same shortest path tree. In our representation, every sidetrack e in a sidetrack sequence may be associated with a different shortest path tree Te . The path represented by a sidetrack sequence (e1 , . . . , er ) can then be reconstructed as follows. Starting in s, we follow the initial SP tree T0 until we reach the tail of e1 . After reaching the tail of ei , we traverse ei and follow Tei until we reach the tail of ei+1 , or, in case i = r, until we reach t. Note that Eppstein’s representation is the special case where Te = T0 for each e in a sidetrack sequence, and that both Eppstein’s sidetrack sequences and our generalized ones may represent non-simple paths. The distance from a node v to t in a shortest path tree Te associated with a sidetrack e is denoted by de (v).

3

Basic Algorithm

In this section, we propose a rather simple way to enumerate the k shortest simple paths. We will describe later how to modify this algorithm to achieve our proclaimed running time guarantee. We initialize an empty priority queue Q that is going to manage candidate paths. The key of a path in Q is its length. We compute the initial shortest path tree T0 and push its unique s-t path, represented by an empty sidetrack sequence, to Q. We now process the paths in Q in order of increasing length until we found k simple paths. 4

() v1

v3 (a)

(c)

c s

t (a, b) a

v2

v4

(a, c)

(c0 )

d

b

(a, b, c)

(a) Example graph

(a, b, d)

(b) Sidetrack sequences

Fig. 1: Example for the basic algorithm. In Figure 1(a), the thick, solid edges belong to T0 . In Figure 1(b), every sidetrack is associated with T0 except for c0 , which is associated with the SP tree T1 comprising the edges b and d. An arrow from sequence p to sequence p0 indicates that p is the parent path of p0 . Let (e1 , . . . , er ) be a sidetrack sequence extracted from Q, and p the path that is represented by this sequence. Although the first path that is pushed to Q is always simple, we will eventually push non-simple paths to Q, too. Therefore, we first have to determine whether p is simple in a pivot step. This check can be done by simply walking p and marking every visited node. We first describe how to handle the simple case. We start by outputting p. For every sidetrack e = (u, v) with u ∈ suff(p), we discover a new path p0 represented by the sequence (e1 , . . . , er , e). We set dev(p0 ) = u, and push p0 to Q. The length of p0 can easily be computed as c(p) + δTe (e). If der (v) is undefined because Ter does not contain a v-t path, we simply ignore e. By choosing Te = Ter , we simply reuse the shortest path tree that is also associated with the last sidetrack in the sequence representing p. Note that sidetracks emanating from t can safely be ignored. Consider the example in Figure 1. The sidetrack sequence (a) with Ta = T0 represents a simple path p that passes the nodes s, v2 , v1 , v3 , t in this order. The suffix of this path is its v2 -t sub-path, and the sidetracks b and c have tails on this suffix. Therefore, when (a) is extracted from Q, p is output and the sequences (a, b) and (a, c) with Ta = Tb = Tc are pushed to Q. Now assume we extracted a non-simple path p represented by the sidetrack sequence (e1 , . . . , er ). We try to extend the concatenation of the prefix path of p and er to a simple s-t path. Let er = (v, w). Any valid extension has to avoid the nodes of pref(p) after v, and we are only interested in shortest extensions. Therefore, we compute a new SP tree T and distances d, but in G − pref(p) instead of G to make sure that nodes of the prefix path of p are not used again. If w ∈ / T , pref(p) cannot be extended to a simple s-t path, and we simply discard p. Otherwise, we push the sequence (e1 , . . . , er ) to Q again. In this new sequence, however, we associate T with er instead of Ter from the old sequence. The sequence represents a path p0 obtained by concatenating the simple prefix path 5

of p, the edge er , and the w-t path in T that, by construction, avoids all nodes of pref(p). The suffix itself is simple, too, because it is a shortest path in a subgraph of G. Hence, p0 is simple. The length of this path is c(pref(p)) + c(er ) + d(w). Consider again the example in Figure 1. The sidetrack sequence (a, c) with Ta = Tc = T0 represents a non-simple path p that visits the nodes s, v2 , v1 , v3 , v2 , v1 , v3 , t in this order. The deviation node of p is v3 , its deviation edge c, and its prefix path is (a, (v2 , v1 ), (v1 , v3 )). We compute a new SP tree T in G − pref(p), which only consists of the edge d. Therefore, T does not contain a v2 -t path, and p is discarded. In contrast, assume the sequence (c) with Tc = T0 was just extracted from Q. It represents almost the same path as the sequence above, but it skips the first visit of v2 . Again, v3 is the deviation node and c the deviation edge. The prefix path comprises the nodes s, v1 and v3 . After removing them temporarily, a new shortest path tree T1 is computed, consisting only of the edges b and d. The sequence (c) with Tc = T1 is pushed to Q. This new sequence represents the simple path ((s, v1 ), (v1 , v3 ), c, b, d), i.e., the concatenation of the prefix path of p, the last sidetrack c in the extracted sequence, and the unique v2 -t path in T1 . Finally, when (c) with Tc = T1 is extracted, the represented path is output. The sidetracks emanating from its prefix are (v2 , v1 ) and (v4 , v3 ). Since v1 , v3 ∈ / T1 , these sidetracks are ignored and no new path is pushed to Q. Lemma 1. The above algorithm computes the k shortest simple s-t paths of a weighted, directed graph G = (V, E). Proof. The algorithm uses the same idea of shortest deviations as existing kSSP algorithms or Eppstein’s kSP algorithm. We only have to show that a non-simple path p is processed before its simple enhancement p0 , resulting from the suffix repair in the non-simple case, is actually needed. The set of nodes that are forbidden when the SP tree for p is computed is a proper subset of the node set that the SP tree for p0 may not use. The suffix of p is therefore not longer than that of p0 , and p is extracted from Q (and subsequently, p0 is pushed) before we need to extract p0 . t u In terms of running time, the above algorithm requires too many computations of SP trees. Lemma 2. The running time of the above algorithm is O(km(m + n log n)). Proof. While processing a non-simple path, at most one new path is pushed to Q, which is always simple. Thus, the parent of a non-simple path is always simple. We have to process at most k simple paths, each of which requires O(m + n) running time. Every simple path may have O(m) sidetracks extensions. In the worst case, all of them represent non-simple paths, yielding O(km) SP tree computations with a total running time of O(km(m + n log n)) if Dijkstra’s algorithm with Fibonacci heaps is used. The running time for the non-simple cases clearly dominates. For every subset of E, there is at most one permutation of this subset that represents a simple s-t path. The maximum number of paths enumerated by the 6

algorithm is therefore k 0 := min{k, 2m }. We can limit the size of Q efficiently to k 0 using a double-ended priority queue [12]. We push O(k 0 m) paths to Q and extract O(k 0 m) paths from it; both operations require O(log k 0 ) time on interval heaps. The total time spent on processing Q is then O(k 0 m log k 0 ) ⊂ O(km2 ). The pivot step requires O(n) running time for each of the O(k 0 m) extracted paths. t u

4

Improvements

We show how the number of SP tree computations can be reduced to O(kn) in the worst case, and seemingly even further in practice. So far, we were only able to bound the number of SP tree computations by O(m) for each extracted simple path. This stems from the fact that there may be O(m) sidetracks from such a path, each of them requiring a subsequent SP tree computation in the worst case. Consider two sidetrack sequences (e1 , . . . , er , f1 = (u, v)), (e1 , . . . , er , f2 = (u, w)) that were added when a path p represented by (e1 , . . . , er ) was processed. Let p1 , p2 be the paths represented by these sequences, respectively. Assume that both sequences represent non-simple paths, and therefore both require a new SP tree. We assume w.l.o.g. that p1 is extracted from Q before p2 . When p1 is extracted from Q, we discover that it contains a cycle. We then have to compute an SP tree T for the graph G − p0 , where p0 is the s-u subpath of p. We push (e1 , . . . , er , f1 ) back to Q, updating Tf1 = T . When p2 is extracted, the basic algorithm computes an SP tree for the exact same graph. This computation may be skipped. We check if an SP tree for this graph has already been computed, and reuse it if it exists. In our case, we simply push (e1 , . . . , er , f2 ) with Tf2 = T to Q. We obtain the following result. Lemma 3. Excluding the time spent on Q, the algorithm proposed in Section 3 in conjunction with SP tree reuse requires O(kn(m + n log n)) time to process non-simple paths. Proof. There are still O(km) many sequences in Q that represent non-simple paths, but only O(kn) of them trigger an SP tree computation. Let p be a nonsimple path extracted from Q. The initial pivot step requires time O(n). If each path in Q manages a pointer to its parent path as well as a pointer to the SP tree for G − p0 for every prefix path p0 , already computed SP trees can be accessed in constant time. t u The total running time of O(km2 ) spent on Q is now no longer dominated. Instead of using a priority queue for the candidate paths, we organize all computed paths in a min-heap in the following way. The shortest path is the root of the min-heap. Whenever a path p0 is computed while a path p is processed, we insert p0 into the min-heap as a child of p. Figure 1(b) shows an example of such a min-heap. 7

We want to extract the km smallest elements from this heap using Fredericksons heap selection algorithm [6]. The heap described above has maximum degree m, again yielding a running time of O(km2 ). Let Pp be the set of paths found during the processing of p. Instead of inserting every p0 ∈ Pp as a heap child of p, we heapify Pp to obtain the heap Hp , using the lengths of the paths for keys again. The root of Hp is then inserted into the global min-heap as a child of p. Note that the parent path of every path in Hp is not its heap parent in Hp , but still p itself. Every simple path p in the min-heap now has at most two heap successors with the same parent path as p, and at most one heap successor whose parent is p itself. Every non-simple path has at most one simple path as heap processor. The maximum degree of the global min-heap is therefore bounded by three and Frederickson’s heap selection can be done in time O(km). Corollary 1. The algorithm proposed in Section 3 in conjunction with SP tree reuse and Frederickson’s heap selection algorithm computes the k shortest simple s-t paths of a weighted, directed graph G = (V, E), s, t ∈ V , in O(kn(m+n log n)) time. We propose two more modifications that do not change the asymptotic worstcase running time. We provide evidence that these changes make the algorithm faster in practice. Consider one of the k simple s-t paths p represented by sidetracks (e1 , . . . , er ) with ei = (vi−1 , vi ), s = v0 and t = vr . When p is processed, we push the set Pp of paths to Q, with |Pp | ∈ O(m). The basic algorithm tests for each p0 ∈ Pp if p0 is simple in time O(n), leading to a total time of O(kmn) for these tests. Let T = Ter . By removing all ei from T , the SP tree decomposes into a set of trees Ti such that Ti is rooted in vi . The block i is the node set of Ti . Observe that the path p0 represented by a sequence (e1 , . . . , er , e), e = (vi , w), with vi , w in block i, j, respectively, is simple iff i < j. If i ≥ j, we follow p until we reach vi , traverse e and follow T to reach vi again. Otherwise, the first node on p we hit after deviating from it via e is vj . Since i < j, the vj -t subpath of p does not contain vi , and therefore, p0 is simple. The partition of V into blocks can be computed in time O(n). We can then collect all sidetracks deviating from p and check for each of them if their heads belong to a smaller block than their tails in O(m) total time. We store this information along with the corresponding sidetrack sequences in Q. The pivot turn is replaced by a constant time lookup. All tests for simplicity then require time O(k(m + n)) instead of O(kmn). Finally, we want to reduce the number of SP tree computations in practice. Let p be a non-simple path represented by the sequence (e1 , . . . , er ) with ei = (vi−1 , vi ). After we discover that p is not simple, the basic algorithm computes an SP tree in G − pref(p). Only then does the algorithm check if vr ∈ T . Obviously, there is a shortest w-t path in G − pref(p), i.e., w ∈ T , iff there is some directed path from w to t in G − pref(p). Latter can be checked by a much simpler reachability check in time O(m + n). A naive approach checks reachability for every combination of some node v ∈ V and one of the O(kn) 8

n NC SB-r SB-o NC 4,000 SB-r SB-o NC 6,000 SB-r SB-o NC 8,000 SB-r SB-o NC 10,000 SB-r SB-o 2,000

m = 2n Med Q.9 0.91 2.24 0.18 0.23 0.09 0.17 0.90 2.63 0.35 0.40 0.16 0.21 2.99 5.88 0.54 0.61 0.23 0.30 1.62 7.16 0.69 0.86 0.28 0.46 1.70 10.86 0.87 0.92 0.35 0.40

m = 4n Med Q.9 0.41 1.06 0.27 0.29 0.07 0.09 0.76 2.39 0.53 0.58 0.12 0.17 0.47 1.53 0.81 0.83 0.21 0.26 0.62 2.52 1.08 1.11 0.23 0.26 1.09 5.00 1.37 1.48 0.30 0.42

m = 10n Med Q.9 0.35 1.14 0.72 0.75 0.09 0.12 0.75 1.79 1.39 1.60 0.13 0.20 0.65 3.06 2.11 2.18 0.20 0.27 1.79 4.27 2.81 2.89 0.32 0.41 2.46 8.77 3.63 3.81 0.33 0.38

m = 30n Med Q.9 0.41 1.39 3.70 3.86 0.11 0.16 1.24 4.61 7.26 7.31 0.17 0.22 1.93 7.14 10.95 11.08 0.27 0.37 3.29 9.23 14.66 15.35 0.34 0.54 5.83 12.19 18.42 18.64 0.37 0.44

m = 50n Med Q.9 1.95 3.41 8.94 9.08 0.16 0.29 1.92 4.92 17.46 17.64 0.27 0.38 2.07 8.37 26.53 26.70 0.31 0.50 2.45 8.94 34.84 35.50 0.39 0.47 8.96 23.97 43.59 43.92 0.45 0.61

Table 1: Median and 90% quantile Q.9 of running times in seconds on random graphs with k = 2000.

nodes on the output simple paths, yielding time O(kn2 (m + n)). Of course, for a fixed prefix pi of some simple path p, we can also check reachability in G − pi in O(m + n) time for every node in G simultaneously, and obtain O(kn(m + n)) total time. Let l be the number of edges in p, and consider sidetracks (vl−1 , w). To determine whether the SP path computation for this sidetrack is necessary, we have to check whether there is a path in Gl−1 , where l is the number of edges of p, and Gi := G − pi . We determine reachability in Gl−1 by starting a reverse depth-first search from t, ignoring every node that lies on p. After this search, w-t reachability can be evaluated in constant time per sidetrack (vl−1 , w). If w turns out to be separated from t, the path represented by (e1 , . . . , er , (vl−1 , w)) is non-simple and cannot be repaired to a simple path. In this case, we discard the sidetrack, which in turn cannot trigger an SP tree computation as it is never extracted from Q. After collecting sidetracks emanating from vi+1 on p, we continue with sidetracks emanating from the predecessor vi . We conduct a depth-first search again, this time starting in vi and reusing the reachability information computed before. This way, we only process nodes that were unreachable before. In other words, we solve an incremental series of reachability instances. This procedure terminates when vi is the deviation node of p, and takes total time O(k(m + n)).

5

Experiments

To demonstrate the effectiveness of our algorithm, we conducted a series of experiments. Feng [5] showed recently that their algorithm is the most efficient one in practice. We therefore only compare our algorithm to Feng’s node classification algorithm (NC). Our implementation of NC does not use express edges and achieves better running times than the implementation of Feng, who appears to have used the same processor as we did. We implemented two variants of 9

k

n

NC SB-r SB-o NC 4000 SB-r SB-o NC 2000 6000 SB-r SB-o NC 8000 SB-r SB-o NC 10 000 SB-r SB-o 2000

m = 2n Dijkstra Polls 22 981 4 402 186 80 158 740 82 162 774 25 130 3 284 822 44 177 236 46 185 205 26 990 12 142 440 43 257 184 44 263 158 26 810 5 478 902 25 199 528 28 223 446 26 633 4 806 279 23 229 591 23 229 634

m = 10n Dijkstra Polls 14 532 613 888 65 129 426 65 129 426 14 580 1 465 676 20 77 830 20 77 830 16 652 719 338 18 110 849 18 110 849 17 316 2 418 042 16 131 850 16 131 850 17 826 3 552 732 15 149 870 15 149 870

m = 50n Dijkstra Polls 14 512 2 105 908 44 87 640 44 87 640 15 604 1 404 670 28 113 790 28 113 790 16 444 1 151 469 21 125 826 21 125 826 17 034 1 109 572 17 135 853 17 135 853 18 186 6 130 731 14 134 892 14 134 892

Table 2: Median number of Dijkstra calls and polls for random graphs.

our algorithm, both of which determine simplicity of all sidetracks of a simple path at once by partitioning the nodes into blocks as discussed above. The first version, SB-r, tries to reduce the number of SP tree computations by solving some reachability problems; the second version, SB-o, spares this measure and is thus more optimistic. None of the implementations uses Frederickson’s heap selection algorithm, resulting in an additional running time O(km log k) for SB-r and SB-o, but not for NC. For space restrictions, we contented ourselves with two graph classes that Feng used in their experiments, including road graphs that are especially relevant in practice. We implemented all algorithms in C++, using forward and reverse star representation for directed graphs. Shortest paths (NC) and SP trees (SB-r, SB-o) are computed using a common implementation of Dijkstra’s algorithm; tentative labels are managed by a pairing heap. Our implementation of Dijkstra’s algorithm stops as soon as the label of t is made permanent if only a single pair shortest path is needed, which is essential for NC. The queue of candidate paths Q is implemented as an interval heap, a form of double-headed priority queues, which allows us to limit its size efficiently to the number of simple paths that have yet to be output. Special care has to be taken here for SB-r and SB-o since Q also has to manage non-simple paths. The experiments ran on an Intel Core i7-3770 @ 3.40GHz with 16GB of RAM on a GNU/Gentoo Linux with kernel version 4.2.5 and TurboBoost turned off. Source code was compiled using the GNU C++ compiler g++-4.9.3 and -O3 optimization. 5.1

Random Graphs

We first considered random graphs generated by the sprand generator provided on the website of the ninth DIMACS implementation challenge [1]. The generator draws at random a fixed amount of edges, possibly resulting in a multigraph. For each combination of graph size n ∈ {2000, 4000, 6000, 8000, 10 000} and linear density m/n ∈ {2, 3, 4, 7, 10, 20, 30, 40, 50}, we generated 20 random graphs, and enumerated k ∈ {200, 500, 1000, 2000} simple paths. 10

Graph n m

NY 264 345 733 846

BAY 321 270 800 172

COL 345 666 1 057 066

FLA 1 070 376 2 712 798

(a) Sizes of four road graphs.

NY

BAY

COL

FLA

NC SB-r SB-o NC SB-r SB-o NC SB-r SB-o NC SB-r SB-o

k = 100 Med Q.9 Polls 2.06 12.14 3 918 630 1.38 3.60 528 614 0.55 2.77 528 614 5.11 17.09 15 215 868 1.76 8.43 963 245 0.79 7.62 963 245 6.53 25.13 16 459 836 2.10 18.06 435 666 0.80 18.43 435 666 30.15 67.43 56 950 588 5.53 9.68 1 070 376 2.54 6.81 1 070 376

Med 3.77 2.64 0.97 9.27 3.15 1.49 11.65 4.01 1.42 58.13 10.60 4.65

k = 200 Q.9 Polls 24.11 6 969 812 7.62 792 628 5.91 792 628 33.81 28 784 851 18.52 1 761 843 17.67 1 761 843 44.08 30 058 602 38.40 435 666 37.98 435 666 126.24 107 818 950 33.13 1 070 376 27.72 1 070 376

Med 5.40 3.88 1.38 13.84 4.77 1.99 15.98 6.00 2.02 83.00 15.73 6.78

k = 300 Q.9 Polls 35.17 9 803 818 11.43 1 055 890 9.10 1 055 890 49.54 38 731 878 28.18 1 922 360 24.97 1 922 360 58.83 42 430 752 62.28 435 666 60.62 435 666 188.03 151 950 959 55.29 1 070 376 47.00 1 070 376

(a) Median and 90% quantile Q.9 of running times in seconds, median number of polls.

Table 5: Sizes and metrics for four large TIGER road graphs.

In Table 1, the median and 90% quantile Q.9 (90% of the running times were at most Q.9 ) of execution times for some densities are summarized. For small densities, we observe that SB-r is faster than NC, but becomes much slower as the density grows. The running time of SB-r increases by a factor of 50 between m = 2n and m = 50n, but only by a factor of 2 for NC. SB-o is about twice as fast as SB-r for very small densities, and is even more robust against density changes than NC. SB-o is thus the fastest of the three algorithms for all graph sizes and densities. Also note the very low dispersion of SB running times. For NC, the 90% quantile of the running time is regularly three times the median running time, and even exceeds a factor of 6 for n = 10 000 and m = 20 000. In contrast, this quotient is always much closer to 1 for SB-o and assumes its maximum of 1.87 for n = 2000, m = 100 000. We can therefore predict the running time of SB-o much more accurate than that of NC. The Q.9 running time of SB-o is still well below the median running time of both SB-r and NC. The dispersion of SB-r is even lower. Table 2 shows the median number of Dijkstra calls. The numbers are relatively stable across the various densities, but the Dijkstra counts for the SB algorithms is orders of magnitudes smaller than the count for the NC algorithm. Note, however, that SB needs to compute the complete SP tree every time. In contrast, NC only solves single pair shortest path problems on rather small subgraphs. We also provide the number of polls, i.e., the total number of nodes that were extracted from Dijkstra’s priority queue, for comparability. NC still requires an order of magnitude more polls compared to SB. Further, there are only small differences in the median number of polls for SB-r and SB-o for m = 2n. Solving incremental reach for each of the 2000 output simple paths only reduced the number of SP tree computations by at most three for n = 8000, m = 2n. The 11

extra effort involved to reduce SP tree computations ranges from 50% of the total running time on very sparse graphs to 98% on very dense graphs, and does not pay off. Therefore, SB-o is clearly the fastest algorithm on random graphs. 5.2

Road Graphs

We considered road graphs of various areas in the USA called TIGER graphs, again provided by the DIMACS website [1]. In particular, we used the road networks of New York (NY), the San Francisco Bay Area (BAY), Colorado (COL), and Florida (FLA). The sizes of these graphs are shown in Table 3(a). We drew 20 s-t pairs at random and enumerated k ∈ {100, 200, 300} paths. The resulting running times are summarized in Table 4(a), along with the median number of polls. The median running time of NC is clearly dominated by both SB variants. SB-o achieves a minimum speedup around 4 on NY across all values of k; on FLA, the speedup is roughly 12. SB-r takes approximately twice the time of SB-o and is still much faster than NC. The ratio of Q.9 and median running time is worse for SB, but the Q.9 time itself of SB-o is still better than that of NC except for 300 paths on COL. On FLA, the largest graph, the 90% quantile of SB-o is much better than the median running time of NC.

References 1. The ninth dimacs implementation challenge: 2005-2006. http://www.dis. uniroma1.it/challenge9/, accessed: 2015-11-12 2. Clarke, S., Krikorian, A., Rausen, J.: Computing the N best loopless paths in a network. J. SIAM 11(4), 1096–1102 (1963) 3. Eppstein, D.: Finding the k shortest paths. SIAM J. Comput. 28(2), 652–673 (1998) 4. Eppstein, D.: K -best enumeration (2014), arXiv: 1412.5075v1 [cs.DS] 5. Feng, G.: Finding k shortest simple paths in directed graphs: A node classification algorithm. Networks 64(1), 6–17 (2014) 6. Frederickson, G.N.: An optimal algorithm for selection in a min-heap. Inf. Comput. 104(2), 197–214 (1993) 7. Gotthilf, Z., Lewenstein, M.: Improved algorithms for the k simple shortest paths and the replacement paths problems. Inf. Process. Lett. 109(7), 352–355 (2009) 8. Hershberger, J., Maxel, M., Suri, S.: Finding the k shortest simple paths: A new algorithm and its implementation. ACM Transactions on Algorithms 3(4) (2007) 9. Katoh, N., Ibaraki, T., Mine, H.: An efficient algorithm for K shortest simple paths. Networks 12(4), 411–427 (1982) 10. Pascoal, M.M.B.: Implementations and empirical comparison of K shortest loopless path algorithms (2006), 9th DIMACS Implementation Challenge Workshop: Shortest Paths 11. Pettie, S.: A new approach to all-pairs shortest paths on real-weighted graphs. Theor. Comput. Sci. 312(1), 47–74 (2004) 12. Sahni, S.: Data Structures, Algorithms, and Applications in C++. McGraw-Hill Pub. Co., 1st edn. (1999) 13. Vassilevska Williams, V., Williams, R.: Subcubic equivalences between path, matrix and triangle problems. In: FOCS 2010. pp. 645–654 (2010) 14. Yen, J.Y.: Finding the k shortest loopless paths in a network. Networks 17(11), 712–716 (1971)

12

13

m = 3n Med Q.9 0.58 1.18 0.22 0.26 0.07 0.11 0.58 2.04 0.44 0.49 0.13 0.19 0.49 3.78 0.66 0.68 0.20 0.22 1.98 4.97 0.88 0.96 0.26 0.34 1.56 7.75 1.09 1.25 0.30 0.47

m = 3n Med Q.9 0.58 1.18 0.04 0.09 0.58 2.04 0.07 0.12 0.49 3.78 0.10 0.16 1.98 4.97 0.12 0.21 1.56 7.75 0.10 0.26

m = 2n Med Q.9 0.91 2.24 0.18 0.23 0.09 0.17 0.90 2.63 0.35 0.40 0.16 0.21 2.99 5.88 0.54 0.61 0.23 0.30 1.62 7.16 0.69 0.86 0.28 0.46 1.70 10.86 0.87 0.92 0.35 0.40

m = 2n Med Q.9 0.91 2.24 0.07 0.15 0.90 2.63 0.09 0.13 2.99 5.88 0.15 0.22 1.62 7.16 0.11 0.33 1.70 10.86 0.12 0.20

m = 4n Med Q.9 0.41 1.06 0.05 0.07 0.76 2.39 0.06 0.10 0.47 1.53 0.10 0.14 0.62 2.52 0.10 0.14 1.09 5.00 0.10 0.23

m = 4n Med Q.9 0.41 1.06 0.27 0.29 0.07 0.09 0.76 2.39 0.53 0.58 0.12 0.17 0.47 1.53 0.81 0.83 0.21 0.26 0.62 2.52 1.08 1.11 0.23 0.26 1.09 5.00 1.37 1.48 0.30 0.42

m = 4n Med Q.9 0.10 0.27 0.07 0.07 0.02 0.02 0.19 0.60 0.13 0.15 0.03 0.05 0.12 0.39 0.20 0.21 0.05 0.06 0.15 0.62 0.27 0.28 0.06 0.07 0.28 1.24 0.34 0.38 0.07 0.11

m = 7n Med Q.9 0.30 0.81 0.05 0.09 0.72 1.78 0.10 0.18 1.20 2.93 0.13 0.17 1.05 6.24 0.12 0.15 1.56 9.06 0.12 0.17

m = 7n Med Q.9 0.30 0.81 0.46 0.50 0.06 0.11 0.72 1.78 0.95 1.02 0.16 0.23 1.20 2.93 1.36 1.45 0.22 0.27 1.05 6.24 1.85 1.87 0.24 0.27 1.56 9.06 2.33 2.39 0.30 0.37

m = 7n Med Q.9 0.06 0.20 0.11 0.13 0.02 0.03 0.20 0.55 0.24 0.26 0.04 0.06 0.30 0.75 0.34 0.36 0.05 0.06 0.27 1.53 0.46 0.47 0.06 0.07 0.38 2.25 0.58 0.60 0.08 0.09

m = 20n Med Q.9 0.11 0.39 0.49 0.51 0.03 0.05 0.71 1.21 0.96 0.99 0.04 0.06 0.28 1.87 1.44 1.56 0.06 0.11 0.36 2.06 1.94 2.18 0.07 0.09 0.80 2.35 2.44 2.67 0.10 0.13

m = 20n Med Q.9 0.46 1.54 1.97 2.04 0.10 0.17 2.88 4.78 3.85 3.95 0.17 0.26 1.18 7.44 5.79 6.46 0.22 0.35 1.40 8.22 7.77 8.57 0.30 0.34 3.18 9.55 9.83 10.28 0.37 0.47

m = 20n Med Q.9 0.46 1.54 0.10 0.16 2.88 4.78 0.14 0.21 1.18 7.44 0.16 0.28 1.40 8.22 0.17 0.23 3.18 9.55 0.19 0.29

k = 500 m = 10n Med Q.9 0.08 0.27 0.18 0.19 0.02 0.03 0.18 0.43 0.35 0.38 0.03 0.05 0.17 0.78 0.52 0.54 0.05 0.07 0.44 1.06 0.71 0.75 0.07 0.09 0.61 2.09 0.91 0.99 0.08 0.10 k = 2000 m = 10n Med Q.9 0.35 1.14 0.72 0.75 0.09 0.12 0.75 1.79 1.39 1.60 0.13 0.20 0.65 3.06 2.11 2.18 0.20 0.27 1.79 4.27 2.81 2.89 0.32 0.41 2.46 8.77 3.63 3.81 0.33 0.38 k = 2000 m = 10n Med Q.9 0.35 1.14 0.08 0.12 0.75 1.79 0.08 0.14 0.65 3.06 0.12 0.19 1.79 4.27 0.13 0.21 2.46 8.77 0.14 0.19

m = 30n Med Q.9 0.41 1.39 0.11 0.16 1.24 4.61 0.15 0.18 1.93 7.14 0.21 0.29 3.29 9.23 0.20 0.33 5.83 12.19 0.21 0.29

m = 30n Med Q.9 0.41 1.39 3.70 3.86 0.11 0.16 1.24 4.61 7.26 7.31 0.17 0.22 1.93 7.14 10.95 11.08 0.27 0.37 3.29 9.23 14.66 15.35 0.34 0.54 5.83 12.19 18.42 18.64 0.37 0.44

m = 30n Med Q.9 0.10 0.35 0.92 0.95 0.03 0.04 0.31 1.14 1.81 1.83 0.04 0.06 0.48 1.78 2.74 2.76 0.07 0.09 0.82 2.20 3.64 3.72 0.08 0.12 1.45 3.01 4.61 4.62 0.09 0.12

m = 40n Med Q.9 0.99 2.22 0.15 0.20 1.19 5.00 0.19 0.28 1.58 7.28 0.21 0.26 3.04 10.36 0.26 0.36 5.92 15.18 0.26 0.33

m = 40n Med Q.9 0.99 2.22 6.10 6.32 0.14 0.20 1.19 5.00 11.87 12.56 0.21 0.34 1.58 7.28 17.62 17.84 0.26 0.32 3.04 10.36 23.78 24.95 0.38 0.47 5.92 15.18 29.76 29.86 0.41 0.49

m = 40n Med Q.9 0.23 0.54 1.50 1.53 0.03 0.04 0.29 1.24 2.96 3.60 0.05 0.08 0.39 1.86 4.40 4.44 0.07 0.09 0.76 2.56 5.93 6.21 0.09 0.13 1.48 3.81 7.42 7.59 0.10 0.13

m = 50n Med Q.9 1.95 3.41 0.17 0.29 1.92 4.92 0.22 0.32 2.07 8.37 0.25 0.41 2.45 8.94 0.28 0.35 8.96 23.97 0.33 0.44

m = 50n Med Q.9 1.95 3.41 8.94 9.08 0.16 0.29 1.92 4.92 17.46 17.64 0.27 0.38 2.07 8.37 26.53 26.70 0.31 0.50 2.45 8.94 34.84 35.50 0.39 0.47 8.96 23.97 43.59 43.92 0.45 0.61

m = 50n Med Q.9 0.47 0.84 2.22 2.26 0.04 0.07 0.48 1.27 4.34 4.74 0.05 0.08 0.51 2.09 6.56 6.59 0.07 0.11 0.62 2.46 8.72 9.69 0.09 0.11 2.25 5.89 10.88 10.93 0.11 0.15

Table 6: Median and 90% quantile Q.9 of running times in seconds of NC, SB-r and SB-o on random graphs of all tested sizes for k ∈ {200, 500, 1000}.

2000

NC SB-o NC 4000 SB-o NC 6000 SB-o NC 8000 SB-o NC 10 000 SB-o

n

2000

NC SB-r SB-o NC 4000 SB-r SB-o NC 6000 SB-r SB-o NC 8000 SB-r SB-o NC 10 000 SB-r SB-o

n

2000

NC SB-r SB-o NC 4000 SB-r SB-o NC 6000 SB-r SB-o NC 8000 SB-r SB-o NC 10 000 SB-r SB-o

n

m = 3n Med Q.9 0.14 0.29 0.05 0.07 0.02 0.03 0.15 0.51 0.11 0.13 0.03 0.05 0.13 0.95 0.16 0.17 0.05 0.06 0.51 1.21 0.22 0.24 0.06 0.09 0.36 1.97 0.27 0.31 0.08 0.11

m = 2n Med Q.9 0.22 0.56 0.04 0.06 0.02 0.04 0.22 0.65 0.09 0.10 0.04 0.04 0.74 1.47 0.14 0.15 0.06 0.07 0.40 1.78 0.17 0.20 0.07 0.10 0.42 2.73 0.21 0.23 0.08 0.10

45 40 35 30

NC SB-r SB-o

4

7

10

20

linear density m/n

30

40

14

25 20 15 10 5 0

2

Fig. 2: Boxplots of running times in seconds for random graphs with 10 000 nodes and k = 2000. Plus signs represent outliers. A red square marks the mean, but was omitted for the SB algorithms as they would completely cover their corresponding boxes.

time[s]

50

15

100

NC SB-r SB-o

100

NC SB-r SB-o

200 enumerated paths

COL

200 enumerated paths

NY

300

300

0

50

100

150

200

250

0

50

100

150

200

250

100

NC SB-r SB-o

100

NC SB-r SB-o

200 enumerated paths

FLA

200 enumerated paths

BAY

Fig. 3: Boxplots of running times in seconds for TIGER graphs. Plus signs represent outliers. A red square marks the mean. The interquartile range of SB on FLA for k ∈ {100, 200} is so small that the boxes appear as lines.

160 140 120 100 80 60 40 20 0

time [s]

time [s]

40 35 30 25 20 15 10 5 0 time [s] time [s]

300

300