Floyd–Warshall algorithm

Kokshetau and Petropavl are 163 kilometres apart. Yet the shortest path between them by air is 2628 kilometres: from Kokshetau to Almaty and from there back north. The database has no direct flights between them, and the plane goes round half the map. 420 such answers are needed, one for each ordered pair of the 21 airports, and together they make up a table. Below, the Floyd–Warshall algorithm fills it in: at each step it allows transfers at one more airport and goes over the table row by row. Tap a cell to see the route, and an airport on the map to close it.

At first the table has only direct flights, 108 pairs out of 420. In the “hubs first” order the first step allows transfers at Almaty, and 382 pairs get a route; after Astana, all of them. The table becomes final after the seventh step already: from then on the airports improve nothing. Switch the order to “hubs last” and for thirteen steps in a row not a single cell changes, and the last two decide everything. The table at the end is the same, but the way to it differs: the strip at the bottom shows how many cells are already final after each step.

The routes come from the open database of Virtual Radar Server. It is filled by users of this program who receive signals from aircraft: for each callsign it records where the flight goes from and to. This is not a timetable, and some flights may be seasonal or out of date, so the table answers the question “how to fly by the flights of this database”, and the length of a route is the sum of great-circle distances. On the ground, of course, those 163 kilometres are driven without any Almaty.

The idea

The shortest route from i to j is either a direct flight or a chain with transfers. Going over chains is pointless; there are too many of them. Floyd–Warshall goes over something else: not routes but the set of airports where transfers are allowed. Number the airports and let d_k[i][j] be the length of the shortest route from i to j all of whose transfers are at airports numbered less than k. With k = 0 no transfers are allowed anywhere, and d₀ is the table of direct flights. When all are allowed, d_n is the answer.

Allow transfers at one more airport k. The best route from i to j either does not go through k, and then it is already in the table, or goes through k exactly once, and then it is made of the best route from i to k and the best from k to j, both already in the table. So the new table is obtained from the old one by a single formula for each cell.

This is dynamic programming in its pure form: the subproblem is “a route with transfers only at the first k airports”, the transition is one new airport. In code k is the outer loop, and i and j run over the whole table:

for k in range(n):
    for i in range(n):
        for j in range(n):
            d[i][j] = min(d[i][j], d[i][k] + d[k][j])

In the scene this is seen literally. The row and column of airport k are highlighted: at step k each row i takes one number from its column k, how far it is to fly to k, adds it to the whole of row k and keeps the smaller of the two in each cell.

Why the table is right

There remains a question interviewers like to ask: the formula speaks of two tables, d_k and d_k+1, while the code writes into one. Won’t a new number in d[i][k] spoil the next cells of the row? It won’t, because row k and column k do not change at step k: d[i][k] = min(d[i][k], d[i][k] + d[k][k]), and d[k][k] = 0. All the cells of a step read the same numbers from row and column k, in whatever order they are visited. So the memory needed is O(n²), not O(n³).

The trap: the order of the loops

The three loops look symmetrical, and the hand writes them in the habitual order by itself: for i, for j, for k. Then each cell (i, j) is computed once, by going over all transfers k, which seems even more natural. But the proof above relied on the cells d[i][k] and d[k][j] being final for all earlier transfers by the time of step k. In the order i, j, k this is not so. Here are six airports of the west, numbered as in the data, by ICAO code. The scene starts from the end of the first pass of the wrong loop.

Turkistan’s row is computed first. When the turn comes to the cell “Turkistan → Aktau”, the cell “Turkistan → Aktobe” is still empty, since its turn in the same row comes later, and Kostanay’s row has not been computed at all. And the real route goes just through Kostanay and Aktobe: 1160 + 547 + 847 = 2554 km. The pass ends, and three cells remain empty although the routes exist. The right order needs 27 comparisons for the same table, counting only those where both halves of the route are already known, while the wrong one is mistaken after 31 comparisons, and only a second pass fixes everything.

On the whole network the error is quieter: the hubs save almost everything, and the wrong loop is mistaken in two cells out of 420: from Turkistan to Aktau it counted 2604 km via Astana instead of 2554 via Kostanay and Aktobe. The order i, k, j is mistaken in seven cells. That is exactly why the error is easy to miss: on a small test the answer agrees, on a large one it almost agrees.

The wrong loop can be rescued by repetition. In 2019 Ikumi Hide, Soh Kumabe and Takanori Maehara proved that the orders i, j, k and i, k, j repeated three times in a row always give the right table; the paper is called just that: “Incorrect implementations of the Floyd–Warshall algorithm give correct solutions after three repeats”. Two are sometimes not enough: the path 1 → 5 → 3 → 2 → 6 → 4 → 0 from the listing still does not know the way from 1 to 0 after the second pass. But three times the work for the sake of an error in one line is a bad trade. The rule is easier to remember: k goes outside.

Code

Besides the distances the program keeps a second table, nxt[i][j], the first flight of the best route from i to j. When a route through k turns out shorter, its first flight is the same as that of the route from i to k, so nxt[i][j] = nxt[i][k]. To recover the whole route, follow these links from i until you reach j. There is one table for all pairs, and the path between any two airports is read in as many steps as it has flights.

108 pairs have a direct flight, 299 manage with one transfer, 13 with two: there are never more than three flights. Detours can be huge: from Astana to Karaganda, 200 km away, the shortest route by the flights of the database goes through Almaty and stretches to 1708 km. The shortest of the farthest paths belongs to Shymkent: from it to any airport is at most 1935 km. Such a vertex is called the centre of the graph, and it is neither Almaty nor Astana: Almaty’s farthest path is 2097 km, Kokshetau’s 3023 km.

The last three lines are the scene in numbers: how many cells are final after each step in the three orders. In Python the inner loop is sped up by two techniques: the rows d[i] and d[k] are put in local variables, and a row where d[i][k] is infinite is skipped entirely, since there is no flying through k from it anyway. Faster still is to hand the whole step k to numpy, as in the course chapter: np.minimum(d, d[:, k, None] + d[None, k, :], out=d).

The same loop without numbers: who depends on whom

Replace the numbers with “yes” and “no”: reach[i][j] says whether j can be reached from i. The minimum turns into “or”, the sum into “and”, and the step formula reads: j can be reached from i if it could be before, or if k can be reached from i and j from k. This is Warshall’s algorithm of 1962, which computes the transitive closure. A row of the table is conveniently stored as bits: then “if there is a path from i to k, add the whole of row k to row i” is a single OR over machine words, 64 cells at a time.

An example is the Go standard library, the same graph as in the article on topological sorting. There we were interested in the build order; here in a simpler and more frequent question: what does a package need at all, not just directly.

The fmt package needs 58 packages, net/http 186, and the command go list -deps prints exactly as many plus the package itself. 331 packages out of 379 depend on unsafe. The order of k again changes only the way: going from lower packages to upper ones, by the fiftieth step 12,218 pairs out of 21,504 are known, and in the reverse order 2980.

Other pairs of operations fit as well. Maximum instead of minimum and minimum instead of sum give the “widest” path, the route whose narrowest point is as wide as possible. With operations on regular expressions, union, concatenation and star, the same triple loop builds a regular expression from a finite automaton: this method is called Kleene’s algorithm and traced back to Stephen Kleene’s work of 1956. The general framework for all the variants is semirings: one needs an operation “choose among options” and an operation “continue the path” that agree with each other.

Negative edges and cycles

The proof required one thing: cycles are non-negative. Negative edges do not get in its way, and this is the difference from Dijkstra. If there is a negative cycle, shortest paths through it do not exist, and the algorithm notices this by itself: a negative number appears on the diagonal, a path from i back to i shorter than zero. So after the algorithm it is enough to check d[i][i] < 0. Such a cycle itself and its meaning, a round of exchanges that returns more than there was, are taken apart in the article on Bellman–Ford. With a negative cycle the numbers in the table can grow in absolute value very quickly, and with fixed-width integers they are best cut off as soon as the diagonal goes negative.

What it costs

Three loops over n are O(n³) time, the table O(n²) memory, and the number of edges affects neither. The alternative is Dijkstra from each vertex: O(n (n + m) log n) with a heap. For a sparse graph, where m is of order n, this is O(n² log n), noticeably better. For a dense one, where m is of order n², it is O(n³ log n), worse. Dijkstra cannot handle negative edges, but they can be removed with potentials: that is Johnson’s algorithm. A benchmark on random graphs: four edges per vertex (a sparse graph) and n edges per vertex (a dense one).

The slope of Floyd–Warshall in Go is close to 3: twice the vertices, eight times the time; in Python the slope in the benchmark is lower, about 2.3. In Python on a sparse graph of 200 vertices Dijkstra from each vertex is almost 10 times faster, and on a dense one of 100 vertices they are nearly level: Floyd–Warshall is only 1.4 times slower. In Go the boundary shifts: Floyd–Warshall has no heap and no pointers, only a pass along rows in order, and even on a sparse graph it is faster up to four hundred vertices, and on a dense one of 400 vertices 6 times faster. The numpy version on the same 200 vertices is 42 times faster than pure Python: step k is one operation over the whole table.

Hence the rule of choice. If all pairs are needed and there are hundreds or a couple of thousand vertices, Floyd–Warshall: it is shorter, does not go wrong on negative edges and suits vector operations and the cache well. With tens of thousands of vertices and more, the n² table by itself will not fit in memory, and pairs are needed on demand, that is, Dijkstra, A* or precomputation such as contraction hierarchies. For the streets of Almaty with 7909 intersections the table would take 62.5 million cells.

Where it works

• Libraries. SciPy has scipy.sparse.csgraph.floyd_warshall, NetworkX has floyd_warshall and floyd_warshall_numpy: when the whole distance table of a small graph is needed, these are used.
• Closures of relations. Who depends on whom, who reports to whom, which states are reachable: wherever a relation is given by a matrix and its closure is needed. The bit version handles 64 pairs per operation.
• Automata and regular expressions. Kleene’s algorithm, by which a regular expression is built from an automaton, is the same triple loop over “allowed intermediate states”.
• Real-time systems. In model checking with timers, constraints of the form “x − y ≤ c” are kept in a matrix, and to bring it to canonical form with the tightest constraints, it is gone over with Floyd–Warshall.

Three papers of 1959–1962

In 1959 Bernard Roy published a note “Transitivité et connexité” in the proceedings of the Academy of Sciences, with the same way of computing the closure. In January 1962 Stephen Warshall’s “A theorem on Boolean matrices” came out in the Journal of the ACM. The story goes that Warshall bet a colleague a bottle of rum on who would be first to find out whether this algorithm always works, and brought a proof the next morning; the rum, it is said, was shared with the loser. In June of the same year Robert Floyd published “Algorithm 97: Shortest Path” in Communications of the ACM, one page with the same scheme for shortest paths. The same year Peter Ingerman, in “Algorithm 141: Path Matrix”, wrote the closure as three nested loops, in the form in which it is written today.

Hence the variety of names: the algorithm of Floyd, of Warshall, of Roy–Warshall, of Roy–Floyd. And the idea itself, going over not paths but the set of allowed intermediate vertices, goes back to Stephen Kleene’s 1956 work on events in nerve nets and finite automata, where regular expressions appeared.

Problems

In the first problem airports open one at a time, and after each the sum of all shortest distances is needed: recomputing the table from scratch takes too long, and the order of k can be chosen freely. In the second the shortest loop along streets is needed, and it has to be checked at the moment of step k when vertex k does not yet take part in paths.