TSP via Dynamic Programming, in depth: implementing Held-Karp with tour recovery, memory sizing, numpy layers and constraints

The recurrence, precisely

Number the cities 0 to n-1 and fix city 0 as the start; a tour is a cycle, so any start gives the same optimum. Let S be a subset of {1, ..., n-1} and j a city in S. Define D(S, j) as the length of the shortest path that starts at 0, visits every city in S exactly once, and ends at j. Then:

D({j}, j)  = w(0, j)                                   for every j in 1..n-1
D(S, j)    = min over i in S - {j} of  D(S - {j}, i) + w(i, j)     for |S| >= 2
answer     = min over j in 1..n-1 of  D({1..n-1}, j) + w(j, 0)

The recurrence is correct because the best path to j through S must arrive from some last city i, and the part before i is itself a best path through S without j that ends at i. That is the optimal-substructure property every DP needs, and it holds because path length is a sum of edge weights. The order of the visits inside S - {j} does not matter for what comes next, which is why the state can forget it. Nothing assumes symmetry or the triangle inequality, so the same code solves the asymmetric problem.

Fixing city 0 outside the mask is not just tidiness. It halves the table compared with masks over all n cities, and it removes the classic bug of a base case D({0}, 0) = 0 combined with masks that may or may not contain 0.

Worked example: four cities

Take four cities with an asymmetric cost matrix, where w(a,b) is the cost from a to b (think one-way streets or uphill legs). Row by row: from 0 the costs to 1, 2, 3 are 4, 9, 6; from 1 to 0, 2, 3 they are 7, 3, 8; from 2 to 0, 1, 3 they are 5, 6, 2; from 3 to 0, 1, 2 they are 3, 9, 4. Base cases: D({1},1) = 4, D({2},2) = 9, D({3},3) = 6.

Size two has only one choice per cell: D({1,2},2) = D({1},1) + w(1,2) = 7, D({1,2},1) = 9 + 6 = 15, D({1,3},3) = 4 + 8 = 12, D({1,3},1) = 6 + 9 = 15, D({2,3},3) = 9 + 2 = 11 and D({2,3},2) = 6 + 4 = 10.

Size three: D({1,2,3},1) = min(D({2,3},2) + w(2,1), D({2,3},3) + w(3,1)) = min(16, 20) = 16. D({1,2,3},2) = min(D({1,3},1) + 3, D({1,3},3) + 4) = min(18, 16) = 16. D({1,2,3},3) = min(D({1,2},1) + 8, D({1,2},2) + 2) = min(23, 9) = 9.

Closing: min(16 + 7, 16 + 5, 9 + 3) = 12. Following the choices back from j=3 gives 3 from 2 from 1 from 0, so the tour is 0, 1, 2, 3, 0 with cost 4 + 3 + 2 + 3 = 12. Driving it the other way, 0, 3, 2, 1, 0, costs 23, which is why the asymmetric case must never assume a tour and its reverse are equal.

Implementation with tour recovery

A clean implementation stores masks over cities 1..n-1, so bit b of the mask stands for city b+1. Masks are processed in increasing numeric order, which is enough: removing a bit always gives a smaller number, so every dependency is already computed. A parent table records the argmin so the tour can be rebuilt.

import math

def held_karp(w):
    n = len(w)
    if n == 1:
        return 0, [0, 0]
    m = n - 1                                   # cities 1..n-1 as bits 0..m-1
    FULL = (1 << m) - 1
    INF = math.inf
    D = [[INF] * m for _ in range(1 << m)]
    P = [[-1] * m for _ in range(1 << m)]
    for j in range(m):
        D[1 << j][j] = w[0][j + 1]
    for S in range(1, FULL + 1):
        for j in range(m):
            if not (S >> j) & 1 or S == 1 << j:
                continue
            prev = S ^ (1 << j)
            best, arg = INF, -1
            for i in range(m):
                if (prev >> i) & 1:
                    v = D[prev][i] + w[i + 1][j + 1]
                    if v < best:
                        best, arg = v, i
            D[S][j], P[S][j] = best, arg
    best, last = min((D[FULL][j] + w[j + 1][0], j) for j in range(m))
    tour, S, j = [], FULL, last
    while j != -1:                              # walk parents back to the start
        tour.append(j + 1)
        S, j = S ^ (1 << j), P[S][j]
    return best, [0] + tour[::-1] + [0]

On the example matrix it returns 12 and [0, 1, 2, 3, 0], exactly the tour found by hand. Pure Python manages about n = 15 to 17 in seconds; beyond that use the vectorised form below or a compiled language.

Sizing time and memory

Size the table before running. With the start fixed there are 2^n-1 masks and n-1 end cities, so (n-1) 2^n-1 cells, and about (n-1)² 2^n-1 inner-loop steps. Storing distances as 4-byte floats and parents as 1-byte integers gives 5 bytes per cell:

So n = 20 is easy anywhere, n = 25 fits on a workstation and takes minutes in compiled code, and n around 28 to 30 needs a large-memory server and patience. Pure Python's per-cell object overhead is tens of bytes rather than five, which is why it hits the wall far earlier. Memory, not time, is usually the first limit. Each layer depends only on the previous one, so with care you can keep just two layers of distances, but you then need to store parents or recompute them to rebuild the tour.

A vectorised version

The inner loop vectorises well if you process subsets layer by layer, grouping masks by how many bits they have. Within a layer every mask depends only on the previous layer, so you can update all of them with one array operation per (i, j) pair:

import numpy as np

def held_karp_np(w):
    w = np.asarray(w, dtype=np.float64)
    m = len(w) - 1
    masks = np.arange(1 << m)
    pop = np.array([bin(x).count("1") for x in masks])
    D = np.full((1 << m, m), np.inf)
    for j in range(m):
        D[1 << j, j] = w[0, j + 1]
    for k in range(2, m + 1):
        layer = masks[pop == k]
        for j in range(m):
            S = layer[((layer >> j) & 1) == 1]
            prev = S ^ (1 << j)
            best = np.full(len(S), np.inf)
            for i in range(m):
                if i != j:
                    best = np.minimum(best, D[prev, i] + w[i + 1, j + 1])
            D[S, j] = best
    return float(np.min(D[-1] + w[1:, 0]))

This version returns only the cost; add an argmin array per layer if you need the tour. The Python loop runs about n³ times while numpy handles the 2ⁿ factor, so it reaches n = 20 to 22 comfortably. With D laid out as [mask][city], the column read D[prev, i] is a strided gather across many rows; storing D as [city][mask] turns it into a gather within one row, which can be kinder to caches. Measure before committing to either layout.

Variants and constraints

The state (S, j) is flexible, which is why Held-Karp survives in practice for small, heavily constrained problems where heuristics struggle to stay feasible.

• Open path, fixed start. Drop the w(j, 0) term in the answer.
• Open path, any start. Add a dummy city with zero-cost edges to and from everything and solve the cycle; the dummy's two edges mark the ends. This is the shortest Hamiltonian path, which backtracking can only find by search.
• Precedence constraints (pick up before drop off). Skip any transition to j while j's predecessor is not in S. Infeasible states are never filled, which also saves time.
• Time windows. Let D hold the earliest arrival time at j. If waiting is allowed and travel times are fixed, arriving earlier never hurts, so the minimum is still the right value to keep. Transitions that miss j's window are dropped. With costs that differ from time, you need a Pareto set per state instead of one number.
• Several salesmen or capacity limits. These usually break the single-number state and are better solved with integer programming.

Testing against brute force

Exact algorithms are only valuable if they are exactly right, so cross-check against brute force on random small instances, both symmetric and asymmetric, and verify that the returned tour has the returned cost:

import itertools, random

def brute(w):
    n = len(w)
    return min(sum(w[a][b] for a, b in zip((0,) + p, p + (0,)))
               for p in itertools.permutations(range(1, n)))

for trial in range(500):
    n = random.randint(2, 8)
    w = [[0 if i == j else random.randint(1, 100) for j in range(n)] for i in range(n)]
    cost, tour = held_karp(w)
    assert cost == brute(w), (w, cost)
    assert sorted(tour[:-1]) == list(range(n)) and tour[0] == tour[-1] == 0
    assert cost == sum(w[a][b] for a, b in zip(tour, tour[1:]))

Include n = 2 and n = 3 in the range, since off-by-one errors in the base case show up there. Use integer weights in tests so equality is exact.

The other Held-Karp, and bigger instances

One naming trap: the Held-Karp lower bound is a different thing. It comes from the same authors' 1970 and 1971 papers on minimum 1-trees with Lagrangian node penalties, and gives a strong lower bound on tour length for branch-and-bound. When a paper says it uses Held-Karp, check which one it means.

For large instances neither is the main tool. Exact solvers such as Concorde use linear programming with cutting planes inside branch-and-bound and have solved instances with tens of thousands of cities. In day-to-day work, local search heuristics such as 2-opt and Lin-Kernighan get within a few percent quickly, and approximation algorithms give guarantees for metric instances. Held-Karp's role is exact answers for small n, constrained sub-problems inside a larger solver, and ground truth for testing heuristics.

Failure modes

• Base case off by one: D({j}, j) must be w(0, j), not 0.
• Forgetting the return edge, which silently solves the open-path problem.
• Iterating subsets by size incorrectly in hand-written enumeration, reading cells not yet filled. Increasing numeric order is always safe.
• Float infinity mixed with int dtypes in numpy, which raises or overflows; use float arrays or a large sentinel.
• Allocating 2ⁿ x n Python lists for n above about 20 and swapping the machine to death.

What to do next

1. Implement held_karp and run the brute-force test until 500 trials pass.
2. Compute your own table size for your target n and dtype before running anything bigger.
3. Try the numpy version and time it against the pure Python one at n = 14, 16, 18.
4. Add one constraint (precedence or time windows) and extend the tests to cover it.
5. Compare a 2-opt heuristic against the exact answer on 100 random instances and record the gap.
6. Review the dynamic programming introduction for how this state design generalises to other subset DPs.