All-Pairs Shortest Paths, in depth: choosing between BFS, repeated Dijkstra, Johnson, Floyd-Warshall, min-plus squaring and Seidel

The problem and the size of its answer

Take a directed graph with n vertices, m edges and a weight on each edge. The output is an n by n matrix D, where D[u][v] is the minimum total weight of any path from u to v. It is infinity if v is unreachable, and undefined if a negative cycle can be reached on the way. If you need the paths themselves as well as their lengths, you also keep a second n by n matrix, either a predecessor matrix or a next-hop matrix. A path is then rebuilt by walking it one step at a time.

That output size is the first fact to absorb. Whatever the algorithm, the answer has n squared entries. At n = 10,000 that is 10⁸ distances, which is 400 MB as 32-bit integers, and the next-hop matrix doubles it. At n = 1,000,000 it is 10¹² entries, about 4 TB, which no single machine stores comfortably. So for graphs of a million vertices, APSP is a storage problem before it is an algorithm problem.

The family at a glance

Two facts drive most choices. First, n single-source runs are independent of each other, so BFS, Dijkstra and Johnson parallelise across sources with no coordination. Second, when m is close to n squared, repeated Dijkstra loses its advantage and pays heap overhead on top. Floyd-Warshall's tight triple loop then wins, and it is far simpler to block for cache and GPU.

Worked cost example: 10,000 vertices, 50,000 edges

Run the arithmetic for one concrete graph before reaching for any algorithm. Take n = 10,000 vertices and m = 50,000 edges, an average out-degree of 5, which is typical of road segments, dependency graphs and social subgraphs.

Floyd-Warshall performs over a hundred times more basic operations here. Its inner loop is a branch-free, vectorisable min over contiguous memory, while Dijkstra's heap is pointer-chasing and cache-hostile, so the wall-clock gap is smaller than the operation ratio. It does not vanish at this sparsity, though. The counts are equal when n(n + m) log₂ n = n³, that is, near m = n²/log₂ n, about 7.5 x 10⁶ edges here, or 7.5% density. At half density, Floyd-Warshall does about six times fewer operations. That is the regime where Floyd-Warshall is the right tool. Treat these as orders of magnitude and measure on your own graph and hardware before committing.

Repeated Dijkstra, the sparse-graph default

For sparse graphs with non-negative weights, the practical answer is one Dijkstra run per source. Keep a parent row per source so paths can be rebuilt. The version below uses lazy deletion: stale heap entries are skipped rather than decreased in place, which is what Python's heapq supports.

import heapq
INF = float("inf")

def dijkstra(adj, s, n):
    """adj[u] is a list of (v, w) with w >= 0. Returns dist and parent rows."""
    dist = [INF] * n
    parent = [-1] * n
    dist[s] = 0
    pq = [(0, s)]
    while pq:
        d, u = heapq.heappop(pq)
        if d > dist[u]:          # stale entry: u was already settled cheaper
            continue
        for v, w in adj[u]:
            nd = d + w
            if nd < dist[v]:
                dist[v] = nd
                parent[v] = u
                heapq.heappush(pq, (nd, v))
    return dist, parent

def apsp(adj, n):
    D, P = [], []
    for s in range(n):           # independent: farm out across cores or machines
        d, par = dijkstra(adj, s, n)
        D.append(d)
        P.append(par)
    return D, P

def path(P, s, t):
    if P[s][t] == -1 and s != t:
        return None              # unreachable
    out = [t]
    while out[-1] != s:
        out.append(P[s][out[-1]])
    return out[::-1]

Because the source loop has no shared state, the production version is a process pool or a distributed map over source ranges, with each worker writing its rows to disk. If edges can be negative, run Bellman-Ford once from a virtual vertex to get potentials h. Reweight each edge as w(u, v) + h(u) - h(v), which is non-negative, and the same loop then applies. That is Johnson's algorithm. Dijkstra with potentials explains why the reweighting preserves shortest paths. Remember to undo it at the end: D[u][v] = D'[u][v] - h(u) + h(v).

Min-plus products and repeated squaring

There is an algebraic view of the problem. Replace multiplication with addition and addition with min, and the matrix product C = A (x) B becomes C[i][j] = min over k of A[i][k] + B[k][j]. Call the weight matrix W, with 0 on the diagonal and infinity for missing edges. Then W (x) W gives shortest paths that use at most two edges, and squaring repeatedly doubles the hop limit each time. Shortest paths without negative cycles use at most n - 1 edges, so ceil(log₂(n - 1)) squarings give the full answer.

Here is a worked example on four vertices with edges A to B (3), B to C (2), C to D (1), A to D (10) and D to A (4). In W, the direct edge gives A to D a distance of 10. After one squaring, W² covers paths of at most two edges. A to C becomes 5 (through B), B to D becomes 3, C to A becomes 5 and D to B becomes 7. A to D is still 10, because no two-edge path reaches D from A. After the second squaring, W⁴ covers paths of up to four edges. A to D drops to 6 (A to C at 5, plus C to D at 1), B to A becomes 7, C to B becomes 8 and D to C becomes 9. With n = 4, paths need at most 3 edges, so two squarings are enough.

import numpy as np

def min_plus(A, B):
    n = A.shape[0]
    C = np.empty_like(A)
    for i in range(n):                       # row at a time keeps memory at O(n^2)
        C[i] = np.min(A[i][:, None] + B, axis=0)
    return C

def apsp_power(W, k=None):
    """Shortest paths using at most k edges; k = n - 1 gives full APSP.
    W: float array, 0 on diagonal, np.inf for no edge."""
    k = W.shape[0] - 1 if k is None else k
    R = np.full_like(W, np.inf)
    np.fill_diagonal(R, 0.0)                 # (min, +) identity: zero edges
    B = W.copy()
    while k:                                 # binary exponentiation of W
        if k & 1:
            R = min_plus(R, B)
        B = min_plus(B, B)
        k >>= 1
    return R

At O(n³ log n) this is strictly slower than Floyd-Warshall, so why keep it? First, it answers hop-limited questions exactly: pass k to get the cheapest route with at most k edges, such as k transfers or k relays. Binary exponentiation matters here, because plain repeated squaring overshoots to the next power of two and cannot take extra hops back. Second, the semiring framing carries over. Swap (min, +) for (max, min) to get widest paths, or for (or, and) to get transitive closure, and the same code answers bottleneck-bandwidth and reachability questions.

Seidel's algorithm for unweighted undirected graphs

For undirected, unweighted, connected graphs, Seidel's algorithm turns APSP into O(log n) ordinary integer matrix products. It squares the graph: two vertices are adjacent in the new graph if they are within distance 2 in the old one. It solves the squared graph recursively, then recovers each distance from the halved one. Every true distance is either 2T or 2T - 1, where T is the distance in the squared graph, and a degree comparison decides which.

import numpy as np

def seidel(A):
    """A: symmetric 0/1 int matrix, zero diagonal, graph connected."""
    n = A.shape[0]
    Z = A @ A
    B = ((A == 1) | (Z > 0)).astype(np.int64)
    np.fill_diagonal(B, 0)
    if (B + np.eye(n, dtype=np.int64)).min() == 1:   # squared graph is complete
        return 2 * B - A
    T = seidel(B)
    X = T @ A
    deg = A.sum(axis=0)
    return 2 * T - (X < T * deg[np.newaxis, :]).astype(np.int64)

Its value is that it rides on matrix multiplication, the most heavily optimised kernel on any CPU or GPU. One practical tip: integer matmul in NumPy does not use BLAS. Casting to float64 does, and stays exact here because every entry is at most n squared, far below 2⁵³. The catch is the precondition. A disconnected graph makes the base case unreachable, so run it per connected component.

When not to compute all pairs

Many systems that say they need all pairs really need many pairs, chosen at query time. For those, keep the graph and pay per query.

• On-demand single-source plus cache. Run Dijkstra or BFS from each source the first time it is queried and cache the row in an LRU. Query traffic usually concentrates on a small set of hot sources, so a few thousand cached rows can serve most requests.
• Goal-directed search. A* with landmark lower bounds (ALT) precomputes distances from a few dozen landmark vertices. The triangle inequality then gives admissible estimates that prune most of the search. Road-routing engines go further with contraction hierarchies.
• Approximate distance oracles. Thorup and Zwick's construction answers queries with stretch at most 2k - 1, using O(k n^{1 + 1/k}) space. Paying a factor of 3 in accuracy (k = 2) cuts storage from n² to about n^1.5.

Operational guidance

• Store compactly. Use the narrowest integer type that holds the maximum distance, with a sentinel for unreachable. Store next-hop indices as int16 when n is under 32,768. Write rows in source blocks so a reader can memory-map one block.
• Shard by source. Row ranges are the natural unit of parallelism, checkpointing and incremental recomputation.
• Incremental updates. When one edge weight drops, any pair whose path could now go through that edge may improve, and checking every pair against the new edge costs O(n²). Weight increases and deletions are harder. Recompute only the rows whose shortest-path trees used the edge, which you can find from the parent matrix.
• Validate cheaply. For a random sample of sources, recompute the row with an independent implementation and compare. Also check the triangle inequality D[u][v] <= D[u][x] + w(x, v) on sampled edges.
• Prefer integer weights when results must be reproducible; float rounding can resolve ties differently between runs.

Failure modes

• Undetected negative cycles. Dijkstra silently returns wrong answers on negative edges. Floyd-Warshall reveals a reachable negative cycle as a negative diagonal entry, and Johnson's Bellman-Ford phase reports it. Negative-cycle detection covers the details.
• Infinity arithmetic. With integer sentinels, INF + w overflows and wraps negative, producing impossible short paths. Skip the relaxation when either operand is the sentinel, or use a sentinel no larger than half the type's maximum.
• Memory blow-up. Broadcasting a min-plus product across all three indices allocates n³ elements, which at n = 2,000 is 8 x 10⁹ floats. Process one row at a time, as the code above does.
• Seidel on a disconnected graph recurses until the stack overflows, because the squared graph never becomes complete.
• Paths that do not match distances. If a predecessor is updated without its distance, or the other way round, the rebuilt path has a different length from D. Assert that the two agree on sampled pairs.

What to do next

1. Write down n, m, the weight type and the query pattern for your graph, and do the cost and memory arithmetic from this page before choosing.
2. If the weights are non-negative and the graph is sparse, implement the repeated Dijkstra code above and parallelise it across sources. Add Johnson's reweighting if you have negative edges.
3. If the graph is dense, read Floyd-Warshall in depth and use a blocked implementation.
4. Run the min-plus example by hand, then use apsp_power(W, k) to answer a hop-limited routing question on your data.
5. If n squared does not fit, prototype an on-demand plus cache design and measure the hit rate on real query logs before building an oracle.
6. Add the sampled cross-check and triangle-inequality validation to the pipeline that produces the matrix.