Dynamic Programming II

The question

Five depots are joined by ten one-way roads. Each road has a cost, and two of them pay you to use them: a rebate of 3 on the road from depot 1 to depot 2, and a rebate of 2 on the road from depot 4 back to depot 0.

4 9 -3 7 2 3 -2 8 2 6 0 1 2 3 4
The depot network H used throughout this lesson: n = 5 depots, 10 one-way roads. The two roads with negative cost are drawn in colour.

The dispatcher wants the cheapest cost from every depot to every other: 20 answers, one for each ordered pair. This is the all-pairs shortest paths problem (APSP): a directed graph on vertices 0,…,n−1 with integer edge weights, some negative, and for every pair the least weight δ(i,j) of a path from i to j.

Negative roads are fine. A negative cycle is not: if a loop paid you to drive round it, "cheapest" would mean nothing. H has no such loop: the cycle 1→2→3→1 costs −3+2+2=1, and no cycle costs less.

This module assumes module 09, which describes every dynamic program by seven decisions, and Dijkstra's and Bellman–Ford's single-source algorithms from a first algorithms course.

import itertools, math

INF = float("inf")
H = [(0, 1, 4), (0, 2, 9), (1, 2, -3), (1, 3, 7), (2, 3, 2), (3, 4, 3), (4, 0, -2),
     (2, 4, 8), (3, 1, 2), (4, 2, 6)]
W = {(u, v): w for u, v, w in H}

def weight(path):
    return sum(W[a, b] for a, b in zip(path, path[1:]))

def simple_paths(n, s, t):
    """Every path from s to t that repeats no vertex."""
    for size in range(n - 1):
        for middle in itertools.permutations([v for v in range(n) if v not in (s, t)], size):
            path = (s,) + middle + (t,)
            if all(edge in W for edge in zip(path, path[1:])):
                yield path

# No negative cycle: the cheapest simple cycle is 1 -> 2 -> 3 -> 1, of cost 1.
cycles = [weight(p) + W[p[-1], s] for s in range(5) for t in range(5)
          if t != s and (t, s) in W for p in simple_paths(5, s, t)]
assert min(cycles) == 1 == weight((1, 2, 3, 1))

DELTA = [[0 if s == t else min((weight(p) for p in simple_paths(5, s, t)), default=INF)
          for t in range(5)] for s in range(5)]
assert DELTA == [[0, 4, 1, 3, 6], [0, 0, -3, -1, 2], [3, 4, 0, 2, 5],
                 [1, 2, -1, 0, 3], [-2, 2, -1, 1, 0]]

From depot 0, depot 2 costs 1 (through depot 1 and its rebate), and from depot 1 you can get back to depot 0 for nothing. The table was found by trying every simple path: fine for five depots, hopeless for fifty.

The naive approaches, and where they waste work

Try every path. A complete graph has more than (n−2)! simple paths between two vertices. That is exponential.

Run a single-source algorithm from every source. Dijkstra's algorithm is out: it settles a vertex for good once it is the closest unsettled one, and a negative road found later can undercut it. Bellman–Ford handles negative weights. It is itself a dynamic program: let dm(v) be the cheapest walk from the source to v with at most m edges. Then

dm(v)=min(dm−1(v), min(u,v)dm−1(u)+w(u,v)),

and with no negative cycle a cheapest path has at most n−1 edges, so dn−1=δ.

def bellman_ford(n, edges, s, count):
    d = [INF] * n
    d[s] = 0
    for m in range(1, n):                     # d holds d_{m-1}; build d_m
        nd = list(d)
        for u, v, w in edges:
            count[0] += 1                     # one relaxation
            if d[u] + w < nd[v]:
                nd[v] = d[u] + w
        d = nd
    return d

count = [0]
assert [bellman_ford(5, H, s, count) for s in range(5)] == DELTA
assert count[0] == 5 * 4 * 10 == 200

From every source that is n·(n−1)·E relaxations: 200 here, and Θ(n4) on a dense graph. The waste is easy to point at. The cheapest way from 1 to 3 (through 2, for −1) is rediscovered inside the run from 0, the run from 4, and the run from 1 itself. The five runs never share what they learn.

Choosing the state

A dynamic program needs a parameter that grows. Bellman–Ford's is the number of edges a path may use. Floyd and Warshall's idea (1962) is which vertices a path may pass through.

Number the vertices 0,…,n−1. The intermediate vertices of a path are all of its vertices except the two ends. Then:

  • State: Dk[i][j] is the least weight of an i→j path whose intermediate vertices all lie in {0,…,k−1}.
  • Base case: D0 allows no intermediates at all, so it is the road table: 0 on the diagonal, w(i,j) for a road, ∞ otherwise.
  • Answer location: Dn[i][j], where every vertex is allowed.
  • Recurrence: Dk+1[i][j]=min(Dk[i][j], Dk[i][k]+Dk[k][j]).

Invariant

After round k (the round that allows vertex k), D[i][j] is the least weight of an i→j path whose intermediate vertices all lie in {0,…,k}.

Proof of the recurrence. Take a cheapest i→j path P whose intermediates lie in {0,…,k}. With no negative cycle, cutting a cycle out never makes a path dearer, so P can be taken simple: it visits k at most once. If P avoids k, it is one of the paths Dk[i][j] ranges over. If P passes through k, cut it there. The part i→k and the part k→j are paths whose intermediates lie in {0,…,k−1}, so they cost at least Dk[i][k] and Dk[k][j]. So Dk+1[i][j] is at least the minimum on the right. It is also at most that minimum, because each term is the weight of a real walk that uses only allowed vertices, and a walk can be shortened to a path without getting dearer. By induction on k, the invariant holds after every round, and after round n−1 it says D=δ.

The step not to skip is "visits k at most once": that is where "no negative cycle" enters. With a negative cycle the loops still run, but they no longer compute path weights.

Floyd–Warshall on the depots

The road table, then one table per round; a shaded cell improved in that round.

Floyd–Warshall on H1D(0)2Via 03Via 14Via 25Via 36Via 4
D(0) is the road table; the panel "Via k" is the table after round k, when paths may also pass through vertex k. Round by round, 1, 5, 5, 4 and 3 cells improve: 18 improvements in all.
Predict: why does round 0 improve only one cell, (4, 1)?

A path that newly uses vertex 0 must enter and leave it. Only one road enters 0, from 4, so only row 4 can change. The detour 4→0→1 costs −2+4=2, better than no road; 4→0→2 costs 7, worse than the direct road's 6.

Follow cell (0,4). It is ∞ until round 2 allows 0→1→2→4 (cost 9). Round 3 lowers it to 6, the path 0→1→2→3→4, by adding two entries the table already had: D[0][3]=3 and D[3][4]=3. That is the saving over Bellman–Ford: the cheapest 0→3 trip is computed once and reused by every pair that passes through 3.

The invariant can be tested: find by brute force the cheapest paths with intermediates in {0,…,k}, and compare consecutive tables.

def restricted(n, k):
    """Brute force: cheapest paths whose intermediate vertices are all below k + 1."""
    return [[0 if s == t else min((weight(p) for p in simple_paths(n, s, t)
                                   if all(v <= k for v in p[1:-1])), default=INF)
             for t in range(n)] for s in range(n)]

tables = [restricted(5, k) for k in range(-1, 5)]    # k = -1: no intermediates, the road table
changes = [[(i, j, before[i][j], after[i][j]) for i in range(5) for j in range(5)
            if after[i][j] != before[i][j]] for before, after in zip(tables, tables[1:])]
assert [len(c) for c in changes] == [1, 5, 5, 4, 3]
assert changes[0] == [(4, 1, INF, 2)] and (0, 2, 9, 1) in changes[1]
assert (0, 4, INF, 9) in changes[2] and (0, 4, 9, 6) in changes[3]
assert changes[4] == [(1, 0, INF, 0), (2, 0, INF, 3), (3, 0, INF, 1)]
assert tables[-1] == DELTA

One table, filled in place

The recurrence names n+1 tables; the algorithm keeps one.

Evaluation order: k outermost. Round k reads only Dk values, so the whole of round k must finish before round k+1 starts. Inside a round, i and j may go in any order.

Space optimization: overwrite in place. In round k, a cell in row k is offered D[k][k]+D[k][j]=D[k][j], because D[k][k]=0. So row k and column k do not change during round k, and they are the only other cells the round reads. Overwriting early is therefore harmless, and n2 numbers do the work of n3.

Predict: the same three loops with k innermost (for i: for j: for k:). Same table?

No, and the code runs without complaint. Each cell is finished before the cells it needs are. On H, D[1][0] stays ∞, although 1→2→3→4→0 costs −3+2+3−2=0, and D[2][0] comes out 6 instead of 3.

Reconstruction. The dispatcher also wants routes. Keep a second matrix, nxt[i][j]: the vertex after i on the best i→j path so far (j itself for a road). When round k improves (i,j), the new path starts the way the best i→k path starts, so nxt[i][j] becomes nxt[i][k]. To read a path, follow nxt[·][j] from i until you reach j. On H this gives 0→1→2→3→4 for the pair (0,4) and 1→2→3→4→0 for (1,0):

assert weight((0, 1, 2, 3, 4)) == DELTA[0][4] == 6
assert weight((1, 2, 3, 4, 0)) == DELTA[1][0] == 0

That completes the seven decisions:

Decision Floyd–Warshall
State Dk[i][j]: cheapest i→j path with intermediates in {0..k−1}
Recurrence avoid k, or pass through it once: min(Dk[i][j],Dk[i][k]+Dk[k][j])
Base case the road table, 0 on the diagonal
Evaluation order k outermost; i,j in any order
Answer location Dn[i][j], every pair
Reconstruction nxt[i][j] = nxt[i][k] on every improvement
Space optimization one n×n table, updated in place

Negative cycles show up on the diagonal

Run the same loops on a graph that has a negative cycle and they still finish. They just no longer compute path weights. Take three vertices with edges 0→1 (cost 1), 1→2 (cost −2) and 2→0 (cost 0): the cycle costs 1−2+0=−1.

Theorem. After Floyd–Warshall, some diagonal entry D[i][i] is negative if and only if the graph has a negative cycle.

Sketch. Every entry the algorithm writes is the weight of a real walk, so a negative D[i][i] is a negative closed walk, which contains a negative cycle. Conversely, take a negative cycle, its highest-numbered vertex k, and another vertex i on it. The arc of the cycle from i to k and the arc from k back to i pass only through vertices below k, so just before round k the table already holds values no larger than their weights, and round k sets D[i][i] to at most the cycle's weight, below 0. On the three-vertex example the diagonal ends negative. So the test is one line after the loops: if any D[i][i]<0, report the cycle instead of distances.

Two other programs for the same table

Bellman–Ford's state counts edges, and you can count edges for all pairs at once with the min-plus product:

(A⊗B)[i][j]=mink(A[i][k]+B[k][j]),

matrix multiplication with (min,+) in place of (+,×). If D(m) holds the cheapest walks of at most m edges, then D(2m)=D(m)⊗D(m): split a walk after its m-th edge. Squaring ⌈log2(n−1)⌉ times reaches n−1 edges, and each product costs n3 steps.

assert math.ceil(math.log2(5 - 1)) * 5**3 == 250
Method on H steps growth
Bellman–Ford from every source 200 relaxations n2E, up to n4
min-plus squaring 250 inner steps n3logn
Floyd–Warshall 125 relaxations n3 exactly

The three programs compute the same numbers and differ only in their state. Doubling the edge count takes logn full products. Adding one allowed vertex costs n2, and n of them finish the job. On H that is a factor of 2; at n=1000 it is ⌈log2999⌉=10.

The cost of all pairs

Floyd–Warshall makes exactly n3 relaxations: n rounds of n2 cells. That is Θ(n3) time and Θ(n2) space in the worst case, on a word-RAM where adding two weights costs O(1). Reading a path out of nxt costs its length.

Is n3 optimal? Nobody knows. For general weights no O(n3−ε) algorithm is known, for any fixed ε>0, and a well-studied conjecture says none exists. Sparse graphs are different: with non-negative weights, Dijkstra from every source with a binary heap costs O(nElogn), which wins when E is far below n2. Module 13's vertex potentials remove negative weights without changing which paths are shortest, and open that route too.

Tours through every city

Now a different question, with the same kind of answer. A courier leaves city 0, must visit cities 1 to 5 once each, and must come back. The costs of the legs are symmetric:

0 1 2 3 4 5
0 0 12 19 8 23 7
1 12 0 9 14 17 21
2 19 9 0 11 7 16
3 8 14 11 0 13 6
4 23 17 7 13 0 10
5 7 21 16 6 10 0

This is the travelling salesperson problem (TSP). Brute force fixes city 0 first and tries the (n−1)!=120 orders of the rest. Going to the nearest unvisited city each time is fast, but here it builds 0,5,3,2,4,1 for a cost of 60. The best tour costs 52.

T = [[0, 12, 19, 8, 23, 7], [12, 0, 9, 14, 17, 21], [19, 9, 0, 11, 7, 16],
     [8, 14, 11, 0, 13, 6], [23, 17, 7, 13, 0, 10], [7, 21, 16, 6, 10, 0]]

def tour_cost(tour):
    return sum(T[a][b] for a, b in zip(tour, tour[1:] + tour[:1]))

tours = [[0] + list(p) for p in itertools.permutations(range(1, 6))]
assert len(tours) == 120
assert min(map(tour_cost, tours)) == 52 == tour_cost([0, 1, 2, 4, 5, 3])

here, greedy = 0, [0]
while len(greedy) < 6:
    here = min((c for c in range(6) if c not in greedy), key=lambda c: T[here][c])
    greedy.append(here)
assert greedy == [0, 5, 3, 2, 4, 1] and tour_cost(greedy) == 60

Where does brute force waste work? Compare the tours that start 0,1,2,3 with those that start 0,2,1,3. Both prefixes have visited {1,2,3} and stand at city 3, so every completion costs the same after either. Only the cheaper prefix can begin a best tour, yet brute force prices both families in full.

Held–Karp: the set, and where you stand

That observation is the whole algorithm (Bellman, and independently Held and Karp, 1962):

  • State: dp[S][j] is the cheapest path that starts at city 0, visits exactly the cities in S⊆{1,…,n−1} in some order, and ends at j∈S.
  • Recurrence: the city before j is some k∈S⧵{j}, and the path up to k is itself a path from 0 through exactly S⧵{j} that ends at k. So dp[S][j]=mink∈S⧵{j}(dp[S⧵{j}][k]+d(k,j)).
  • Base case: dp[{j}][j]=d(0,j).
  • Evaluation order: increasing S as a bitmask: every proper subset is a smaller number.
  • Answer location: minj(dp[{1..n−1}][j]+d(j,0)).
  • Reconstruction: store each minimizing k; walk back from the best last city, removing it from S at each step.
  • Space optimization: two layers of values suffice, but then the stored k's, and the tour, are gone.

Invariant

Once the sets of size s are filled, every dp[S][j] with |S|≤s is the cheapest way to leave city 0, visit exactly S, and stand at j. The rest of the tour depends only on (S,j), never on the order inside S.

Here is part of the table for the six cities: the sets inside {1,2,3}.

1 2 3 1,2 1,3 2,3 1,2,3 j=1 j=2 j=3 12 – – 28 22 – 28 – 19 – 21 – 19 31 – – 8 – 26 30 32
dp[S][j] for the sets S inside {1, 2, 3}: each column is a set S, each row the end city j, and a dash means j is not in S. The shaded cell dp[{1,2,3}][3] = 32 comes from the column {1, 2}: the cheaper of 28 + d(1,3) = 42 and 21 + d(2,3) = 32.

Column {1,2} holds 28 (via 2, ending at 1) and 21 (via 1, ending at 2). The cell dp[{1,2,3}][3] arrives at 3 from one of them: 28+14=42 or 21+11=32. It keeps 32 and remembers "came from 2". Checked by brute force:

def cheapest_path(S, j):
    """Brute force: the cheapest path 0 -> (all of S, in some order) ending at j."""
    return min(sum(T[a][b] for a, b in zip((0,) + p + (j,), p + (j,)))
               for p in itertools.permutations(sorted(set(S) - {j})))

assert [cheapest_path({1, 2}, j) for j in (1, 2)] == [28, 21]
assert cheapest_path({1, 2, 3}, 3) == min(28 + T[1][3], 21 + T[2][3]) == 32
Predict: why not make the state just the set S, the cities visited so far?

The next leg's cost d(k,j) depends on where you stand. Two paths through the same S ending at different cities can't be compared: one may be cheaper so far and dearer to continue. The state needs (S,j): 2n−1(n−1) entries, not 2n−1.

Complexity: n22n against n!

One transition reads dp[S⧵{j}][k] and adds d(k,j). A set of size s has s choices of j and s−1 choices of k, so the number of transitions is

∑s=2n−1(n−1s)s(s−1)=(n−1)(n−2)2n−3≤n22n8.

The table has (n−1)2n−1 entries. So Held–Karp takes O(n22n) time and O(n2n) space in the worst case, on a word-RAM. For six cities that is 160 transitions against 120 tours, so brute force still wins. Not for long:

def transitions(n):
    return sum(math.comb(n - 1, s) * s * (s - 1) for s in range(2, n))

assert transitions(6) == 160
assert all(transitions(n) == (n - 1) * (n - 2) * 2 ** (n - 3) <= n * n * 2 ** n / 8
           for n in range(3, 25))
assert [transitions(n) for n in (8, 10, 12, 14)] == [1344, 9216, 56320, 319488]
assert math.factorial(13) == 6227020800
cities n Held–Karp transitions tours (n−1)!
8 1,344 5,040
10 9,216 362,880
12 56,320 39,916,800
14 319,488 6,227,020,800

It is still exponential: TSP is NP-hard, so no polynomial algorithm is expected, and the table alone has about 400 million entries at n=25. But within a billion steps, brute force handles 13 cities and Held–Karp handles 23.

assert 24 * 2**24 == 402_653_184                             # (n - 1) 2^(n - 1) entries, n = 25
assert math.factorial(12) < 10**9 < math.factorial(13)       # tours for 13 and 14 cities
assert transitions(23) < 10**9 < transitions(24)

A problem that looks different

A trader swaps between six currencies at posted rates: 1 unit of A buys 0.9 of B, 1 of B buys 1.2 of C, and so on for every pair. Can some sequence of trades start and end with A and make money? Nothing here mentions roads, yet one table from this lesson answers it once each rate is replaced by a different number. Which table, and which number? The lab's last problem is a different one.

Practise

In the lab you write Floyd–Warshall on a new graph and watch each round's table, every frame checked against a brute-force search. You predict a Held–Karp table layer by layer on five new cities, then write Held–Karp over bitmasks under an operation budget that trying every order cannot meet. You measure min-plus squaring against n3, and finish with one problem that doesn't say what it is.

Recap

  • You can now: compute all-pairs cheapest paths with negative edges and read the paths; detect a negative cycle from the diagonal; compare three dynamic programs for the same table by their state; and solve the travelling salesperson problem exactly for about 20 cities.
  • Invariant: after round k, D[i][j] is optimal over paths whose intermediates lie in {0,…,k}; dp[S][j] is optimal over paths from 0 through exactly S to j.
  • Complexity achieved: exactly n3 relaxations for all pairs, against n3logn by squaring and n2E by Bellman–Ford from every source; (n−1)(n−2)2n−3 transitions for TSP, against (n−1)! tours. All worst case, on a word-RAM.
  • Failure mode: the k loop not outermost: the code runs, and some distances are wrong.
  • In real software: SciPy's scipy.sparse.csgraph.floyd_warshall and NetworkX's floyd_warshall implement this algorithm; NetworkX's floyd_warshall_predecessor_and_distance also returns the predecessor matrix used to read paths.
  • Next: module 11 turns to a different design tool, network flow: how much can get from s to t, and the cut that proves no more can.

Check yourself

After the lab, the tutor will ask you to defend your work out loud:

  1. State Floyd–Warshall's invariant, and prove that one round preserves it. Where exactly does "no negative cycle" enter the proof?
  2. Add a road from depot 2 to depot 0 with cost −4 to H. What does the algorithm output, what should it report, and which pairs are affected?
  3. For a dense graph on 1,000 vertices with some negative weights, compare Floyd–Warshall, min-plus squaring and Bellman–Ford from every source by their step counts.

Topic list modeled on the public schedule of CMU 15-451/651 Algorithm Design and Analysis (Fall 2025). All lecture text, proofs, examples, code, and exercises are original to SciMigo; no CMU course materials are reproduced. This course is independent and not affiliated with or endorsed by Carnegie Mellon University.

Your turn

Now try it yourself

Write the code, press Run, and scrub through the frames your program draws. The checks tell you when it works.

Loading the lab…