Floyd-Warshall, in depth: the all-pairs recurrence, path reconstruction, negative cycles, semirings, blocked versions and incremental updates

The recurrence

Number the vertices 0 to n-1. Define d_k(i, j) as the length of the shortest path from i to j whose intermediate vertices all come from the set {0, ..., k-1}. With k = 0 no intermediates are allowed, so d₀(i, j) is the direct edge weight, 0 on the diagonal and infinity where there is no edge. With k = n every vertex is allowed, so d_n is the answer.

Going from k to k+1 adds one new allowed vertex, k. A shortest path from i to j using intermediates in {0, ..., k} either avoids k, in which case its length is d_k(i, j), or passes through k exactly once, in which case it splits into a path i to k and a path k to j, each using only {0, ..., k-1}. Without negative cycles a shortest path never needs to visit k twice. So:

d[k+1][i][j] = min( d[k][i][j],  d[k][i][k] + d[k][k][j] )

Two consequences shape the implementation. First, k must be the outermost loop: the whole matrix for stage k has to be finished before stage k+1 reads it. Putting k innermost is the most common bug and gives wrong answers on many graphs. Second, one matrix is enough. During stage k, row k and column k do not change, because d(i, k) can only improve through k itself, and d(i, k) + d(k, k) is not smaller than d(i, k) when d(k, k) is 0 or more. So updating in place reads the same values the three-dimensional version would.

Implementation with path reconstruction

The implementation keeps a nxt matrix alongside distances: nxt[i][j] is the first hop on the best known path from i to j. When the path through k wins, the first hop toward j becomes the first hop toward k. Reconstruction then follows first hops.

INF = float("inf")

def floyd_warshall(n, edges):
    d = [[INF] * n for _ in range(n)]
    nxt = [[None] * n for _ in range(n)]
    for i in range(n):
        d[i][i] = 0
        nxt[i][i] = i
    for u, v, w in edges:
        if w < d[u][v]:             # keep the cheapest parallel edge
            d[u][v] = w
            nxt[u][v] = v
    for k in range(n):              # k MUST be outermost
        dk = d[k]
        for i in range(n):
            dik = d[i][k]
            if dik == INF:
                continue            # skip rows that cannot reach k
            di = d[i]
            for j in range(n):
                nd = dik + dk[j]
                if nd < di[j]:
                    di[j] = nd
                    nxt[i][j] = nxt[i][k]
    return d, nxt

def path(nxt, u, v):
    if nxt[u][v] is None:
        return []                   # unreachable
    p = [u]
    while u != v:
        u = nxt[u][v]
        p.append(u)
    return p

Skipping rows where d[i][k] is infinite is more than an optimisation on sparse graphs. It also avoids computing infinity plus a negative number, which is still infinity in floating point but becomes a wrong finite number if you replace infinity with a large integer.

Worked example

Take the graph in the figure, with edges 0 to 1 (3), 0 to 3 (7), 1 to 0 (8), 1 to 2 (2), 2 to 0 (5), 2 to 3 (1) and 3 to 0 (2). Running the code above and printing the matrix after each stage gives the table below; changed entries from the previous stage are in bold.

Follow one entry. d(1, 0) starts at 8, the direct edge. At k=2 it drops to 7 via 1 to 2 to 0 (2 + 5). At k=3 it drops to 5 via 1 to 2 to 3 to 0 (2 + 1 + 2), because vertex 3 is now an allowed intermediate and d(1, 3) had already become 3 at k=2. Path reconstruction confirms it: path(nxt, 1, 0) returns [1, 2, 3, 0], path(nxt, 0, 3) returns [0, 1, 2, 3] and path(nxt, 3, 2) returns [3, 0, 1, 2].

Negative cycles

If the graph has a cycle of negative total weight, "shortest path" is undefined for any pair whose path can touch the cycle, because going around once more always helps. Floyd-Warshall detects this cheaply: after the main loop, vertex v lies on or can reach and return through a negative cycle exactly when d(v, v) is below 0. On the three-vertex cycle 0 to 1 (1), 1 to 2 (-3), 2 to 0 (1), with total weight -1, the code ends with diagonal values -1, -1 and -2. Vertex 2 shows -2 because its row was improved twice, which also shows that once a negative cycle exists the numbers are not meaningful distances, only evidence.

Checking the diagonal is not enough if you need to know which pairs are affected. Pair (i, j) has no shortest path if some k with d(k, k) below 0 is reachable from i and can reach j. Mark those pairs as minus infinity:

def mark_negative(d):
    n = len(d)
    for k in range(n):
        if d[k][k] < 0:
            for i in range(n):
                if d[i][k] == INF:
                    continue
                for j in range(n):
                    if d[k][j] != INF:
                        d[i][j] = -INF
    return d

On a five-vertex graph that adds an edge from the cycle to vertex 3 and a separate vertex 4 pointing at 3, this marks every pair among 0, 1, 2 and 3 that starts in the cycle as minus infinity, leaves d(4, 3) = 1 untouched, and leaves d(3, 3) = 0, since 3 cannot get back into the cycle. For the single-source version and cycle extraction, see Detecting Negative Cycles.

The semiring view and transitive closure

Nothing in the recurrence depends on min and plus specifically. Replace them with any pair of operations that form a closed semiring and the same triple loop solves a different problem:

Transitive closure is worth doing with bitsets. Store each row as one integer whose bit j means i reaches j; then the inner loop over j collapses to a single OR of whole rows, and the algorithm runs in n³/w word operations for word size w.

def closure(n, edges):
    reach = [1 << i for i in range(n)]      # every vertex reaches itself
    for u, v in edges:
        reach[u] |= 1 << v
    for k in range(n):
        bit = 1 << k
        for i in range(n):
            if reach[i] & bit:                # i reaches k ...
                reach[i] |= reach[k]          # ... so i reaches all k reaches
    return reach

The same loop with min-plus over matrices is also the bridge to algebraic methods: repeated min-plus squaring of the adjacency matrix computes all-pairs distances in O(n³ log n), slower here, but the idea behind Matrix Exponentiation for bounded-hop paths.

Making it fast

The cost is n³ relaxations and n² memory. For n = 1,000 that is a billion relaxations: minutes in pure Python, around a second in tight C. For n = 10,000 it is a trillion relaxations and an 800 MB float64 matrix. Two techniques help.

Vectorise each stage. Stage k is an outer sum of column k and row k followed by an element-wise minimum, which NumPy does in one call. This removes the Python overhead and gives a large constant-factor speed-up while keeping the memory at one matrix plus one temporary.

import numpy as np

def fw_numpy(w):                         # w: n x n, np.inf for no edge, 0 on diagonal
    d = np.array(w, dtype=np.float64)
    for k in range(d.shape[0]):
        np.minimum(d, d[:, k, None] + d[None, k, :], out=d)
    return d

Block it for cache and GPUs. Split the matrix into B by B tiles. For each block index kb, first run Floyd-Warshall inside the diagonal tile (kb, kb); then update the tiles in block row kb and block column kb using that tile; then update every remaining tile (i, j) from tiles (i, kb) and (kb, j). The third phase is a min-plus matrix multiply over tiles that fit in cache or GPU shared memory, and it does almost all the work. This is the standard form for GPU implementations and makes the algorithm compute-bound rather than memory-bound.

When to use something else

For the sparse cases, follow Johnson's Algorithm and Dijkstra, in depth; the reweighting trick in Johnson's is the standard way to get Dijkstra speed with negative edges.

Failure modes

• k not outermost. The code still runs and often passes small tests. Write a randomised test against repeated Bellman-Ford or Dijkstra; it catches this immediately.
• Integer infinity overflow. With a sentinel such as 2³¹-1, INF + INF overflows to a negative number in fixed-width integers, inventing a path. Skip infinite terms explicitly, or use a sentinel at most half the maximum value.
• Negative cycles growing values. Once a negative cycle exists, values can keep falling and in integer arithmetic can underflow. Check the diagonal and stop trusting results.
• Undirected graphs with a negative edge. An undirected negative edge is a two-edge negative cycle, so every pair that touches it is undefined. Floyd-Warshall will tell you, but the model is usually what is wrong.
• Parallel edges and self-loops. Keep the minimum parallel edge; a negative self-loop is a negative cycle; a positive self-loop must not overwrite the 0 on the diagonal.
• Floating-point ties. Accumulated rounding can flip the chosen path between equal-cost routes; compare with a tolerance if path identity matters, or use integers.

Operational guidance

In production, all-pairs results are usually precomputed tables: hop costs between warehouse locations, latency between data centres, or a routing table in a small network. Three habits keep them correct and cheap.

First, update incrementally when an edge gets cheaper or is added. A new edge u to v with weight w can only improve paths that use it, so one O(n²) pass suffices: d(i, j) = min(d(i, j), d(i, u) + w + d(v, j)). Copy column u and row v first so the pass does not read values it has just written. This holds only if the new edge does not close a negative cycle, that is unless d(v, u) + w is below 0; check that first. On the worked example, adding 3 to 2 with weight 1 changes only d(3, 2), from 7 to 1, and matches a full recompute. Increases and deletions are harder; recompute, or recompute only rows whose stored path used the edge.

def decrease_edge(d, u, v, w):
    n = len(d)
    if d[v][u] + w < 0:
        raise ValueError("new edge closes a negative cycle; recompute and mark")
    du = [d[i][u] for i in range(n)]   # snapshot column u
    dv = list(d[v])                    # snapshot row v
    for i in range(n):
        if du[i] == INF:
            continue
        for j in range(n):
            c = du[i] + w + dv[j]
            if c < d[i][j]:
                d[i][j] = c

Second, validate every recompute: diagonal all zero, triangle inequality holds on a random sample of triples, and a few pairs agree with an independent single-source run such as Bellman-Ford. Third, store the nxt matrix with the distances, since a distance without a route is rarely actionable.

What to do next

1. Implement the version above with path reconstruction and run it on the worked example.
2. Write a randomised test against Bellman-Ford from every source, including negative edges.
3. Add negative-cycle marking and test it on a graph with a reachable and an unreachable cycle.
4. Rewrite the stage loop in NumPy and time n = 500, 1000 and 2000.
5. Implement the bitset transitive closure and compare memory with the distance matrix.
6. Add the O(n squared) incremental update and check it against full recompute.
7. For your own graph, compute n cubed against n m log n and pick Floyd-Warshall, Johnson or Dijkstra.