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.
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 with integer edge weights, some negative, and for every pair the least weight of a path from to .
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 costs , 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 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 be the cheapest walk from the source to with at most edges. Then
and with no negative cycle a cheapest path has at most edges, so .
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 relaxations: 200 here, and on a dense graph. The waste is easy to point at. The cheapest way from 1 to 3 (through 2, for ) 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 . The intermediate vertices of a path are all of its vertices except the two ends. Then:
- State: is the least weight of an path whose intermediate vertices all lie in .
- Base case: allows no intermediates at all, so it is the road table: 0 on the diagonal, for a road, otherwise.
- Answer location: , where every vertex is allowed.
- Recurrence: .
Invariant
After round (the round that allows vertex ), is the least weight of an path whose intermediate vertices all lie in .
Proof of the recurrence. Take a cheapest path whose intermediates lie in . With no negative cycle, cutting a cycle out never makes a path dearer, so can be taken simple: it visits at most once. If avoids , it is one of the paths ranges over. If passes through , cut it there. The part and the part are paths whose intermediates lie in , so they cost at least and . So 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 , the invariant holds after every round, and after round it says .
The step not to skip is "visits 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.
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 costs , better than no road; costs 7, worse than the direct road's 6.
Follow cell . It is until round 2 allows (cost 9). Round 3 lowers it to 6, the path , by adding two entries the table already had: and . That is the saving over Bellman–Ford: the cheapest 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 , 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 tables; the algorithm keeps one.
Evaluation order: outermost. Round reads only values, so the whole of round must finish before round starts. Inside a round, and may go in any order.
Space optimization: overwrite in place. In round , a cell in row is offered , because . So row and column do not change during round , and they are the only other cells the round reads. Overwriting early is therefore harmless, and numbers do the work of .
Predict: the same three loops with 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, stays , although costs , and comes out 6 instead of 3.
Reconstruction. The dispatcher also wants routes. Keep a second matrix, nxt[i][j]: the
vertex after on the best path so far ( itself for a road). When round
improves , the new path starts the way the best path starts, so nxt[i][j]
becomes nxt[i][k]. To read a path, follow nxt[·][j] from until you reach . On H this
gives for the pair
and for :
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 | : cheapest path with intermediates in |
| Recurrence | avoid , or pass through it once: |
| Base case | the road table, 0 on the diagonal |
| Evaluation order | outermost; in any order |
| Answer location | , every pair |
| Reconstruction | nxt[i][j] = nxt[i][k] on every improvement |
| Space optimization | one 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 (cost 1), (cost ) and (cost 0): the cycle costs .
Theorem. After Floyd–Warshall, some diagonal entry 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 is a negative closed walk, which contains a negative cycle. Conversely, take a negative cycle, its highest-numbered vertex , and another vertex on it. The arc of the cycle from to and the arc from back to pass only through vertices below , so just before round the table already holds values no larger than their weights, and round sets 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 , 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:
matrix multiplication with in place of . If holds the cheapest walks of at most edges, then : split a walk after its -th edge. Squaring times reaches edges, and each product costs steps.
assert math.ceil(math.log2(5 - 1)) * 5**3 == 250
| Method on H | steps | growth |
|---|---|---|
| Bellman–Ford from every source | 200 relaxations | , up to |
| min-plus squaring | 250 inner steps | |
| Floyd–Warshall | 125 relaxations | exactly |
The three programs compute the same numbers and differ only in their state. Doubling the edge count takes full products. Adding one allowed vertex costs , and of them finish the job. On H that is a factor of 2; at it is .
The cost of all pairs
Floyd–Warshall makes exactly relaxations: rounds of cells. That is
time and space in the worst case, on a word-RAM where adding two
weights costs . Reading a path out of nxt costs its length.
Is optimal? Nobody knows. For general weights no algorithm is known, for any fixed , 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 , which wins when is far below . 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 orders of the rest. Going to the nearest unvisited city each time is fast, but here it builds 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 with those that start . Both prefixes have visited 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: is the cheapest path that starts at city 0, visits exactly the cities in in some order, and ends at .
- Recurrence: the city before is some , and the path up to is itself a path from 0 through exactly that ends at . So .
- Base case: .
- Evaluation order: increasing as a bitmask: every proper subset is a smaller number.
- Answer location: .
- Reconstruction: store each minimizing ; walk back from the best last city, removing it from at each step.
- Space optimization: two layers of values suffice, but then the stored 's, and the tour, are gone.
Invariant
Once the sets of size are filled, every with is the cheapest way to leave city 0, visit exactly , and stand at . The rest of the tour depends only on , never on the order inside .
Here is part of the table for the six cities: the sets inside .
Column holds 28 (via 2, ending at 1) and 21 (via 1, ending at 2). The cell arrives at 3 from one of them: or . 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 , the cities visited so far?
The next leg's cost depends on where you stand. Two paths through the same ending at different cities can't be compared: one may be cheaper so far and dearer to continue. The state needs : entries, not .
Complexity: against
One transition reads and adds . A set of size has choices of and choices of , so the number of transitions is
The table has entries. So Held–Karp takes time and 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 | Held–Karp transitions | tours |
|---|---|---|
| 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 . 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 , 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 , is optimal over paths whose intermediates lie in ; is optimal over paths from 0 through exactly to .
- Complexity achieved: exactly relaxations for all pairs, against by squaring and by Bellman–Ford from every source; transitions for TSP, against tours. All worst case, on a word-RAM.
- Failure mode: the loop not outermost: the code runs, and some distances are wrong.
- In real software: SciPy's
scipy.sparse.csgraph.floyd_warshalland NetworkX'sfloyd_warshallimplement this algorithm; NetworkX'sfloyd_warshall_predecessor_and_distancealso returns the predecessor matrix used to read paths. - Next: module 11 turns to a different design tool, network flow: how much can get from to , and the cut that proves no more can.
Check yourself
After the lab, the tutor will ask you to defend your work out loud:
- State Floyd–Warshall's invariant, and prove that one round preserves it. Where exactly does "no negative cycle" enter the proof?
- Add a road from depot 2 to depot 0 with cost to H. What does the algorithm output, what should it report, and which pairs are affected?
- 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.