NP-Hardness and Approximation: When Exact Is Too Slow

#Easy Problems and Hard Problems

At n = 60 the polynomial columns finish in a blink and the last column outlasts the machine. A faster computer only moves the wall: one extra item doubles the work, so a machine a thousand times faster buys only ten more items.

This lesson is about problems that seem to sit on the wrong side of that wall, and what to do when you meet one.

#P, NP and Checking an Answer

The name NP comes from "nondeterministic polynomial", an old machine model. It does not mean "not polynomial". Every problem in P is also in NP, because if you can solve a problem fast, you can check an answer fast by solving it again and comparing. So P sits inside NP.

The big open question is whether they are the same set. Does being easy to check always imply being easy to solve? Nobody has proved it either way. Most experts believe the answer is no, but belief is not proof.

from itertools import product

# Every clause must have at least one true literal.
clauses = [(1, 2, -3), (-1, 3), (-2, -3), (2, 3)]


def satisfied(clauses, assign):
    """Check a certificate. Time O(L) for L literals in the formula. Space O(1)."""
    return all(any(assign[abs(lit)] == (lit > 0) for lit in clause) for clause in clauses)


def sat_brute(clauses, n):
    """Search all 2^n assignments. Returns (assignment or None, assignments tried)."""
    tries = 0
    for bits in product([False, True], repeat=n):
        tries += 1
        assign = {i + 1: b for i, b in enumerate(bits)}
        if satisfied(clauses, assign):
            return assign, tries
    return None, tries


print(satisfied(clauses, {1: False, 2: True, 3: False}))
print(sat_brute(clauses, 3))
print(sat_brute([(1,), (-1,)], 1))
# Output:
# True
# ({1: False, 2: True, 3: False}, 3)
# (None, 2)

Checking costs one pass over the literals, so SAT is in NP. Searching costs up to 2^n checks, because each of the n variables has two choices. The last line shows the sad case: the formula "x1 and not x1" has no answer, and the search had to try everything to be sure.

#Reductions: Turning One Problem into Another

A reduction works in two directions at once, and mixing them up is the number one beginner error. If you have a fast solver for B, the recipe gives you a fast solver for A: translate, then solve. So B is at least as hard as A. If A is known to be hard, then B must be hard too, since a fast B would break A.

Hardness therefore flows along the arrow. To show your new problem is hard, you reduce a known hard problem to it, never the other way.

The pattern generalizes: S is independent exactly when everything outside S is a cover. The code checks this for every subset of a five-vertex graph, then finds the smallest cover by brute force.

from itertools import combinations

edges = [(0, 1), (1, 2), (2, 3), (3, 4), (1, 3)]
V = range(5)


def is_cover(edges, S):
    return all(u in S or v in S for u, v in edges)


def is_independent(edges, S):
    return all(not (u in S and v in S) for u, v in edges)


for r in range(6):
    for S in combinations(V, r):
        S = set(S)
        assert is_independent(edges, S) == is_cover(edges, set(V) - S)

smallest = min((set(S) for r in range(6) for S in combinations(V, r) if is_cover(edges, set(S))), key=len)
print("complement rule holds for all 32 subsets")
print(sorted(smallest), len(smallest))
# Output:
# complement rule holds for all 32 subsets
# [1, 3] 2

So "an independent set of size k exists" is the same question as "a cover of size n - k exists". The recipe is one subtraction, which is as fast as a reduction can be. The two problems are equally hard.

Give each pair of cities a distance: 1 if the graph has that edge, 2 if it does not. A round trip of total length n exists exactly when the graph has a Hamiltonian cycle, because any use of a non-edge adds at least 1 extra.

How do we know any problem is hard in the first place? In the early 1970s, Cook and Levin showed that every problem in NP reduces to SAT. After that, hardness spread by chains of reductions.

#A Short Zoo of NP-Complete Problems

Here are five you should know by name. The decision version of each is NP-complete.

All five have a certificate you can check in one pass. All five have an exact algorithm that tries every subset or order, costing about 2^n or n! steps. And all five have fast cousins that look almost the same, which is the topic of a later section.

#Subset Sum and Knapsack: Pseudo-Polynomial Time

def subset_sum_brute(nums, target):
    """Try all 2^n subsets. Time O(n * 2^n), space O(1) beyond the list."""
    ops = 0
    n = len(nums)
    for mask in range(1 << n):
        ops += 1
        if sum(nums[i] for i in range(n) if mask >> i & 1) == target:
            return True, ops
    return False, ops


def subset_sum_dp(nums, target):
    """reach[s] says whether sum s can be made. Time O(n * T), space O(T)."""
    reach = [False] * (target + 1)
    reach[0] = True
    ops = 0
    for x in nums:
        for s in range(target, x - 1, -1):    # high to low: each number used once
            ops += 1
            if reach[s - x]:
                reach[s] = True
    return reach[target], ops


nums = [3, 34, 4, 12, 5, 2]
print(subset_sum_brute(nums, 9), subset_sum_dp(nums, 9))
print(subset_sum_brute(nums, 100), subset_sum_dp(nums, 100))
for n in (10, 15, 20):
    even = list(range(2, 2 * n + 1, 2))        # all even, so an odd target is impossible
    t = 2 * n * n + 1
    print(n, subset_sum_brute(even, t)[1], subset_sum_dp(even, t)[1])
# Output:
# (True, 21) (True, 26)
# (False, 64) (False, 546)
# 10 1024 1910
# 15 32768 6540
# 20 1048576 15620

On the six numbers, brute force tried 21 of the 64 subsets before it found 4 + 5 = 9. The DP needed 26 table updates. With target 100 the picture flips: brute force must try all 64, the DP does 546 updates, and the DP wins only when T is small. The loop at the bottom shows the growth: brute force doubles with each added number, the DP grows slowly.

for digits in (3, 4, 5, 6):
    target = 10 ** (digits - 1)
    print(digits, target, subset_sum_dp([7, 11, 13, 17, 19], target)[1])
# Output:
# 3 100 438
# 4 1000 4938
# 5 10000 49938
# 6 100000 499938

This is why "knapsack in O(n * W)" and "NP-complete" are both true. Knapsack has the same DP, with the table indexed by capacity W. If W is a few thousand, the DP is excellent. If W is a 60-digit number, it is hopeless.

Some problems stay hard even when all their numbers are small. TSP is one: its distances can be 1 and 2, as in the reduction above, and it is still NP-hard.

#How to Recognise a Hard Problem

You will not meet a label that says NP-hard. You meet a story, and you have to see the shape. These signs usually mean trouble.

• You must choose a subset, an order, or an assignment from an exponentially large space, and every choice interacts with the others.
• A greedy rule or a simple DP keeps failing on small counterexamples, and your DP state would need to remember a whole set.
• The problem asks for the best of something over all subsets, orders or colorings, with a global limit such as "at most k" or "total at most W".
• It looks like a known hard problem in disguise. Scheduling with conflicts is graph coloring. Choosing teams with no feuds is independent set.

The next table is the most useful trick here: a small change in the wording often moves a problem across the wall.

Two quick tests help. First, try to prove a greedy rule works with the exchange argument from the greedy lesson. If you cannot, be suspicious. Second, check sizes: if n is at most about 20, even a brute-force 2^n method is fine.

#Plan A: Exact but Exponential

The first plan is to accept exponential time and make it as small as you can. This is the right plan when n is modest or when you need the exact optimum.

import math

for n in (5, 10, 15, 20):
    print(n, math.factorial(n - 1), n * n * 2 ** n)
# Output:
# 5 24 800
# 10 362880 102400
# 15 87178291200 7372800
# 20 121645100408832000 419430400

At n = 5 brute force wins, because the DP has overhead. By n = 10 the DP is ahead, and at n = 20 brute force is hopeless while the DP is just feasible. It is still exponential, so n = 40 is out of reach.

The search tree has depth at most k and two branches per node, so it has at most 2^(k+1) - 1 nodes. Each node costs O(m) to filter m edges. The total is O(2^k * m).

#Plan B: Approximation

#Vertex cover in two steps per edge

Scan the edges. When you meet an edge with neither endpoint in the cover so far, put both endpoints in. Otherwise skip the edge. The whole run is one pass over the edges with set lookups, so it takes O(m) time and O(n) space for the cover.

Here is the code with a ratio check against brute force. The brute force tries all subsets by increasing size and takes the first cover it finds, which is slow and exact.

from itertools import combinations


def vc_approx(edges):
    """2-approximation. Time O(m) for m edges, space O(n)."""
    cover = set()
    ops = 0
    for u, v in edges:
        ops += 1
        if u not in cover and v not in cover:
            cover.add(u)
            cover.add(v)
    return cover, ops


def vc_exact(edges):
    """Brute force. Time O(2^n * m), the answer to compare against."""
    verts = sorted({x for e in edges for x in e})
    for size in range(len(verts) + 1):
        for c in combinations(verts, size):
            if all(u in c or v in c for u, v in edges):
                return set(c)


graphs = {
    "path": [(0, 1), (1, 2), (2, 3)],
    "star": [(0, i) for i in range(1, 6)],
    "triangle": [(0, 1), (1, 2), (0, 2)],
    "mix": [(0, 1), (1, 2), (2, 3), (3, 4), (1, 3)],
}
for name, g in graphs.items():
    cover, ops = vc_approx(g)
    best = vc_exact(g)
    print(name, sorted(cover), len(best), len(cover) / len(best), ops)
# Output:
# path [0, 1, 2, 3] 2 2.0 3
# star [0, 1] 1 2.0 5
# triangle [0, 1] 2 1.0 3
# mix [0, 1, 2, 3] 2 2.0 5

The ratio never passes 2, and on the path, the star and the mixed graph it hits 2 exactly. Now the proof.

The picked edges are pairwise disjoint: both endpoints of a picked edge were new to the cover, so no later picked edge can share one. Call the picked edges M. Any cover, even the best, must contain at least one endpoint of each edge in M, and those endpoints are all different, so the best cover has at least |M| vertices. The algorithm returns exactly 2|M| vertices. So its answer is at most twice the best.

Despite decades of effort, nobody has found a polynomial algorithm whose ratio beats 2 by a fixed amount. A greedy rule that takes the vertex touching the most uncovered edges looks smarter, but its guarantee is worse: the ratio can grow with the graph.

#Greedy set cover

def greedy_set_cover(universe, sets):
    """Each round scans every set. At most len(universe) rounds, so O(rounds * sets * |universe|)."""
    left = set(universe)
    chosen = []
    while left:
        name = max(sorted(sets), key=lambda k: len(sets[k] & left))
        chosen.append(name)
        left -= sets[name]
    return chosen


U = set(range(1, 15))
sets = {"S1": set(range(1, 8)), "S2": set(range(8, 15)),
        "T1": {1, 2, 3, 4, 8, 9, 10, 11}, "T2": {5, 6, 12, 13}, "T3": {7, 14}}
print(greedy_set_cover(U, sets))
# Output:
# ['T1', 'T2', 'T3']

S1 and S2 together cover all 14 elements with 2 sets. Greedy takes T1 first because it covers 8 elements, more than the 7 of S1 or S2. After that, it needs T2 and T3 as well, so it uses 3 sets for a best of 2. Building bigger traps makes the gap grow.

The guarantee is a ratio of about ln n for n elements. Here is the idea of the proof. If the best answer uses k sets, then at any moment some set covers at least 1/k of what is left, since k sets cover all of it. Greedy covers at least that much, so each round shrinks the leftover by a factor of at least (1 - 1/k). After about k ln n rounds, nothing is left.

That is a sketch, not a full proof. It is also known that no polynomial algorithm does much better than ln n unless P equals NP, so the greedy rule is close to the best possible.

#Metric TSP in two steps

def dist(p, q):
    return abs(p[0] - q[0]) + abs(p[1] - q[1])          # grid distance: a metric


def mst_edges(pts):
    """Prim's algorithm, simple version. Time O(n^2), space O(n)."""
    inside = {0}
    best = {j: (dist(pts[0], pts[j]), 0) for j in range(1, len(pts))}
    edges = []
    while len(inside) < len(pts):
        j = min(best, key=lambda k: (best[k][0], k))
        edges.append((best[j][1], j))
        inside.add(j)
        del best[j]
        for k in best:
            d = dist(pts[j], pts[k])
            if d < best[k][0]:
                best[k] = (d, j)
    return edges


def tsp_2approx(pts):
    """MST, then first-visit order of a depth-first walk. Time O(n^2)."""
    kids = {i: [] for i in range(len(pts))}
    for a, b in mst_edges(pts):
        kids[a].append(b)
    order = []

    def walk(u):
        order.append(u)                                  # first visit only: the shortcut
        for v in kids[u]:
            walk(v)

    walk(0)
    return order


def tour_len(pts, order):
    return sum(dist(pts[order[i]], pts[order[(i + 1) % len(order)]]) for i in range(len(order)))


pts = [(0, 0), (4, 0), (4, 3), (0, 3), (2, 1), (6, 2), (1, 5), (5, 5)]
order = tsp_2approx(pts)
mst_weight = sum(dist(pts[a], pts[b]) for a, b in mst_edges(pts))
print(order)
print(tour_len(pts, order), mst_weight)
# Output:
# [0, 3, 6, 4, 1, 2, 5, 7]
# 34 21

The tour has length 34, and the MST weighs 21. The proof of the ratio has three claims.

First, the MST weighs at most the best tour: delete one edge from the best tour and you have a path through all cities, which is a spanning tree, and the MST is the lightest of those. Second, a walk along every tree edge, down and back up, has length exactly twice the MST. Third, the shortcuts never add length. So our tour is at most 2 * MST, which is at most 2 times the best tour.

The brute force check would try all (n - 1)! orders. For these eight points it finds 26 as the best, so the ratio here is 34 / 26, about 1.3, well inside the promise. Christofides' algorithm improves the ratio of metric TSP to 1.5 by being more careful about the walk.

Some problems are gentler still. For knapsack, a polynomial algorithm gets within any chosen fraction of the best value, at a cost that grows as the fraction tightens. That is a fully polynomial-time approximation scheme (FPTAS).

#Plan C: Heuristics, Local Search and Annealing

Heuristics need a safety net. Compare against brute force on small inputs, as the code did above, or run a bound beside them, so you know how far from the best you might be.

#Cautions: What Hardness Does and Does Not Say

NP-hardness is easy to over-read. Keep these five points.

A good interview script: name the known hard problem yours reduces from, give the exact method for small n, then give an approximation with its ratio.

#Where This Shows Up

• Choosing the smallest set of tests, servers or sensors that covers all requirements is set cover. Engineers often use the greedy rule, with its ln n guarantee.
• Delivery routing, drilling holes in a circuit board and ordering stops all generalize TSP. Practical tools use local search methods such as 2-opt, often after a first tour built by a simpler rule.
• Dependency resolution in package managers can be written as SAT. Some tools hand the problem to a SAT solver, which works well on real inputs even though SAT is NP-complete.
• Timetables, exam schedules and register allocation in compilers are graph coloring problems in disguise.
• In machine learning, picking the best subset of k features and finding the best k-means clustering are both NP-hard in the worst case. The standard k-means loop is a local search heuristic, and different random starts can land in different local optima.

#Common Mistakes

#Try It: Five Problems

The ladder runs from easy to adversarial. Each function returns an operation count so you can see the growth. Problem 1 checks a certificate, problems 2 and 3 are the two subset-sum methods, problem 4 checks the 2-approximation against the truth, and problem 5 is the fixed-parameter search tree.

#Recap

A polynomial cost grows gently and an exponential one hits a wall. P holds the decision problems solvable in polynomial time. NP holds those whose yes-answers can be checked in polynomial time from a certificate. P sits inside NP, and whether they are equal is open.

A reduction turns problem A into problem B, so hardness flows from A to B. NP-hard means everything in NP reduces to it, and NP-complete adds that it is in NP. SAT, subset sum, TSP, vertex cover and knapsack are all NP-complete.

The subset-sum DP costs O(n * T) but T can be exponential in the bit length of the input, so it is pseudo-polynomial. Brute force costs O(n * 2^n).

Three plans cover most real cases. Plan A is exact search: bitmask DP in O(2^n * n^2) for TSP, branch and bound, and the fixed-parameter tree for vertex cover in O(2^k * m).

Plan B is approximation: vertex cover has a 2-approximation in O(m), greedy set cover has a ratio of about ln n, and metric TSP has a 2-approximation through the MST in O(n^2). Plan C is heuristics such as 2-opt and simulated annealing, with no guarantee.