Dynamic programming: Held-Karp

Read first

• Post 2: Why is the TSP hard? Brute force checks all (n − 1)!/2 tours. That is 2,520 for eight stops and about 7.76 × 10⁶⁵ for berlin52. The TSP is NP-hard, so no known exact method is polynomial.
• Primer B: algorithms and complexity Big-O notation, recursion, and the difference between n! and 2ⁿ. Both are exponential, but n! grows much faster.

The question this post answers: can we find the optimal tour, with proof, without checking every tour?

Intuition

Brute force wastes most of its time repeating itself. Take two partial tours on our stops that start at A, visit B, C and D in some order, and end at D:

• A → B → C → D costs 41 + 36 + 41 = 118
• A → C → B → D costs 71 + 36 + 63 = 170

Whatever happens after D is the same for both: the rest of the tour has to visit the same remaining stops and return to A. So the second partial tour can never be part of an optimal tour. Brute force does not know this, and it extends both, along with every other order of B, C and D.

The fix is to remember only what the future depends on: which stops have been visited, and where we are now. The order inside the visited set does not matter. That is the whole idea of the Held-Karp algorithm, published independently by Bellman (1962) and by Held and Karp (1962).

I will work it on six stops, A, B, C, D, G and H, taken from the eight-stop example in earlier posts. A is the start. Coordinates are in kilometres; distances are in units of 100 m, rounded, as before.

Step through the diagram below to compare the two partial tours from the start of this section and see the full tour that Held-Karp eventually returns. Its lengths are unrounded kilometres, a few tenths off the matrix values.

Write C(S, k) for the cheapest way to leave A, visit every stop in the set S exactly once, and finish at stop k, where k is in S. We fill in these values layer by layer, from small sets to large ones.

Layer 1: one stop. The only way to visit {k} and end at k is to go there directly, so C({k}, k) = d(A, k).

Layer 2: two stops. To visit {j, k} and end at k, you must have visited j just before, so C({j, k}, k) = C({j}, j) + d(j, k). For example, C({B, C}, C) = 41 + 36 = 77 and C({B, C}, B) = 71 + 36 = 107. All twenty entries:

Layer 3 onward: the real choice. Now there is more than one way to arrive. To finish at k, the stop before k was some m in S other than k, and the cost up to m is already in the previous layer. Take the best m:

C({B,C,D}, D) = min( C({B,C}, C) + d(C,D),   C({B,C}, B) + d(B,D) )
              = min( 77 + 41,                107 + 63 )
              = min( 118, 170 ) = 118        (best predecessor: C)

This is the comparison from the start of the section, done once and stored. Every later partial tour that needs “B, C, D visited, standing at D” reuses 118 and never looks at the 170 again.

Layer 3 has 30 entries, layer 4 has 20, and layer 5, the full set {B, C, D, G, H}, has 5. The diagram below shows the shape of the computation: each layer reads only the layer before it.

flowchart TD
  L1["Layer 1: 5 entries"] --> L2["Layer 2: 20 entries"]
  L2 --> L3["Layer 3: 30 entries"]
  L3 --> L4["Layer 4: 20 entries"]
  L4 --> L5["Layer 5: 5 entries"]
  L5 --> END["Add leg back to A"]
  END --> OPT["Minimum: 247"]

Closing the tour. Each full-set entry is a path that ends somewhere; add the leg back to A and take the minimum:

The tie between B and G is not a coincidence: they are the same tour driven in opposite directions.

Recovering the tour. Each entry also stores its best predecessor: the stop it came from. Start at the winning entry and follow those pointers back to A:

Reading the “Came from” column bottom to top gives the tour forwards:

That is the optimum for these six stops, and it comes with a proof, because every possible tour was accounted for by some chain of minima. We filled 80 table entries. Brute force on six stops checks 60 tours, so at this size there is no saving. The saving appears as n grows, and it is enormous.

The formal idea

Number the cities 1 to n, with city 1 as the start. Let d(i, j) be the cost from i to j. The algorithm does not need d(i, j) = d(j, i), so it solves the asymmetric TSP as well.

For every subset S ⊆ {2, …, n} and every k ∈ S, define C(S, k) as the minimum cost of a path that starts at city 1, visits every city in S exactly once, and ends at k.

Base:        C({k}, k) = d(1, k)                                  for k = 2, …, n

Recurrence:  C(S, k)   = min   [ C(S ∖ {k}, m) + d(m, k) ]         for |S| ≥ 2, k ∈ S
                        m ∈ S∖{k}

Answer:      OPT       = min   [ C({2, …, n}, k) + d(k, 1) ]
                        k ∈ {2,…,n}

The diagram below shows one application of the recurrence, for the entry C({B, C, D}, D) from the hand example: each state in the layer below offers one candidate, and the cheaper one wins.

flowchart TD
  S1["C({B,C}, C) = 77"] -->|"+ d(C,D) = 41"| M{"min(118, 170)"}
  S2["C({B,C}, B) = 107"] -->|"+ d(B,D) = 63"| M
  M --> R["C({B,C,D}, D) = 118"]

S ∖ {k} means S with k removed. The recurrence is correct because the last step of any such path goes from some m to k, and the part before it must itself be the cheapest path that covers S ∖ {k} and ends at m. If it were not, swapping in the cheaper one would improve the whole path. This “optimal solutions are built from optimal pieces” property is what makes dynamic programming work.

Counting the table. There are 2ⁿ⁻¹ subsets of the n − 1 non-start cities, and each subset S has |S| possible end cities. Summing gives (n − 1) · 2ⁿ⁻² entries. For n = 6 that is 5 × 16 = 80, matching the hand count.

Time. Each entry C(S, k) scans |S| − 1 predecessors. Summed over all entries, that is (n − 1)(n − 2) · 2ⁿ⁻³ additions and comparisons, which is O(n² · 2ⁿ).

Memory. To compute layer s you need layer s − 1, and to recover the tour you need the predecessor of every entry. Keeping them all costs O(n · 2ⁿ) memory.

The table below compares brute force with Held-Karp, assuming 4 bytes per table entry (a 32-bit cost, before storing predecessors).

At twenty cities, Held-Karp does about 45 million operations where brute force would price 61 quadrillion tours. Bellman’s paper, first issued as a RAND report in 1961, says the method could solve problems of up to seventeen cities on the computers of the day.

Two things stand out in the last column. First, Held-Karp is still exponential: every extra city slightly more than doubles both the table and the work. Second, memory hits the wall before time does. At 30 cities, 1.1 × 10¹¹ operations is about two minutes at an assumed billion simple operations per second, but the 31 GB table is more RAM than most laptops have.

Has anyone beaten it? For sixty years, O(n² · 2ⁿ) was the best known worst-case time for the general TSP. In 2024 Mihail Stoian posted an algorithm that runs in 2ⁿ · n² / 2^Ω(√log n) time by recasting each layer as a min-plus matrix product. That trims the polynomial factor, not the 2ⁿ. Whether the general TSP can be solved in O(1.9999ⁿ) time is still open.

Common misconceptions

“Dynamic programming makes the TSP polynomial.” It makes it less exponential. O(n² · 2ⁿ) beats O(n!) by an astronomical factor, but it still doubles with every city. NP-hardness has not gone anywhere.

“Held-Karp is brute force with caching.” Caching complete permutations would save nothing, because no two permutations are the same. The saving comes from choosing a state that forgets the order: (set visited, current stop). Many orders collapse onto one state, and only the best one survives.

“Running time is what stops it.” Memory runs out first. At 30 cities the operation count is manageable but the table needs about 31 GB, and storing predecessors to recover the tour adds more.

“It only works for the plain symmetric TSP.” The recurrence never assumes d(i, j) = d(j, i), so it handles the asymmetric TSP unchanged. Variants with time windows can be solved by dynamic programming too, by adding arrival time to the state and discarding states that violate a window. Dumas, Desrosiers, Gelinas and Solomon (1995) built an exact TSPTW algorithm that way, using the windows to cut down the state space.

Try it in Python

Here is Held-Karp in plain Python, including the walk back through the stored predecessors. It stores each subset as a bitmask, an integer whose bit i is set when city i is in the set:

import math
from itertools import combinations

STOPS = {"A": (2, 2), "B": (6, 1), "C": (9, 3), "D": (8, 7),
         "E": (5, 9), "F": (1, 8), "G": (4, 5), "H": (10, 9)}  # km

def dist_100m(p, q):  # straight-line distance in units of 100 m, rounded
    return int(10 * math.hypot(p[0] - q[0], p[1] - q[1]) + 0.5)

def held_karp(names, coords):
    n = len(names)
    d = [[dist_100m(coords[a], coords[b]) for b in names] for a in names]
    # C[(mask, k)] = cheapest path from city 0 through the cities in mask, ending at k.
    # mask is a bitmask over cities 1..n-1; city 0 is the fixed start.
    C, parent = {}, {}
    for k in range(1, n):
        C[(1 << k, k)] = d[0][k]
    for size in range(2, n):
        for subset in combinations(range(1, n), size):
            mask = sum(1 << c for c in subset)
            for k in subset:
                prev = mask & ~(1 << k)
                C[(mask, k)], parent[(mask, k)] = min(
                    (C[(prev, m)] + d[m][k], m) for m in subset if m != k)
    full = (1 << n) - 2
    cost, last = min((C[(full, k)] + d[k][0], k) for k in range(1, n))
    tour, mask, k = [], full, last          # walk the parent pointers back
    while k != 0:
        tour.append(names[k])
        mask, k = mask & ~(1 << k), parent.get((mask, k), 0)
    return cost, names[0] + "".join(reversed(tour)), len(C)

print(held_karp(list("ABCDGH"), STOPS))
print(held_karp(list("ABCDEFGH"), STOPS))

It prints (247, 'AGDHCB', 80) for the six stops and (315, 'AGFEHDCB', 448) for all eight: the optimal cost, the tour, and the number of table entries. 315 matches the brute-force optimum from Post 2, reached with 448 table entries instead of 2,520 tours. Python dictionaries are far heavier than 4 bytes per entry, so this version runs out of memory well before the table above says it should.

Keep this function. From Part 4 on, the series builds heuristics that are fast but carry no proof, and on small instances Held-Karp is the referee: run both on the same random instances and check how often the heuristic finds the true optimum, and how far off it is when it doesn’t.

Exercise

Run held_karp on random instances of 12, 14, 16 and 18 cities. Use Python’s tracemalloc module to record peak memory for each. Before running 18, predict its peak memory from the 16-city result.

Hint: the table has (n − 1) · 2ⁿ⁻² entries. How much does that grow when n goes up by two?

What this method can’t do yet

Held-Karp turned 2,520 tours into 448 table entries for eight stops, and it gives a proof. But it fills its whole table no matter what, including entries for partial tours that are obviously terrible, like A → H → B. Memory explodes past a few dozen cities. For berlin52 the table would have about 5.7 × 10¹⁶ entries, roughly 230 petabytes at 4 bytes each, before storing a single predecessor. Even 30 cities needs 31 GB. Post 5, “Branch and bound”, takes the opposite approach: explore only the parts of the search space that could still beat the best tour found so far, and use lower bounds to throw the rest away.

References

• Bellman, R. (1962). “Dynamic Programming Treatment of the Travelling Salesman Problem.” Journal of the ACM 9(1): 61–63.
• Held, M., and Karp, R. M. (1962). “A Dynamic Programming Approach to Sequencing Problems.” Journal of the Society for Industrial and Applied Mathematics 10(1): 196–210.
• Stoian, M. (2024). “TSP Escapes the O(2ⁿn²) Curse.” arXiv:2405.03018.
• Dumas, Y., Desrosiers, J., Gelinas, E., and Solomon, M. M. (1995). “An Optimal Algorithm for the Traveling Salesman Problem with Time Windows.” Operations Research 43(2): 367–371.
• Cook, W. J. (2012). In Pursuit of the Traveling Salesman: Mathematics at the Limits of Computation. Princeton University Press.

See the whole series · ← All writing