07 - Shortest Path Algorithms: Dijkstra, Bellman-Ford, Floyd-Warshall

In the previous article breadth-first search found the path with the fewest edges, and a counterexample showed that it is not the cheapest: a direct link of 10 beats, by number of edges, a three-leg trip of 1 each. When edges have a cost (kilometres, minutes, euros of toll) the question becomes "what is the minimum-cost path?", and the answer depends on a single thing: are there negative weights or not. Dijkstra is fast but does not admit them; Bellman-Ford admits them and recognizes negative cycles; Floyd-Warshall answers for all pairs at once. Here you implement all three, compare them with networkx on thousands of graphs, and build the case in which Dijkstra gets it wrong. You need the adjacency lists and the BFS from article 6.

What is computed, and the gesture that unites them

Given a directed graph with a weight on every edge, the cost of a path is the sum of the weights of its edges, and the distance between two vertices is the minimum of these costs (infinite if there is no path). The gesture common to the three algorithms is the relaxation of an edge (u, v) of weight w: if reaching v through u costs less than the best estimate you have, that is dist[u] + w < dist[v], you update dist[v]. Estimates start from infinity (0 at the source) and only go down; the algorithms differ in the order in which they relax the edges, and it is that order that decides when the result is guaranteed.

The sp.py module, which opens the code of the whole article, contains the three functions. You will read them in pieces in the sections that follow; for now it is enough to know that each one counts its own operations (heap insertions, passes over the edges, steps of the inner loop) because complexity is stated with a count, not with seconds.

Dijkstra: always the nearest frontier, as long as weights are non-negative

Dijkstra keeps a heap of (distance, vertex) pairs and always pops the vertex with the smallest estimate: at that point its estimate is final, because any other path would go through vertices that are already farther away and, with non-negative weights, could only add cost. Then it relaxes the edges leaving that vertex. The comparison is a fire spreading at a speed equal to the weight of the edges: every vertex catches fire in the order of its true distance, and one that has caught fire does not catch fire again.

In Python the heap is heapq, and two details of the module matter. There is no "lower the priority" operation: you insert a new entry and the old ones stay, so on extraction you have to recognize and discard them, as the heapq documentation suggests with the marked-entry technique. Moreover pairs are compared element by element, so at equal distance the vertex decides: with numeric vertices that is fine, with objects that cannot be compared a counter would be needed, as in the example on the same page.

The dijkstra function in sp.py has two modes. With definitivo=True (the classic version) a vertex popped once is put in chiusi and is never expanded again; with False only the stale entries are discarded (d > dist[u]) and a vertex improved after its extraction is expanded again. The difference does not matter with non-negative weights and matters enormously with negative weights, as you will see shortly. First the correctness test: 3,000 random directed graphs with integer weights from 0 to 20, compared with networkx (3.7, in a trial virtual environment), for both modes. Then the counters on three large graphs, and finally the worst case.

import heapq

INF = float("inf")


def dijkstra(adj, s, definitivo=True):
    """adj[u] = list of (v, weight). Returns distances and operation counters.

    definitivo=True: a vertex popped from the heap is never reopened
    (the classic version, correct only with weights >= 0).
    """
    dist = [INF] * len(adj)
    dist[s] = 0
    heap = [(0, s)]
    chiusi = set()
    push, pop, rilassati = 1, 0, 0
    while heap:
        d, u = heapq.heappop(heap)
        pop += 1
        if definitivo:
            if u in chiusi:
                continue
            chiusi.add(u)
        elif d > dist[u]:  # stale entry: u has already been improved
            continue
        for v, w in adj[u]:
            if d + w < dist[v]:
                dist[v] = d + w
                rilassati += 1
                heapq.heappush(heap, (dist[v], v))
                push += 1
    return dist, {"push": push, "pop": pop, "relaxations": rilassati}


def bellman_ford(n, archi, s=None):
    """archi = list of (u, v, weight). If s is None it starts from 0 on every vertex
    (finds a negative cycle anywhere). Returns dist, negative_cycle, passes, checks."""
    dist = [INF] * n if s is not None else [0] * n
    if s is not None:
        dist[s] = 0
    controlli = 0
    for passata in range(1, n + 1):
        cambiato = False
        for u, v, w in archi:
            controlli += 1
            if dist[u] + w < dist[v]:
                dist[v] = dist[u] + w
                cambiato = True
        if not cambiato:
            return dist, False, passata, controlli
    return dist, True, n, controlli  # still changing at pass n: negative cycle


def floyd_warshall(n, peso):
    """peso[i][j] = weight of the edge or INF; diagonal 0. Returns matrix and steps."""
    d = [riga[:] for riga in peso]
    passi = 0
    for k in range(n):
        dk = d[k]
        for i in range(n):
            dik = d[i][k]
            di = d[i]
            for j in range(n):
                passi += 1
                if dik + dk[j] < di[j]:
                    di[j] = dik + dk[j]
    return d, passi

import random
import networkx as nx
from sp import dijkstra, INF

rng = random.Random(31)
for prova in range(3000):
    n = rng.randint(1, 30)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 3)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v, weight=rng.randint(0, 20))
    adj = [[(v, G[u][v]["weight"]) for v in G[u]] for u in range(n)]
    s = rng.randrange(n)
    ref = nx.single_source_dijkstra_path_length(G, s)
    for modo in (True, False):
        dist, _ = dijkstra(adj, s, definitivo=modo)
        assert {v: d for v, d in enumerate(dist) if d < INF} == ref
print("3000 digraphs with weights 0..20: both variants equal to networkx", nx.__version__)

print("V       E        push     pop      relaxations   push <= E+1")
for n in (1000, 10000, 100000):
    e = 4 * n
    adj = [[] for _ in range(n)]
    for _ in range(e):
        u, v = rng.randrange(n), rng.randrange(n)
        adj[u].append((v, rng.randint(1, 100)))
    _, c = dijkstra(adj, 0)
    print(f"{n:<7} {e:<8} {c['push']:<8} {c['pop']:<8} {c['relaxations']:<13} {c['push'] <= e + 1}")

print()
print("worst case: k vertices u_i, each linked to all m vertices v")
print("k=m   E        push     push/E")
for k in (30, 60, 120):
    n = 1 + k + k  # s, u_1..u_k, v_1..v_k
    adj = [[] for _ in range(n)]
    for i in range(1, k + 1):
        adj[0].append((i, i))  # dist(u_i) = i
        for j in range(k + 1, n):
            adj[i].append((j, 2 * k - 2 * i))  # dist(u_i) + weight = 2k - i: decreases with i
    e = sum(len(x) for x in adj)
    dist, c = dijkstra(adj, 0)
    print(f"{k:<5} {e:<8} {c['push']:<8} {c['push'] / e:.3f}")

# Output:
# 3000 digraphs with weights 0..20: both variants equal to networkx 3.7
# V       E        push     pop      relaxations   push <= E+1
# 1000    4000     1268     1268     1267          True
# 10000   40000    12878    12878    12877         True
# 100000  400000   129181   129181   129180        True
#
# worst case: k vertices u_i, each linked to all m vertices v
# k=m   E        push     push/E
# 30    930      931      1.001
# 60    3660     3661     1.000
# 120   14520    14521    1.000

Both variants coincide with the reference on all 3,000 cases. In the random graphs of 1,000, 10,000 and 100,000 vertices (with four edges per vertex) the heap insertions are about 1.27 per vertex (1,268, 12,878, 129,181): an insertion happens only when an estimate improves, and in a random graph few estimates improve more than once. The theoretical limit is E + 1, because every edge can improve an estimate at most once (when its starting vertex is popped), and the line push <= E+1 confirms it on these instances.

The worst case is constructed: k vertices u_i with distance i, each linked to each of k vertices v with weight 2k - 2i, so that every newly popped u_i improves all the v. Here the insertions reach E + 1 (931 on 930 edges for k = 30, 14,521 on 14,520 for k = 120): every edge costs an insertion. Every heap operation costs O(log of how many elements there are), and the elements are at most E + 1, so the total cost is O(E log V): the page on shortest paths in Algorithms, 4th edition states for Dijkstra with non-negative weights a time proportional to E log V in the worst case, and the count above shows its origin, one insertion per edge.

Where Dijkstra goes wrong: a negative edge

Dijkstra's argument ("a longer path can only cost more") falls apart as soon as an edge can cost less than zero. The minimal counterexample: from s you go to a at cost 1 and to b at cost 4; from b to a the edge costs -5, so the true cost of a is 4 - 5 = -1; from a you go to c at cost 1, so the true cost of c is 0. Dijkstra closes a at cost 1 before examining b, and when it discovers the better path it does not reopen a: the cost of c stays 2.

import random
from sp import dijkstra, bellman_ford, INF

# s -> a costs 1, s -> b costs 4, b -> a costs -5, a -> c costs 1:
# the true cost of a is 4 - 5 = -1, that of c is 0
archi = [(0, 1, 1), (0, 2, 4), (2, 1, -5), (1, 3, 1)]
adj = [[] for _ in range(4)]
for u, v, w in archi:
    adj[u].append((v, w))
print("Dijkstra classic      :", dijkstra(adj, 0, definitivo=True)[0])
print("Dijkstra with recheck :", dijkstra(adj, 0, definitivo=False)[0])
print("Bellman-Ford          :", bellman_ford(4, archi, 0)[0])

# how often is it wrong? acyclic graphs (hence no negative cycles), weights from -5 to 10
rng = random.Random(8)
errori = {True: 0, False: 0}
for prova in range(3000):
    n = rng.randint(2, 12)
    archi = [(u, v, rng.randint(-5, 10)) for u in range(n) for v in range(u + 1, n) if rng.random() < 0.4]
    adj = [[] for _ in range(n)]
    for u, v, w in archi:
        adj[u].append((v, w))
    ref = bellman_ford(n, archi, 0)[0]
    for modo in (True, False):
        if dijkstra(adj, 0, definitivo=modo)[0] != ref:
            errori[modo] += 1
print("on 3000 acyclic graphs with negative weights, wrong answers:")
print("  classic:", errori[True], " with recheck:", errori[False])

# Output:
# Dijkstra classic      : [0, -1, 4, 2]
# Dijkstra with recheck : [0, -1, 4, 0]
# Bellman-Ford          : [0, -1, 4, 0]
# on 3000 acyclic graphs with negative weights, wrong answers:
#   classic: 105  with recheck: 0

The classic version gives c = 2 against the true value 0; the variant with recheck and Bellman-Ford give 0. The second number says how often it happens: on 3,000 random acyclic graphs (acyclic, so that no negative cycle makes the problem unsolvable) with weights from -5 to 10, the classic version is wrong 105 times, the one with recheck never. This does not make the variant with recheck a solution: the Sedgewick and Wayne page ties the E log V guarantee to the rule that every vertex enters only once, and without that rule I have not measured a bound on its worst case. The honest choice with negative weights is another algorithm.

Bellman-Ford: relax all the edges, V - 1 times

Bellman-Ford does not choose a clever order: it relaxes all the edges, in any order, and repeats. After the first pass the estimates of all vertices reachable with a single edge are right; after the second, those reachable with two; a minimum path without cycles has at most V - 1 edges, so after V - 1 passes everything is right. If at the V-th pass an estimate still changes, there is a reachable negative-cost cycle, and the distances make no sense because going round the cycle lowers them forever. The function in sp.py stops early if a pass changes nothing, and returns the passes and checks made.

The listing tests two things. On 3,000 random digraphs with weights from -4 to 9 it compares the recognition of the negative cycle with networkx.negative_edge_cycle (starting from 0 on every vertex, to find cycles anywhere) and, where there is no cycle, the distances with single_source_bellman_ford_path_length. Then it measures the cost on a chain of V vertices, with the edges listed in the worst order (from last to first, so that each pass fixes only one vertex) and in random order.

import random
import networkx as nx
from sp import bellman_ford, INF

rng = random.Random(2)
con_ciclo = 0
for prova in range(3000):
    n = rng.randint(2, 14)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 3)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v, weight=rng.randint(-4, 9))
    archi = [(u, v, G[u][v]["weight"]) for u, v in G.edges]
    _, ciclo, _, _ = bellman_ford(n, archi)  # start from every vertex
    assert ciclo == nx.negative_edge_cycle(G)
    if ciclo:
        con_ciclo += 1
        continue
    s = rng.randrange(n)
    dist, _, _, _ = bellman_ford(n, archi, s)
    ref = nx.single_source_bellman_ford_path_length(G, s)
    assert {v: d for v, d in enumerate(dist) if d < INF} == ref
print("3000 digraphs with weights -4..9:", con_ciclo, "with a negative cycle, recognized as by networkx", nx.__version__, "; in the others the distances coincide")

print("chain of V vertices, edges listed in the worst order or in random order")
print("V     E      passes(worst)      checks       V*E       passes(random)")
for n in (50, 100, 200):
    catena = [(i, i + 1, 1) for i in range(n - 1)]
    _, _, p1, c1 = bellman_ford(n, catena[::-1], 0)
    mescolata = catena[:]
    rng.shuffle(mescolata)
    _, _, p2, _ = bellman_ford(n, mescolata, 0)
    print(f"{n:<5} {n - 1:<6} {p1:<18} {c1:<12} {n * (n - 1):<9} {p2}")

# Output:
# 3000 digraphs with weights -4..9: 682 with a negative cycle, recognized as by networkx 3.7 ; in the others the distances coincide
# chain of V vertices, edges listed in the worst order or in random order
# V     E      passes(worst)      checks       V*E       passes(random)
# 50    49     50                 2450         2450      28
# 100   99     100                9900         9900      51
# 200   199    200                39800        39800     101

682 of the 3,000 graphs have a negative cycle and all are recognized as networkx does; in the others the distances coincide. In the chain with the worst order the passes are V (the V - 1 needed plus the one that changes nothing) and the checks are exactly V times E: 2,450 for V = 50, 9,900 for V = 100, 39,800 for V = 200. This is the case that gives O(VE). With the same edges in random order the passes drop to about half of V (28, 51, 101), because a pass propagates far when the order is favourable; a linear dependence on V remains, and no guarantee better than the worst case.

Floyd-Warshall: all pairs with three loops

If you need the distances between all pairs of vertices, repeating an algorithm from every source works, but Floyd-Warshall (Floyd, Algorithm 97: Shortest path, 1962, building on Warshall's theorem on Boolean matrices, 1962) uses a more compact idea: it considers the vertices one at a time as possible intermediate stops. After treating vertices 0 to k, d[i][j] is the minimum cost from i to j using as stops only vertices between 0 and k. To go from k - 1 to k one question per pair is enough: is it worth going through k? That is d[i][k] + d[k][j] < d[i][j]. The order of the loops matters: k must be the outermost, otherwise stop k is not yet "available" when it is needed.

The starting matrix has 0 on the diagonal, the weight where there is an edge and infinity elsewhere. If at the end an element of the diagonal is negative, a negative cycle exists: the vertex returns to itself at a cost below zero. The listing verifies both things against networkx (matrices equal to floyd_warshall, cycle recognized as negative_edge_cycle) on 1,500 digraphs with weights from -3 to 9, and counts the steps of the innermost loop on graphs of 20, 40 and 80 vertices.

import random
import networkx as nx
from sp import floyd_warshall, INF

def matrice(n, G):
    m = [[INF] * n for _ in range(n)]
    for i in range(n):
        m[i][i] = 0
    for u, v, w in G.edges(data="weight"):
        m[u][v] = min(m[u][v], w)
    return m

rng = random.Random(4)
ciclo_neg = 0
for prova in range(1500):
    n = rng.randint(1, 12)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 3)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v, weight=rng.randint(-3, 9))
    d, _ = floyd_warshall(n, matrice(n, G))
    negativo = any(d[i][i] < 0 for i in range(n))  # a vertex that returns to itself at cost < 0
    assert negativo == nx.negative_edge_cycle(G)
    if negativo:
        ciclo_neg += 1
        continue
    ref = dict(nx.floyd_warshall(G))
    assert all(d[i][j] == ref[i][j] for i in range(n) for j in range(n))
print("1500 digraphs: matrices equal to networkx", nx.__version__, "and", ciclo_neg, "negative cycles recognized from the diagonal")

print("V     inner-loop steps          V^3")
for n in (20, 40, 80):
    G = nx.gnp_random_graph(n, 0.2, seed=3, directed=True)
    nx.set_edge_attributes(G, 1, "weight")
    _, passi = floyd_warshall(n, matrice(n, G))
    print(f"{n:<5} {passi:<25} {n ** 3}")

# Output:
# 1500 digraphs: matrices equal to networkx 3.7 and 207 negative cycles recognized from the diagonal
# V     inner-loop steps          V^3
# 20    8000                      8000
# 40    64000                     64000
# 80    512000                    512000

The steps are exactly V cubed: 8,000, 64,000 and 512,000. The cost does not depend on how many edges there are, because the triple loop scans every pair for every stop: on a graph with few edges Floyd-Warshall costs the same as on a dense one, which is why it is worthwhile only when V is small or the graph is dense. In memory it occupies V squared cells, the same adjacency matrix as in article 6.

Which one to choose

The three measured complexities fit in a table; the costs are those of the counts above, not times.

The correct reading of the table starts from the question about the weights: if they are all non-negative (distances, times, real costs) Dijkstra is the choice, because it is the only one of the three whose cost depends on E log V rather than on V times E. If they can be negative (a gain in a chain of exchanges, a discount), you need Bellman-Ford for one source, and its negative-cycle check tells you whether the question makes sense. If you need all pairs and V is small, Floyd-Warshall is the shortest to write and to check.

I have neither implemented nor measured A*, contraction hierarchies or bidirectional search: I name them because faster variants of Dijkstra exist for very large road graphs, but any claim about how much faster they are would require a source and a measurement that are not here.

What you can do now

You can relax an edge, you can write Dijkstra with heapq without the "lower the priority" operation, and you can build the four-vertex graph on which it fails because a weight is negative. You can use Bellman-Ford to recognize a negative cycle, and Floyd-Warshall to get all the distances with a triple loop. You can state each of the three complexities with a count that the reader can redo. The next and last article, on dynamic programming, takes up the idea of Floyd-Warshall (a table filled in an order that guarantees the partial answers are already there) and applies it to problems that are not graphs. To go back, article 5 is the one on trees.

The same series also includes complexity and Big-O, sorting algorithms, linear search, binary search and hash tables and lists, stacks, queues and deques.