LP Rounding for Approximation, in depth: integrality gaps, threshold rounding, extreme points, unrelated-machine scheduling and iterative rounding

Relax, solve, round

Every LP rounding algorithm has the same four steps, and the analysis lives in the gap between them.

1. Formulate. An integer program whose optimum OPT is the answer you want.
2. Relax. Replace x in {0, 1} with 0 ≤ x ≤ 1. The LP optimum, LP*, satisfies LP* ≤ OPT for minimisation.
3. Solve. Polynomial in theory, fast in practice with HiGHS, Gurobi or similar.
4. Round. Produce an integral solution of cost at most ρ times LP*. Then cost ≤ ρ LP* ≤ ρ OPT.

The rounding step must preserve feasibility, which is usually the hard part, and it must lose a bounded factor of cost. Each technique below is a different way of reading enough structure from x to do both.

The integrality gap is the ceiling

The integrality gap is the worst ratio OPT / LP* over all instances of a formulation. No rounding analysed against LP* can beat it, because on the worst instance even a perfect algorithm pays the gap. That makes it the first number to compute when choosing a formulation.

For vertex cover, the complete graph on n vertices gives LP* = n/2 with every x = 1/2, while OPT = n - 1, so the gap tends to 2. For set cover the gap is logarithmic in the number of elements, which is why no LP rounding beats roughly ln n there. When a gap is too large, the remedy is a stronger formulation, with extra valid inequalities that cut off the bad fractional points, not a cleverer rounding of the weak one.

Threshold rounding

Threshold rounding is the simplest deterministic scheme. In set cover where every element lies in at most f sets, each element's constraint sums at most f variables to at least 1, so at least one of them is at least 1/f. Keep every set with x ≥ 1/f: every element is covered, and since each kept set had x ≥ 1/f, its cost is at most f times its LP contribution. The result is an f-approximation, and vertex cover is the case f = 2 with threshold 1/2.

import numpy as np
from scipy.optimize import linprog

def threshold_cover(costs, sets, n_elems):
    """Set cover by LP threshold rounding. If every element lies in at most f
    sets, keeping each set with x >= 1/f costs at most f * LP*."""
    m = len(sets)
    A = np.zeros((n_elems, m))
    for j, s in enumerate(sets):
        for e in s:
            A[e, j] = 1.0
    f = int(A.sum(axis=1).max())
    res = linprog(costs, A_ub=-A, b_ub=-np.ones(n_elems),
                  bounds=[(0, 1)] * m, method="highs")
    chosen = [j for j in range(m) if res.x[j] >= 1.0 / f - 1e-9]
    return res.fun, res.x, chosen, f

def prune(chosen, costs, sets, n_elems):
    """Drop redundant sets, most expensive first. Never raises cost."""
    keep = sorted(chosen, key=lambda j: -costs[j])
    for j in list(keep):
        rest = [k for k in keep if k != j]
        if set().union(*(sets[k] for k in rest)) == set(range(n_elems)):
            keep = rest
    return keep

The prune pass is not in the proof and costs nothing in the bound, but in practice it is where most of the quality comes from, as the next example shows. The 1e-9 tolerance matters: solvers return 0.4999999999 for 1/2, and an exact comparison silently drops a set and breaks feasibility.

Worked example: the five-cycle

Take the five-cycle as a vertex cover instance with unit costs: five sets, one per vertex, each covering its two incident edges, so f = 2. Running the code gives LP* = 2.5 with every x = 0.5. Threshold rounding keeps all five vertices, cost 5, exactly twice LP*, so the factor-2 analysis is tight on this instance. OPT is 3, an odd cycle needs three vertices. The prune pass removes two redundant vertices and returns a cover of cost 3, which is optimal.

Two lessons follow. The worst case of the analysis is reached on a tiny, symmetric instance, and symmetry is exactly what makes LP solutions fractional. And the reported guarantee and the delivered cost can differ by the full factor, so always report both the cost and LP*: the ratio cost / LP* is an instance-specific certificate of quality that is usually far better than the worst-case bound.

Extreme points carry structure

Threshold rounding ignores where the LP solution lies. Extreme-point rounding uses it. A basic feasible solution, a vertex of the polytope, is determined by tight constraints, and a vertex of an LP with k constraints other than variable bounds has at most k variables strictly between their bounds; with only x ≥ 0, that means at most k non-zero variables. That counting argument limits how fractional a vertex can be, and for structured problems the fractional part has a shape you can round exactly.

Vertex cover is the cleanest case: every vertex of its LP is half-integral, with values in {0, 1/2, 1}, a result due to Nemhauser and Trotter; the vertex cover article uses it. The bigger payoff is scheduling. Note the practical requirement: you need a vertex, which simplex methods return, while interior-point methods return a point in the middle of the optimal face unless you run crossover. In SciPy, ask for method="highs-ds".

Scheduling on unrelated machines

Scheduling n jobs on m unrelated machines, where job j takes p_ij on machine i, to minimise the makespan is NP-hard, and Lenstra, Shmoys and Tardos showed in 1990 that no polynomial algorithm achieves a ratio below 3/2 unless P = NP. Their algorithm gets 2.

The plain LP is useless: it can split a long job across machines and report a makespan far below OPT. The fix is to guess a target T and remove every pair with p_ij > T, because no schedule with makespan T can use them. That turns the optimisation into a feasibility LP, and bisection finds the smallest feasible T, which is at most OPT.

Now count. The LP has n + m constraints besides bounds, so a vertex has at most n + m non-zero variables. Every job has at least one, so at most m jobs are split. The finer argument shows the graph joining split jobs to their machines is a pseudoforest, each component having at most one cycle, so it has a matching that covers every split job. Assign integral jobs as the LP says, at most T per machine, and give each split job its matched machine. Each machine gains at most one job, which costs at most T because the pair survived pruning. Makespan ≤ 2T ≤ 2 OPT.

def lp_schedule(p, T):
    """Feasibility LP for makespan T; pairs with p[i, j] > T are removed."""
    m, n = p.shape
    pairs = [(i, j) for i in range(m) for j in range(n) if p[i, j] <= T]
    if any(all(p[i, j] > T for i in range(m)) for j in range(n)):
        return None, pairs
    k = len(pairs)
    A_eq, A_ub = np.zeros((n, k)), np.zeros((m, k))
    for c, (i, j) in enumerate(pairs):
        A_eq[j, c] = 1.0          # each job fully assigned
        A_ub[i, c] = p[i, j]      # each machine loaded at most T
    res = linprog(np.zeros(k), A_ub=A_ub, b_ub=np.full(m, T), A_eq=A_eq,
                  b_eq=np.ones(n), bounds=[(0, None)] * k, method="highs-ds")
    return (res.x if res.status == 0 else None), pairs

def match(frac_jobs, edges):
    """Give every fractional job its own machine (augmenting paths)."""
    owner = {}
    def try_job(j, seen):
        for i in edges[j]:
            if i not in seen:
                seen.add(i)
                if i not in owner or try_job(owner[i], seen):
                    owner[i] = j
                    return True
        return False
    for j in frac_jobs:
        if not try_job(j, set()):
            raise RuntimeError("no matching: x is not a vertex")
    return {j: i for i, j in owner.items()}

Worked example: twelve jobs on four machines

Twelve jobs on four machines with integer times drawn from 2 to 19 (NumPy seed 7). Bisection brackets the threshold between 17.66, where LP(T) is infeasible, and 18.06, where it is feasible. Infeasibility is the certificate: no schedule of makespan 17.66 exists, and with integer times OPT is at least 18. The vertex solution at 18.06 splits three jobs; the matching puts them on distinct machines, giving loads of 22, 20, 24 and 7 and a makespan of 24, well inside the 2T bound of 36. An exact mixed-integer solve gives OPT = 21, and a greedy rule that places each job where it finishes earliest gives 23.

So the guarantee algorithm lost to greedy on this instance. That is normal: worst-case guarantees protect you on adversarial inputs, not on average ones. The practical recipe is to run both, take the better, and use the largest infeasible T as the lower bound to report how far from optimal you might be: here at most 23 / 18, about 1.28. A short local search that moves single jobs off the most loaded machine usually closes more of the gap.

Iterative rounding

Iterative rounding pushes the extreme-point idea further. Solve the LP, find a variable the vertex structure guarantees is large, fix it, update the residual LP and re-solve. Jain's 2001 algorithm for survivable network design, choosing edges so each pair of nodes has a required number of disjoint paths, rests on a theorem that every vertex of the residual LP has some edge with x ≥ 1/2. Rounding that edge up and repeating loses at most a factor 2 over all iterations.

The same pattern gives the Shmoys and Tardos result for generalised assignment, cost at most the LP cost with every machine overloaded by at most one job, and the Singh and Lau result for minimum-cost spanning trees with degree bounds, cost at most OPT with every degree bound exceeded by at most one. Iterative rounding costs one LP solve per fixed variable, so warm-start the solver from the previous basis.

Filtering for facility location

Filtering handles problems where a fractional solution spreads each client over many facilities, some of them far away. In uncapacitated facility location you open facilities at a cost and pay each client's distance to its facility. Filtering, due to Lin and Vitter, discards the far assignments for each client, the ones beyond a constant multiple of its fractional distance, and rescales the remainder, losing a constant factor in opening cost. Clients whose filtered neighbourhoods overlap are then clustered, and one facility opened per cluster. Shmoys, Tardos and Aardal gave the first constant factor this way in 1997, about 3.16; later work reached 1.488, against a hardness bound of about 1.463.

Choosing a technique

Failure modes

• Exact float comparisons. 0.4999999999 is not 1/2. Use a tolerance and verify feasibility after rounding.
• Interior-point output. Extreme-point arguments need a vertex; use simplex or crossover.
• Unpruned pairs. Skipping the p > T filter makes the scheduling LP a lie, and the 2T bound fails.
• Weak formulations. A large integrality gap cannot be rounded away; add valid inequalities.
• Reporting only the bound. Report cost / LP* per instance; it is the real certificate.
• LP size. Assignment LPs have m times n variables; prune pairs before building the matrix.

Trade-offs

Rounding gives a guarantee and a per-instance lower bound, at the cost of an LP solve and some code. Greedy is faster and often better on typical data but has no certificate. A MILP solver gives OPT on small instances and a gap estimate on large ones, and is often the right production answer; LP rounding is how you seed it, bound it, or replace it when instances outgrow it.

What to do next

1. Write your problem as an integer program and compute the LP bound on real instances.
2. Find a small instance with a large OPT / LP* ratio to estimate the integrality gap.
3. Try threshold rounding plus pruning first and report cost / LP*.
4. For assignment or scheduling, prune impossible pairs and round from a vertex.
5. Run greedy and a MILP with a time limit alongside; keep the best and the bound.
6. If the gap is too large, strengthen the formulation before tuning the rounding.