Here we go, my first blog entry, ever. Let’s start with two short quizzes.
1. The well known Dijkstra’s algorithm is:
[a] A greedy algorithm
[b] A dynamic programming algorithm
[c] A primal-dual algorithm
[d] It was discovered by Dantzig
2. Which is the best C++ implementation of Dijkstra’s algorithm among the following?
[a] The Boost Graph Library (BGL)
[b] The COIN-OR Lemon Graph Library
[c] The Google OrTools
[d] Hei dude! We can do better!!!
What is your answer for the first question? … well, the answers are all correct! And for the second question? To know the correct answer, sorry, you have to read this post to the end…
If you are curious to learn more about the classification of the Dijkstra’s algorithm proposed in the first three answers, please consider reading [1] and [2]. Honestly, I did not know that the algorithm was independently discovered by Dantzig [3] as a special case of Linear Programming. However, Dantzig is credited for the first version of the bidirectional Dijkstra’s algorithm (should we called it Dantzig’s algorithm?), which is nowadays the best performing algorithm on general graphs. The bidirectional Dijkstra’s algorithm is used as benchmark to measure the speed-up of modern specialized shortest path algorithms for road networks [4,5], those algorithms that are implemented, for instance, in our GPS navigation systems, in yours smartphones (I don’t have one, argh!), in Google Maps Directions, and Microsoft Bing Maps.
Why a first blog entry on Dijkstra’s algorithm? That’s simple.
- Have you ever implemented an efficient version of this well-known and widely studied algorithm?
- Have you ever used the version that is implemented in well-reputed graph libraries, such as, the Boost Graph Library (BGL), the COIN-OR Lemon, and/or Google OrTools?
I did while programming in C++, and I want to share with you my experience.
The Algorithm
The algorithm is quite simple. First partition the nodes of the input graph G=(N,A) in three sets: the sets of (1) scanned, (2) reachable, and (3) unvisited nodes. Every node has a distance label and a predecessor vertex . Initially, set the label of the source node , while set for all other nodes. Moreover, the node s is placed in the set of reachable nodes, while all the other nodes are unvisited.
The algorithm proceedes as follows: select a reachable node i with minimum distance label, and move it in the set of scanned nodes, it will be never selected again. For each arc (i,j) in the forward star of node i check if node j has distance label ; if it is the case, update the label and the predecessor vertex . In addition, if the node was unvisited, move it in the set of reachable nodes. If the selected node i is the destination node t, stop the algorithm. Otherwise, continue by selecting the next node i with minimum distance label.
The algorithm stops either when it scans the destination node t or the set of reachable nodes is empty. For the nice properties of the algorithm consult any textbook in computer science or operations research.
At this point it should be clear why Dijkstra’s algorithm is greedy: it always select a reachable node with minimum distance label. It is a dynamic programming algorithm because it maintains the recursive relation for all . If you are familiar with Linear Programming, you should recognize that the distance labels play the role of dual variable of a flow based formulation of the shortest path problem, and the Dijkstra’s algorithm costructs a primal solution (i.e. a path) that satisfies the dual constraints .
Graphs and Heaps
The algorithm uses two data structures: the input graph G and the set of reachable nodes Q. The graph G can be stored with an adjacency list, but be sure that the arcs are stored in contiguous memory, in order to reduce the chance of cache misses when scanning the forward stars. In my implementation, I have used a std::vector to store the forward star of each node.
The second data structure, the most important, is the priority queue Q. The queue has to support three operations: push, update, and extract-min. The type of priority queue used determines the worst-case complexity of the Dijkstra’s algorithm. Theoretically, the best strongly polynomial worst-case complexity is achieved via a Fibonacci heap. On road networks, the Multi Bucket heap yields a weakly polynomial worst-case complexity that is more efficient in practice [4,5]. Unfortunately, the Fibonacci Heap is a rather complex data structure, and lazy implementations end up in using a simpler Binomial Heap.
The good news is that the Boost Library from version 1.49 has a Heap library. This library contains several type of heaps that share a common interface: d-ary-heap, binomial-heap, fibonacci-heap, pairing-heap, and skew-heap. The worst-case complexity of the basic operations are summarized in a nice table. Contrary to text-books, these heaps are ordered in non increasing order (they are max-heap instead of min-heap), that means that the top of the heap is always the element with highest priority. For implementing Dijkstra, where all arc lengths are non negative, this is not a problem: we can store the elements with the distance changed in sign (sorry for the rough explanation, but if you are really intrested it is better to read directly the source code).
The big advantage of boost::heap is that it allows to program Dijkstra once, and to compile it with different heaps via templates. If you wonder why the Boost Graph Library does not use boost::heap, well, the reason is that BGL was implemented a few years ago, while boost::heap appeared this year.
Benchmarking on Road Networks
Here is the point that maybe interests you the most: can we do better than well-reputed C++ graph libraries?
I have tried three graph libraries: Boost Graph Library (BGL) v1.51, COIN-OR Lemon v1.2.3, and Google OrTools cheked out from svn on Sep 7th, 2012. They all have a Dijkstra implementation, even if I don’t know the implementation details. As a plus, the three libraries have python wrappers (but I have not test it). The BGL is a header only library. Lemon came after BGL. BGL, Lemon, and my implementation use (different) Fibonacci Heaps, while I have not clear what type of priority queue is used by OrTools.
Disclaimer: Google OrTools is much more than a graph library: among others, it has a Constraint Programming solver with very nice features for Large Neighborhood Search; however, we are interested here only in its Dijkstra implementation. Constraint Programming will be the subject of another future post.
A few tests on instances taken from the last DIMACS challenge on Shortest Path problems show the pros and cons of each implementation. Three instances are generated using the rand graph generator, while 10 instances are road networks. The test are done on my late 2008 MacBookPro using the apple gcc-4.2 compiler. All the source code, scripts, and even this post text, are available on github.
RAND Graphs
The first test compares the four implementations on 3 graphs with different density d that is the ratio . The graphs are:
- Rand 1: with n=10000, m=100000, d=0.001
- Rand 2: with n=10000, m=1000000, d=0.01
- Rand 3: with n=10000, m=10000000, d=0.1
For each graph, 50 queries between different pairs of source and destination nodes are performed. The table below reports the average of query times (total time divided by query numbers). The entries in bold highlight the shortest time per row.
| Graph | MyGraph | BGL | Lemon | OrTools |
|---|---|---|---|---|
| Rand 1 | 0.0052 | 0.0059 | 0.0074 | 1.2722 |
| Rand 2 | 0.0134 | 0.0535 | 0.0706 | 1.6128 |
| Rand 3 | 0.0705 | 0.5276 | 0.7247 | 4.2535 |
In these tests, it looks like my implementation is the winner… wow!
Although, the true winner is the boost::heap library, since the nasty implementation details
are delegated to that library.
… but come on! These are artificial graphs: who is really interested in shortest paths on random graphs?
Road Networks
The second test uses road networks that are very sparse graphs. We report only average computation time in seconds over 50 different pair of source-destination nodes. We decided to leave out OrTools since it is not very performing on very sparse graphs.
This table below shows the average query time for the standard implementations that use Fibonacci Heaps.
| Area | nodes | arcs | MyGraph | BGL | Lemon |
|---|---|---|---|---|---|
| Western USA | 6,262,104 | 15,248,146 | 2.7215 | 2.7804 | 3.8181 |
| Eastern USA | 3,598,623 | 8,778,114 | 1.9425 | 1.4255 | 2.7147 |
| Great Lakes | 2,758,119 | 6,885,658 | 0.1808 | 0.8946 | 0.2602 |
| California and Nevada | 1,890,815 | 4,657,742 | 0.5078 | 0.5808 | 0.7083 |
| Northeast USA | 1,524,453 | 3,897,636 | 0.6061 | 0.5662 | 0.8335 |
| Northwest USA | 1,207,945 | 2,840,208 | 0.3652 | 0.3506 | 0.5152 |
| Florida | 1,070,376 | 2,712,798 | 0.1141 | 0.2753 | 0.1574 |
| Colorado | 435,666 | 1,057,066 | 0.1423 | 0.1117 | 0.1965 |
| San Francisco Bay | 321,270 | 800,172 | 0.1721 | 0.0836 | 0.2399 |
| New York City | 264,346 | 733,846 | 0.0121 | 0.0677 | 0.0176 |
From this table, BGL and my implementation are equally good, while Lemon comes after.
What would happen if we use a diffent type of heap?
This second table shows the average query time for the Lemon graph library with a specialized Binary Heap implementation, and my own implementation with generic 2-Heap and 3-Heap (binary and ternary heaps) and with a Skew Heap. Note that in order to use a different heap I just modify a single line of code.
| Area | nodes | arcs | 2-Heap | 3-Heap | Skew Heap | Lemon 2-Heap |
|---|---|---|---|---|---|---|
| Western USA | 6,262,104 | 15,248,146 | 1.977 | 1.934 | 2.104 | 1.359 |
| Eastern USA | 3,598,623 | 8,778,114 | 1.406 | 1.372 | 1.492 | 0.938 |
| Great Lakes | 2,758,119 | 6,885,658 | 0.132 | 0.130 | 0.135 | 0.109 |
| California and Nevada | 1,890,815 | 4,657,742 | 0.361 | 0.353 | 0.372 | 0.241 |
| Northeast USA | 1,524,453 | 3,897,636 | 0.433 | 0.421 | 0.457 | 0.287 |
| Northwest USA | 1,207,945 | 2,840,208 | 0.257 | 0.252 | 0.256 | 0.166 |
| Florida | 1,070,376 | 2,712,798 | 0.083 | 0.081 | 0.080 | 0.059 |
| Colorado | 435,666 | 1,057,066 | 0.100 | 0.098 | 0.100 | 0.064 |
| San Francisco Bay | 321,270 | 800,172 | 0.121 | 0.117 | 0.122 | 0.075 |
| New York City | 264,346 | 733,846 | 0.009 | 0.009 | 0.009 | 0.007 |
Mmmm… I am no longer the winner: COIN-OR Lemon is!
This is likely due to the specialized binary heap implementation of the Lemon library. Instead, the boost::heap library has a d-ary-heap, that for d=2 is a generic binary heap.
So what?
Dijkstra’s algorithm is so beatiful because it has the elegance of simplicity.
Using an existing efficient heap data structure, it is easy to implement an “efficient” version of the algorithm.
However, if you have spare time, or you need to solve shortest path problems on a specific type of graphs (e.g., road networks), you might give a try with existing graph libraries, before investing developing time in your own implementation. In addition, be sure to read [4] and the references therein contained.
All the code I have used to write this post is available on github. If you have any comment or criticism, do not hesitate to comment below.
References
-
Pohl, I. Bi-directional and heuristic search in path problems. Department of Computer Science, Stanford University, 1969. [pdf]
-
Sniedovich, M. Dijkstra’s algorithm revisited: the dynamic programming connexion. Control and cybernetics vol. 35(3), pages 599-620, 2006. [pdf]
-
Dantzig, G.B. Linear Programming and Extensions. Princeton University Press, Princeton, NJ, 1962.
-
Delling, D. and Sanders, P. and Schultes, D. and Wagner, D. Engineering route planning algorithms. Algorithmics of large and complex networks Lecture Notes in Computer Science, Volume 5515, pages 117-139, 2009. [doi]
-
Goldberg, A.V. and Harrelson, C. Computing the shortest path: A-star search meets graph theory. Proc. of the sixteenth annual ACM-SIAM symposium on Discrete algorithms, 156-165, 2005. [pdf]