Algorithms · Dynamic Programming

Subset DP with Bitmask States: Held-Karp and the 2-to-the-n Wall

Encoding a visited set as an integer makes the traveling-salesman recurrence exact in O(2^n n^2), and the same encoding supports submask enumeration and the sum-over-subsets transform.

A bitmask DP exists when the state needs to record a subset of a small universe and nothing else. The subset becomes an integer, the table is indexed by that integer, and the recurrence is a loop over the elements not yet in it. It is the extreme case of the extra-dimension trade described in state sufficiency: the dimension is a set, and its size is exponential in the universe. The technique is worth studying separately from ordinary tabulation because its cost model is different in kind: the table size is 2^n for a universe of n elements, so the practical limit is set by memory bandwidth rather than by instruction count, and the boundary between feasible and infeasible is sharp.

The canonical instance is the traveling salesman problem solved exactly by the Held-Karp recurrence, published in 1962 and still the best known general exact method: O(2^n n²) time and O(2^n n) space, against O(n!) for enumeration.

The state and the recurrence

Fix a starting vertex 0. Let dp[S][v] be the minimum cost of a path that starts at 0, visits exactly the vertices in S, and ends at v, where v ∈ S and 0 ∈ S. Then

dp[{0}][0] = 0
dp[S][v]   = min over u in S \ {v} of ( dp[S \ {v}][u] + w(u, v) )
answer     = min over v != 0 of ( dp[full][v] + w(v, 0) )

The set S is sufficient state, by the argument in state sufficiency: the only constraint on extending a path is which vertices are already used, and the cost of the extension depends on the endpoint. Nothing else about the path matters. Every component contains vertex 0, so there are 2^{n-1} sets rather than 2^n, and at most n endpoints per set: the number of states is 2^{n-1} · n, and each state scans up to n predecessors, so the time is O(2^n n²) with a constant of one half.

Encoding S as a uint32 or uint64 is the standard choice, with bit i set when vertex i is in S. The operations needed are membership (S >> v & 1), insertion (S | 1 << v), removal (S & ~(1 << v)), and iteration over members. For n ≤ 64 a single machine word holds the set, and the table is a flat array of size 2^n · n with index S * n + v. For n > 64 the set needs multiple words and the table size has already passed what any machine can hold, so the restriction costs nothing in practice.

Ordering: clearing a bit always lowers the mask

Every transition goes from S \ {v} to S, so the state graph is layered by |S|. Iterating masks in increasing order of their integer value happens to be a valid topological order, because clearing a bit always decreases the value, and a decrease in value implies the source mask comes first. Numeric order therefore refines the layering by cardinality without having to compute it. This is an accident of the binary encoding, not a general property of subset DPs.

Grouping by popcount is the alternative and is what makes the layered structure explicit. It is required when the transition can keep the set size constant or increase it by more than one, as in some matching and covering formulations. For Held-Karp, ascending integer order suffices and avoids the cost of bucketing masks by popcount.

def tsp(dist):
    n = len(dist)
    full = (1 << n) - 1
    INF = float('inf')
    dp = [[INF] * n for _ in range(1 << n)]
    dp[1][0] = 0                                  # only vertex 0 visited, standing at 0
    for mask in range(1, full + 1):
        if not mask & 1:                          # every path starts at 0
            continue
        for v in range(n):
            if not mask >> v & 1:
                continue
            base = dp[mask][v]
            if base == INF:
                continue
            for u in range(n):
                if mask >> u & 1:
                    continue
                nm = mask | (1 << u)
                cand = base + dist[v][u]
                if cand < dp[nm][u]:
                    dp[nm][u] = cand
    return min(dp[full][v] + dist[v][0] for v in range(1, n))

Two structural notes. The push form used above — extend from a complete state to a new state — is equivalent to the pull form in the recurrence and is usually faster, because it skips unreachable states early and writes to a contiguous target block. And the mask & 1 filter halves the table traffic: no state without vertex 0 is ever reachable.

The memory wall, quantified

The table has 2^n · n entries. At four bytes each:

nentriesbytes
154.9×10^52 MB
202.1×10^784 MB
229.2×10^7370 MB
258.4×10^83.4 GB
287.5×10^930 GB

The progress in each column is a factor of four per added vertex, in entries and in bytes alike; the bytes column runs slightly under that because of rounding. Restricting every mask to contain vertex 0 halves both columns throughout, so a table that already applies that filter reaches one vertex further than the figures suggest. The practical ceiling on a commodity machine is around n = 20 to n = 22, and it is a memory limit, not a time limit: time grows by the same factor of four per vertex, so the two escalate together, and it is the table that has to be resident and written that runs out first. Storing dp as uint16 when distances fit in 16 bits halves the footprint, which buys about half a vertex of additional reach: halving the bytes offsets only part of the factor of four that each added vertex costs. Dropping the last layer is not possible — every layer is consulted by the next.

The Held-Karp lower bound used inside branch-and-bound solvers is a different object: it comes from the solution of a 1-tree relaxation, not from the subset recurrence above, which is why the name appears in two contexts.

Enumerating submasks

A frequent sub-step is iterating over all subsets of a given mask. The idiom

sub = mask
while sub:
    # ... use sub ...
    sub = (sub - 1) & mask

visits every non-empty submask exactly once, in decreasing numeric order. The step (sub - 1) & mask clears the lowest set bit of sub and then refills with all bits of mask below it — the standard trick for counting down through a masked set without leaving it. Termination is guaranteed because the value strictly decreases and is bounded below by zero.

The total cost of enumerating submasks for every mask is Σ_{S} 2^{|S|} = 3^n, since each element is in the outer set, the inner set, or neither. 3^n is a worse base than 2^n but a much better one than 4^n, and it is the complexity of several subset-partition DPs, including the standard exact solution for graph coloring and for partitioning a set into the fewest valid groups.

The sum-over-subsets transform

A different pattern uses the same encoding for a transform rather than a recurrence. Given f[S] for all 2^n subsets, compute g[S] = Σ_{T ⊆ S} f[T]. The direct computation is 3^n. The zeta transform does it in n · 2^n:

def subset_sums(f, n):
    g = f[:]                       # g[S] accumulates over subsets of S
    for bit in range(n):
        step = 1 << bit
        for mask in range(1 << n):
            if mask & step:
                g[mask] += g[mask ^ step]
    return g

The invariant after processing bit i is that g[S] sums f[T] over all T that agree with S on bits above i and are arbitrary below. Each bit is added exactly once, so at the end the sum covers all subsets. The same loop with the direction reversed computes sums over supersets, and the same structure with subtraction instead of addition inverts the transform. This is the subset lattice’s analogue of the fast Walsh-Hadamard and fast Möbius transforms, and it is the workhorse behind counting problems phrased as “how many subsets satisfy a property defined per element”.

The memory access pattern is worth noting: the inner loop strides by step, which for small bit is contiguous and for large bit jumps across the array. The transform is memory-bound at n · 2^n read-modify-writes, and its practical speed is dominated by cache behaviour rather than by the arithmetic.

Where the exponential cannot be avoided

Bitmask DP is exact and its cost is exponential in the universe size. Three responses exist, and choosing among them is a modeling decision rather than an optimization.

Restrict the universe. If the problem’s structure limits which elements can be combined — a graph of bounded pathwidth, a set of items grouped into small independent components — the DP can be run per component and combined, and the total cost is the product of the component costs rather than 2^(total).

Meet in the middle. Split the universe into halves of size n/2, enumerate all subsets of each half, and join. This gives O(2^{n/2}) for problems where the two halves interact through a small interface, at the cost of a matching step between the halves. For subset-sum, sorting one half’s sums and binary-searching the other is the standard form.

Approximate. For metric instances of the TSP, the minimum spanning tree doubled gives a 2-approximation and the Christofides construction gives 3/2, both polynomial. Exactness and polynomial time are not simultaneously available unless P = NP.

The exponential base is the algorithm’s identity, not a defect to be engineered around. What can be engineered is the constant: the table’s layout, the width of the distance values, the early pruning of unreachable masks, and the choice of push versus pull transitions. On a modern machine those choices account for a factor of several, while the base accounts for everything else.

Two implementation hazards

The first is integer width. For n = 31 in C, the mask (1 << n) - 1 does not fit a signed 32-bit int, and shifting a signed value into its sign bit is undefined behaviour; the mask and the table size must be unsigned, and the literal must be written as 1u << n so the shift happens in the unsigned type. The failure mode is a table allocated with a negative or absurd size, or a loop bound that never terminates.

The second is that the answer’s correctness depends on the problem being a simple path. Held-Karp computes the minimum over paths that visit each vertex exactly once. That is the right question for the TSP, whose feasible tours are Hamiltonian cycles by definition, but it is not the shortest walk: if edge weights can be negative, a walk that repeats a vertex can be cheaper than every Hamiltonian cycle, and the recurrence will never consider it. Put differently, the recurrence solves the shortest-Hamiltonian-cycle problem exactly and says nothing about the shortest-path problem, and which one the caller needs is a property of the problem statement, not of the code.

References

  1. A Dynamic Programming Approach to Sequencing Problems — Michael Held and Richard Karp, Journal of the Society for Industrial and Applied Mathematics 10(1), 1962
  2. The Art of Computer Programming, Volume 4A: Combinatorial Algorithms, Part 1 — Section 7.1.3 — Donald Knuth; Addison-Wesley
  3. Introduction to Algorithms, 4th edition — Chapter 15, Dynamic Programming — Cormen, Leiserson, Rivest, Stein; MIT Press