Algorithm Engineering of SSSP with Negative Edge Weights

We initiate the transfer of the theoretical break-through on the problem of shortest paths with negative edge weights [5, 6] to the application domain and provide the first implementation inspired by the work of [6], to the best of our knowledge. In particular, our work primarily aims to demonstrate the potential of this novel approach, suggesting avenues for future research to establish its competitive advantage over classical Bellman-Ford based algorithms. We first present an overview of the algorithm of [6] to make it easily accessible to a broader audience and especially practitioners.⁴⁴4Naturally, explaining the running time analysis of the algorithm of [6] is out of the scope of this paper and we merely focus on the algorithm. We discuss different adaptions that are necessary to obtain an implementation that is competitive with existing algorithms. We also show an improved instance-based upper bound on a crucial parameter that controls the recursion depth of the algorithm. We create a set of benchmark instances that is inspired by the work of [10], adapting their benchmark instances to make them more difficult to solve for our implementation in order to allow for a fair comparison. As the algorithm of [6] contains multiple constants that can be chosen in different ways and our modifications also allow for different parameter settings, we conduct experiments to empirically determine a good parameter choice. We then conduct extensive experiments to understand the running time and especially also the scaling behavior of our implementation compared to one of the state-of-the-art algorithms, namely GOR [21]. On the hardest instance type, we are faster by up to almost two orders of magnitude, suggesting that our implementation particularly thrives on hard instances. Indeed, the experiments also suggest that the running time of our implementation scales near-linearly.

We note that Li and Mowry [25] very recently and independent of our work further simplify the algorithm of [6]. However, they do not provide an implementation and conduct no experimental evaluation.

In this work, we are always given a directed, weighted graph $G=(V,E,w)$ equipped with a weight function $w:E\mapsto \mathbb{Z}$ that assigns (potentially negative) integer weights to the edges. With slight abuse of notation, we use $|G|$ to denote the number of vertices of a graph $G$. Given a subset of vertices $C\subset V$, we denote with $G[C]$ the induced subgraph of $C$ in $G$. We sometimes consider the graph $G_{\ge 0}$, which we define as the graph $G$ with the modified cost function $w_{\ge 0}(e)=\max \{w(e),0\}$ where all negative edge weights are set to zero. We also need access to the reverse graph. We hence denote with $G^{\mathrm{out}}$ the original graph $G$, and with $G^{\mathrm{in}}$ the graph $G$ but where each edge is flipped in its orientation. We use the notation $B_{\ge 0}^{\mathrm{dir}}(v,r)$ to denote the set of vertices that have distance at most $r$ from $v$ in $G_{\ge 0}^{\mathrm{dir}}$, for $\mathrm{dir}\in \{\mathrm{in},\mathrm{out}\}$. We say that $B_{\ge 0}^{\mathrm{dir}}(v,r)$ is an $r$-ball of $v$. Given such a ball $B≔B_{\ge 0}^{\mathrm{dir}}(v,r)$, we denote its boundary by $\partial B$, which comprises of all the edges that lead from nodes inside the ball to nodes outside the ball. Finally, given a direction $\mathrm{dir}\in \{\mathrm{in},\mathrm{out}\}$, we denote by $\overline{\mathrm{dir}}$ the opposite direction.

The goal of this work is to solve the following problem efficiently in practice:

SSSP with Negative Weights reduces to finding a potential function $\varphi :V\mapsto \mathbb{Z}$ such that the modified weights $w(u,v)+\varphi (u)-\varphi (v)$ of all edges $(u,v)\in E$ are nonnegative; we call such a potential function valid. This is also called “Johnson’s Trick” [22]. There are three important facts about such potential functions:

• $\mathrm{■}$

  a valid potential function $\varphi$ always exists iff $G$ does not contain a negative cycle, and

• $\mathrm{■}$

  $\varphi$ preserves the structure of the shortest paths in $G$, and

• $\mathrm{■}$

  given shortest path costs w.r.t. $\varphi$, we can easily reconstruct the original costs.

Especially, if we found a valid potential function $\varphi$, we can just apply Dijkstra’s algorithm [14] on $G$ with the modified weights and hence obtain our shortest path distances in the graph with the original weights.

Consider any cycle $C=e_1,\mathrm{\dots },e_{\mathrm{\ell }}$ in $G$, its mean is defined as ${\sum }_{i\in [\mathrm{\ell }]}w(e_i)/\mathrm{\ell }$. The minimum cycle mean is defined as the minimum mean any cycle in $G$ has. In this work we focus on restricted graphs.

Note that this definition also implies that $G$ has no negative cycle. Furthermore, the original definition also includes the existence of a super source, due to technical reasons. We drop this requirement.

We mostly restrict to these types of graphs in this work (unless mentioned otherwise), due to the following reasons:

• $\mathrm{■}$

  The main contribution of [6] is an algorithm that runs on restricted graphs. How to generalize such an algorithm to the general case is technical, but comparably simple, and was already previously known.⁵⁵5Our implementation can easily be adapted to general negative-cycle-free graphs by the scaling technique of [6]. On the other hand, the negative cycle detection of [6] is impractical as it works via budgeting – a technique that, very roughly speaking, given the theoretical worst-case running time of the algorithm on well-formed inputs, will stop if it is exceeded. We discuss negative cycles more in Section 3. Hence, restricting to these types of graphs enables for a cleaner evaluation of how their techniques perform in practice.

• $\mathrm{■}$

  Most of our benchmark graphs are derived from the hard instances described in [10]. These graphs are restricted or almost restricted and thus, even if we adapt our implementation for the general case, we expected the performance on these graphs to not significantly change.

• $\mathrm{■}$

  The restriction does not have any effect on the correctness of our implementation: even if our implementation is given a non-restricted graph, it still outputs the correct result.⁶⁶6This is the case as Line 9 in Algorithm 4 computes the correct result also on non-restricted graphs even with negative cycles. Our implementation is merely not optimized for these cases and might have higher running times. In particular, if the graph contains a negative cycle, it still will eventually be detected by our implementation.

In this subsection we give a description of [6, Algorithm 2]. For the readers’ convenience, we refer to the statements and algorithms in the freely available full version [7] of [6]. Note that in the algorithm of [6] a designated source $s\in V$ is added and connected to all vertices via zero-weight edges. In our description and our implementation, we replace the explicit source node by an implicit source. Computing a shortest path tree from this implicit source then gives us a valid potential function, and we can subsequently run Dijkstra’s algorithm from any source with the modified weights, see Section 2.1.

As the goal of this work is to understand the practicability of the algorithmic approach of [6] and [5], we use the structure described in Section 2.3, and in particular Algorithm 4, as basis and modify it to obtain a practically faster algorithm.

As a first preprocessing step, we decompose the graph into strongly connected components. This allows us to run our algorithm separately on each component and later obtain a valid potential function that also makes the edges in-between the components positive using FixDAGEdges. To this end, we also have to choose one of the many algorithms to compute strongly connected components. We use Kosaraju’s and Sharir’s algorithm [30], which also gives us a topological sorting of the strongly connected components as a by-product.

A subroutine that is pervasive in [6] is Dijkstra’s algorithm. In [6] the authors use Dijkstra’s algorithm with a priority queue by Thorup [31] with the goal to shave an additional logarithmic factor. We instead use a classical priority queue, more precisely, we use a 4-heap that has proven to be fast in practice in shortest path problems with non-negative edge weights, see e.g. [20].

A critical parameter that controls the size of the components that we recurse on and thereby also the recursion depth is $\kappa$ – hence choosing $\kappa$ well can have significant influence on the practical running time. In Section 4, we show that using $\mathrm{diam}(G_{\ge 0})$ is an upper bound to $\kappa$ on any restricted graph $G$. As $\mathrm{diam}(G_{\ge 0})$ could be greater than $n$, we set $\kappa =\min \{n,\mathrm{diam}(G_{\ge 0})\}$ in our implementation. It is also possible that the graph contains only very few negative edges in total, so we could additionally have the total number of negative edges in $G$ as third argument in the $\min$ above. However, as all the instances that we run our experiments on contain a significant amount of negative edges, we do not use this optimization in the experiments.

The computation of the estimated set of light vertices $L$ in Decompose is an important factor in the algorithm. To reduce the running time, we can reduce the number of randomly sampled vertices that we use to estimate the lightness. While the choice in the work of [6] is aimed at ensuring a high probability of success, it is conceivable that on practical instances we need less vertices for a good estimation. To this end, we introduce a parameter $K\in \mathbb{N}$ to our algorithm such that we choose $1/K$-times as many vertices to estimate the lightness in the Decompose routine.

We also have to decide at which point we want to invoke the base case and call LazyDijkstra. In [6] this is done when $\kappa \le 2$, see Algorithm 4. While this is theoretically guaranteed to occur within a logarithmic number of recursive calls with sufficiently high probability, in practice this can slow down the implementation significantly. Hence, we settled on calling the base case if $n+\kappa \le 300$. We use the sum as we want to only call the base case if both $\kappa$ and $n$ are small. Note that due to how we set $\kappa$, we also always have $\kappa \le n$.

Finally, our implementation also performs negative cycle detection, however, this is not its strong suit. In [6] negative cycle detection is done via budgeting, i.e., when the algorithm exceeds a specific running time, then the existence of a negative cycle is reported; this technique is impractical and we do not use it. For an alternative approach, note that the only cases where negative cycles can occur in Algorithm 4 are the LazyDijkstra calls. Hence, we add negative cycle detection to LazyDijkstra by merely checking the number of distance improvements that any vertex sees. If this number exceeds the number of vertices, then we detected the existence of a negative cycle and we report this. How to adapt our implementation to also perform efficient negative cycle detection is an intriguing task for future work.⁷⁷7We note that in a very recent parallel and independent work that appeared on arXiv, an improved cycle detection is presented, but this work does not present any implementation or experiments [25].

Let us consider what asymptotic running time we can hope for our implementation. Note that we use the wording “hope for” on purpose, as the above modifications can potentially break the theoretical guarantees and our statements should be considered an intuitive sketch. The theoretically expected running time for the version of Algorithm 4 in [6] is $O((m+n\log \log n){\log }^{2}n)$, where $n$ and $m$ are the number of vertices and edges in the input graph, respectively. We perform one crucial modification, that is, we use a $4$-heap instead of the priority queue by [31]. This replaces the $\log \log n$ factor in the above running time by a $\log n$ factor. Hence, we can hope for an $O((m+n\log n){\log }^{2}n)$ running time.

In our implementation, following [6], we use LazyDijkstra as algorithm to compute a valid potential function on graphs where each shortest path should only contain a small number of negative edges. This happens in the base case of Algorithm 4 as well as in the case of obtaining a valid potential function for the cut edges in-between components. Instead of LazyDijkstra we could actually employ any other shortest path algorithm for graphs with negative edge weights that is fast in practice. Especially, we can also use the baseline algorithm that we later compare to in Section 5 as subroutine.

We wrote the entire experimental framework and all the algorithms in C++20. We compiled our code using GCC on Debian. We ran the experiments on a server with 48 Intel Xeon E5-2680 v3@2.50GHz CPUs with 256 GB of RAM. We note that we did not parallelize our implementation as our comparison is with a single-threaded implementation, but parallelization of our implementation is straight-forward (the recursive calls are independent).

We compare our implementation against a state-of-the-art algorithm by Goldberg and Radzik [21] called GOR. More specifically, we use the implementation in SPLIB⁸⁸8See https://www3.cs.stonybrook.edu/~algorith/implement/goldberg/distrib/splib.tar. We modified it such that we can use it with our graph data structure and compile it within our implementation. We ensured that this did not incur a significant slowdown., which was developed by Cherkassky, Goldberg, and Radzik. We also implemented the algorithm called BFCT from [10] (a previous version was described in [11]). However, as our implementation was consistently outperformed by GOR, we only compare to GOR here. Our reason for choosing GOR is twofold:

• $\mathrm{■}$

  In [10], the authors show that GOR has overall the fastest average running times for scans of a single vertex and a good performance on all instances regarding the number of scans required to compute a shortest path tree, without large deviations.

• $\mathrm{■}$

  A highly-tuned implementation of GOR to compare to is publicly available.

In the following, we present a brief overview of all the different types of instances we perform experiments on.

As discussed in Section 3, there are several choices to be made and parameters to be set. We list the interesting and non-obvious choices and parameters here and empirically evaluate them to understand what good choices are. Note that we set and evaluated the parameters in the order described here, i.e., the best choice of the first mentioned parameter is already used in the subsequent experiments. Hence, we chose the order below carefully and on purpose.

As described in Section 3.1, we can change the number of random vertices that we use to estimate the lightness in Algorithm 3. To that end, our implementation has a parameter $K\in [1,\mathrm{\infty }]$ such that we sample $\max \{1,k/K\}$ random vertices, where $k$ is the theoretical number of vertices as used in Algorithm 3. We tested OUR with $K\in \{1,5,10,20,40,\mathrm{\infty }\}$ on several large instances, see Table 1. The speed-up induced by this parameter is significant; it is up to almost a factor of 4. We notice that the worst choice is $K=1$, i.e., the original value that is used in Algorithm 4. There is no clear best choice between $K$ being $20$, $40$, or $\mathrm{\infty }$. To strike a balance, we use $K=40$ in our subsequent experiments.

Next we consider whether we should use the improved upper bound on $\kappa$ in our implementation, see Section 4 and Section 3.1. We call this method the diameter upper bound here. In almost all of our experiments we noticed a similar scaling behavior with the diameter upper bound and without it. On one hand, on small instances computing the diameter sometimes poses a significant running time overhead and makes it a bit slower in total. On the other hand, on the RANDOM RESTRICTED instances, using the diameter upper bound leads to a better scaling behavior and a speed-up of almost an order of magnitude, see Figure 2. Hence, despite small disadvantages on some instances, we decide to use it as it can massively speed-up some cases.

Finally, as discussed in Section 3.1, we can also replace LazyDijkstra by any other algorithm that performs shortest path computations on graphs with negative edge weights. The obvious replacement to evaluate is GOR. While the speed-up is minor when using GOR as replacement, the slowdown is massive on AUG-GOR instances, see Figure 2. Hence, we conclude that using LazyDijkstra as base algorithm is the more robust choice.

In this section we compare the performance of OUR and GOR on several instances. The size of the instances that we use ranges from $5\cdot {10}^5$ edges to $2\cdot {10}^7$ edges. We consider the total running times but also the scaling behavior of both algorithms.

Let us first consider the behavior of OUR and GOR on the RANDOM RESTRICTED instances, see Figure 4. The scaling behavior of OUR is only slightly better than GOR, however OUR is consistently slower by a significant factor. We note, however, that the overall running times for both algorithms on these instances are low compared to our other experiments on graphs with similar sizes. Let us now consider the DIMACS instance. We compare the running time of OUR and GOR on the unmodified instance, which we refer to as USA-0, and on the instance with a potential shift $W$ (see Section 5.3) of $1$, $10$, and $100$ (USA-1, USA-10, USA-100, respectively). On USA-0 OUR is $5.4$ time slower. On USA-1 and USA-10 OUR is $7.1$ times slower, and on USA-100 OUR is $6.1$ time slower.

We now consider experiments on the BAD instances described in Section 5.3. As discussed in Section 5.3, a comparison on these instances is misleading as our implementation explicitly exploits the DAG structure while GOR does not. However, it is still enlightening to see for which instances GOR performs poorly. We show the results of these experiments in Table 2. We can see that OUR has a stable running time on these instances, as expected. Furthermore, we can clearly see the running time overhead of the data structures of OUR over the more lightweight GOR. However, the instances BAD-GOR and BAD-RDB are clearly pathological for GOR and it takes excessive time to solve these two instances, while the running time of OUR is similar to the one for other instances. We believe that an approach similar to ours that reduces the overhead induced by the data structures could become competitive to GOR on the instances where it is faster.

We now perform experiments on the AUG instances, see Figures 3(a) to 3(e). First, note that on all AUG instances the scaling behavior of OUR is better than the scaling behavior of GOR. In particular, the regression results for the running times of OUR on all of the restricted instances (i.e., Figures 3(b) to 3(e)) are compatible with the theoretically expected near-linear running time: As described in Section 3.1, we could expect a running time of $O((m+n\log n){\log }^{2}n)$, which is $O(n{\log }^{3}n)$ for $m\in O(n)$. The same regression that we use for the running times applied to $x{\log }^3(x)$ for sampled $x\in [5\cdot {10}^5,2\cdot {10}^7]$ leads to a polynomial factor of $x^{1.20\pm 0.01}$, which is very close to the experimental results in the restricted cases. On the non-restricted augmented instance AUG-GOR, the running time of OUR is significantly worse than on the restricted instances, even though it is still significantly better than the running time of GOR. When it comes to absolute running times, our experiments show a diverse picture. On the AUG-GOR and AUG-RD instances, OUR is always faster than GOR, see Figures 3(a) and 3(d). However, on the AUG-BFCT and AUG-RDB instances, GOR is always faster than OUR, see Figures 3(b) and 3(e). On AUG-DFS, OUR is slower on the smallest instance and consistently faster on all the other sizes, see Figure 3(c). Summarizing, we want to stress that while GOR indeed is very fast on some instances, it also exhibits fluctuations of multiple orders of magnitudes on graphs of similar size due its adaptiveness; we do not see such fluctuations in our implementation when running it on restricted instances.

We now consider the experiments on the SHIFT-GOR instance, see Figure 3(f). This instance is particularly hard for GOR, as it highlights its quadratic nature. Despite being far from restricted, OUR performs consistently better. OUR outperforms GOR by almost two orders of magnitudes on some of these non-trivial instances. It is to notice, however, that the running time of OUR is strongly unstable, probably because of the strong non-restricted nature of this instance. Each point in the plot in Figure 3(f) is the average of the running time of $5$ different random generations of an instance of that size (and only one run per instance). This only marginally smoothed the unstable trend.