Integer Programming Hardness, in depth: the 3-SAT reduction, why IP is in NP, strong NP-hardness and the structures that make it easy

The problem, stated precisely

An integer program asks for an integer vector x minimizing c·x subject to Ax ≥ b. Complexity is defined for the decision version: given A, b, c and a target k, is there an integer x with Ax ≥ b and c·x ≤ k? The target can be folded into the constraints as one more row, so the decision problem is really feasibility: does the polyhedron contain an integer point? A 0-1 integer program restricts each variable to 0 or 1. A mixed-integer program allows some continuous variables, which never makes the worst case easier.

Input size matters for every statement below. Coefficients are written in binary, so a number like 1,000,000 costs about 20 bits, not a million units. An algorithm whose time grows with the magnitude of the numbers rather than their bit length is not polynomial, and that distinction is the whole difference between knapsack and 0-1 integer programming in general.

The complexity landscape

NP-completeness by reduction from 3-SAT

0-1 integer programming was on Karp's 1972 list of 21 NP-complete problems. Membership in NP is immediate for the 0-1 case: a candidate x is a list of bits, and checking Ax ≥ b takes polynomial time. Hardness comes from a reduction from 3-SAT, which is already known to be NP-complete (see the 3-SAT reduction). Each Boolean variable becomes a 0-1 variable. Each clause becomes one linear row: a positive literal x_i contributes x_i, a negated literal contributes 1 − x_i, and the clause is satisfied exactly when the sum is at least 1. Moving the constants to the right-hand side gives an integer row.

def sat_to_ip(n_vars, clauses):
    """3-CNF -> 0-1 feasibility. Literal +i means x_i, -i means not x_i.
    Returns rows (a, b) meaning a . x >= b."""
    rows = []
    for clause in clauses:
        a, b = [0] * n_vars, 1
        for lit in clause:
            v = abs(lit) - 1
            if lit > 0:
                a[v] += 1
            else:                      # (1 - x_v) contributes -x_v and moves 1 right
                a[v] -= 1
                b -= 1
        rows.append((a, b))
    return rows

Take the formula (x1 ∨ x2 ∨ x3) ∧ (¬x1 ∨ ¬x2 ∨ x3) ∧ (x1 ∨ ¬x3 ∨ x4) ∧ (¬x2 ∨ ¬x3 ∨ ¬x4). The rows are x1 + x2 + x3 ≥ 1, −x1 − x2 + x3 ≥ −1, x1 − x3 + x4 ≥ 0 and −x2 − x3 − x4 ≥ −2. Searching 0-1 vectors in order, the all-zero vector fails the first row, (0,0,0,1) fails it too, (0,0,1,0) fails the third, and (0,0,1,1) satisfies all four; setting x3 and x4 true satisfies every clause of the formula, as the reduction promises.

The reduction runs in linear time and maps satisfying assignments to feasible points and back, so a polynomial algorithm for 0-1 feasibility would solve 3-SAT. Notice what it produces: every coefficient is 0, 1 or −1 and every right-hand side is between −2 and 1. That detail matters two sections from now.

Testing the reduction against brute force

Reductions are easy to get subtly wrong, usually in the sign handling of negated literals. Test one the way you would test any transformation: against brute force on many small random instances, checking both directions and the witness.

import itertools, random

def ip_feasible(n, rows):
    for x in itertools.product((0, 1), repeat=n):
        if all(sum(ai * xi for ai, xi in zip(a, x)) >= b for a, b in rows):
            return x
    return None

def sat_brute(n, clauses):
    for x in itertools.product((False, True), repeat=n):
        if all(any(x[abs(l) - 1] == (l > 0) for l in cl) for cl in clauses):
            return x
    return None

random.seed(1)
for _ in range(2000):
    n = random.randint(3, 8)
    clauses = [tuple(random.choice((1, -1)) * v
                     for v in random.sample(range(1, n + 1), 3))
               for _ in range(random.randint(1, 6 * n))]
    x = ip_feasible(n, sat_to_ip(n, clauses))
    assert (x is None) == (sat_brute(n, clauses) is None)
    if x is not None:              # the IP point must satisfy the formula
        assert all(any((x[abs(l) - 1] == 1) == (l > 0) for l in cl) for cl in clauses)

On 2,000 random formulas with 3 to 8 variables, 1,752 of them satisfiable, the two answers agreed every time and every feasible point decoded to a satisfying assignment.

Why general IP is in NP at all

For 0-1 programs the certificate is short by construction. For general integer programs it is not obvious that a feasible system has any solution you could even write down in polynomial space: the polyhedron may be unbounded, and its integer points may all be astronomically large. If the smallest solution needed exponentially many bits, a yes-answer would have no short certificate and the problem would not obviously be in NP.

Borosh and Treybig, and Papadimitriou in 1981, showed that this does not happen: if Ax ≥ b has an integer solution, it has one whose entries are bounded by a quantity whose bit length is polynomial in the size of A and b. The proof uses the fact that the vertices of the polyhedron are ratios of determinants of submatrices, which are bounded by the coefficients, plus a bounded integer step from a vertex. So general integer programming is NP-complete, not merely NP-hard. The practical echo is that solvers benefit from explicit bounds on every integer variable; an unbounded integer variable is both a modelling smell and a source of numerical trouble.

Strong versus weak NP-hardness

The 0-1 knapsack problem is also NP-hard, yet a dynamic program solves it in time proportional to n times the capacity W. That is polynomial in the magnitude of W, not its bit length, so it is pseudo-polynomial: fast when the numbers are small, exponential when they are written with many digits. Problems with such algorithms are called weakly NP-hard.

0-1 integer programming is strongly NP-hard: it stays NP-hard even when every number in the input is bounded by a polynomial in the input size. The 3-SAT reduction proves it, because its coefficients are 0 and ±1. Two consequences follow, both under the assumption that P differs from NP. There is no pseudo-polynomial algorithm, so no knapsack-style dynamic program over the right-hand sides can rescue the general case. And problems that are strongly NP-hard with polynomially bounded objective values have no fully polynomial-time approximation scheme. When someone proposes a DP over capacities for a general model, the number of constraint rows is what blows up: the state space is the product of all the right-hand sides.

Why approximation guarantees fail

For many NP-hard optimization problems we settle for approximation: a fast algorithm guaranteed to land within a factor of the optimum. General integer programming offers no such guarantee, for a blunt reason: deciding whether any feasible point exists is NP-complete, so no polynomial algorithm can even promise to return a feasible solution, let alone a good one, unless P equals NP.

Guarantees reappear only for specific problem families whose feasibility is trivial, and even there the limits are tight. Set cover, a 0-1 program where taking every set is always feasible, can be approximated within a factor of about ln n by greedy, and Dinur and Steurer showed in 2014 that doing meaningfully better is NP-hard. Techniques that achieve such guarantees mostly go through the LP relaxation; LP rounding covers how integrality gaps bound them.

When integer programs are easy

Hardness is a worst-case statement about the whole class. Some subclasses are polynomial, and recognizing them in your model is often worth more than any solver setting.

Total unimodularity. A matrix is totally unimodular (TU) if every square submatrix has determinant 0, 1 or −1. Hoffman and Kruskal showed that for an integer matrix A, the polyhedron {x ≥ 0 : Ax ≤ b} has integer vertices for every integer b exactly when A is TU. Then the LP relaxation, solvable in polynomial time, already returns an integer optimum. Incidence matrices of bipartite graphs, node-arc incidence matrices of directed graphs (network flow) and interval matrices with consecutive ones are TU, which is why assignment, flow and many interval scheduling models solve instantly. Checking TU by brute force is exponential but fine for spotting the structure in a small slice of a model:

import itertools
from fractions import Fraction

def det(M):
    M = [[Fraction(v) for v in row] for row in M]
    n, d = len(M), Fraction(1)
    for col in range(n):
        piv = next((r for r in range(col, n) if M[r][col] != 0), None)
        if piv is None:
            return 0
        if piv != col:
            M[col], M[piv] = M[piv], M[col]
            d = -d
        d *= M[col][col]
        for r in range(col + 1, n):
            f = M[r][col] / M[col][col]
            for k in range(col, n):
                M[r][k] -= f * M[col][k]
    return int(d)

def is_tu(A):
    """Exponential check: every square submatrix has det in {-1, 0, 1}."""
    m, n = len(A), len(A[0])
    for k in range(1, min(m, n) + 1):
        for rows in itertools.combinations(range(m), k):
            for cols in itertools.combinations(range(n), k):
                sub = [[A[r][c] for c in cols] for r in rows]
                if det(sub) not in (-1, 0, 1):
                    return False, rows, cols
    return True, None, None

The vertex-edge incidence matrix of the complete bipartite graph on two plus two vertices passes. The incidence matrix of a triangle fails with determinant 2 on the full matrix, and that odd cycle is exactly why matching on a triangle has an LP optimum of 1.5, every edge at one half, while the best integer matching is 1.

Fixed dimension. Lenstra showed in 1983 that integer programming with a fixed number of variables is solvable in time polynomial in the input size; later algorithms improved the dependence on the dimension, which remains exponential. Small numbers and few rows. With a constant number of constraints and small coefficients, dynamic programming over right-hand sides is pseudo-polynomial, as in knapsack and its multi-dimensional variant.

What hardness looks like in a solver

Worst-case hardness says some instances of the class are hard; it does not say yours is. In practice hardness shows up in recognizable ways. Weak LP relaxations leave a large gap that branch and bound must close by enumeration. Symmetry, such as interchangeable machines or identical trucks, multiplies equivalent subtrees. Equality constraints with large coefficients, the integer analogue of knapsack equations, defeat both bounding and rounding. And random instances near a phase transition, like random 3-SAT near about 4.27 clauses per variable, are hard for every known method.

When a model stalls, read the solver log before tuning anything: if the best bound barely moves while the incumbent improves, the formulation is weak; if both stall and the node count explodes, look for symmetry; if the solver struggles to find any feasible point, the feasibility problem itself is hard and a heuristic start or a relaxation of hard constraints into penalties may be needed. Always run with a time limit and report the optimality gap with the answer.

Failure modes

• Reducing in the wrong direction. Showing your problem can be written as an IP proves nothing about its hardness; you must reduce a known hard problem to yours.
• Calling a pseudo-polynomial algorithm polynomial. A DP over capacities is fast until someone scales the units to cents and the table grows by 100 times.
• Treating NP-hard as hopeless. Structured instances with tight formulations solve routinely; check for TU or near-TU structure first.
• Leaving integer variables unbounded. Bounds are free information for the solver and guard against numerical blow-up.
• Expecting an approximation ratio from a heuristic on a general model. Without guaranteed feasibility there is none; report the gap instead.

Trade-offs

Exact solution with branch and bound gives proven optimality at unpredictable cost. Exploiting structure, such as TU or decomposition into flow plus side constraints, gives speed at the price of modelling effort and fragility when requirements change. Problem-specific approximation gives guarantees but only for the family it was built for. Heuristics and time-limited solves give answers on time with a measured gap rather than a guarantee. Most production systems combine them: a structured model, a heuristic warm start, and a solver with a time limit whose gap is logged.

What to do next

1. Write down the decision version of your model and identify which variables are 0-1, general integer or continuous.
2. Look for a known hard problem embedded in it; if 3-SAT, set cover or knapsack fits, expect hardness.
3. Check slices of the constraint matrix for total unimodularity before adding cuts or heuristics.
4. Put explicit, tight bounds on every integer variable.
5. If coefficients are small and rows few, try a dynamic program and compare it with the solver.
6. Run the solver with a time limit and log bound, incumbent, gap and node count per run.
7. Test any reduction or reformulation against brute force on small random instances.
8. When the gap stays large, strengthen the formulation before buying a bigger machine.