Many hard optimisation problems become easy once you let the variables take fractional values. Choosing which warehouses to open, which sets cover a universe, which path each network flow takes, or which variables make the most clauses true are all NP-hard as integer programs. Their linear-programming relaxations, however, solve in polynomial time. The catch is that "open 0.37 of a warehouse" is not a decision. Randomized rounding, introduced by Prabhakar Raghavan and Clark Thompson in 1987, turns the fractional solution into an integer one by treating each fractional value as a probability. Probability theory then shows that the result is close to optimal.
This article builds the technique from first principles. It covers the LP relaxation, the expectation argument, concentration bounds and two classic analyses (set cover and MAX-SAT). It gives runnable Python with SciPy, a worked example with measured results, derandomization, and the engineering lessons that the proofs leave out. It assumes basic LP knowledge; integer programming and LP is the prerequisite.
The idea: probabilities from the LP
Write the problem as an integer program with 0/1 variables x_j, then relax each one to 0 <= x_j <= 1 and solve the LP. Call the optimal fractional solution x* and its cost LP*. For a minimisation problem, the relaxation can only do better than the integer problem, so LP* <= OPT. The rounding step is the simplest one possible: independently for each j, set x_j = 1 with probability x*_j.
Linearity of expectation then gives the cost bound immediately. The expected cost is the sum over j of c_j times P(x_j = 1), which is the sum of c_j x*_j, which equals LP*. In expectation, the rounded solution costs exactly what the LP costs, and that is no more than the optimum. The difficulty is never the objective. It is the constraints: a row that the fractional solution satisfied exactly can be violated after rounding. Every randomized-rounding analysis consists of bounding how often that happens and paying for the fix.
Set cover: the canonical analysis
Weighted set cover makes the pattern concrete. You have a universe of n elements and sets S_j with costs c_j, and you want the cheapest collection of sets that covers every element. The LP is: minimise the sum of c_j x_j subject to, for each element e, the sum of x_j over sets containing e being at least 1, with 0 <= x_j <= 1.
Take an element e covered by sets with fractional values x_1 to x_m, whose sum is at least 1. One round of rounding misses e with probability equal to the product of (1 - x_i). Since 1 - x <= e^-x, that product is at most e^-(sum of x_i), which is at most 1/e. One round therefore leaves each element uncovered with probability at most about 37%, which is too high. So repeat: run t independent rounds and take the union. Element e survives all of them with probability at most e^-t. With t = 2 ln n, each element is uncovered with probability at most 1/n^2. A union bound over n elements makes the probability that anything is uncovered at most 1/n. The expected cost is at most t times LP*, and by Markov's inequality the cost exceeds 4t LP* with probability at most 1/4. So with probability at least 3/4 - 1/n, a single run is a feasible O(log n)-approximation.
This is essentially the best possible. The set cover LP has an integrality gap of order log n, and Dinur and Steurer showed in 2014 that approximating set cover within (1 - epsilon) ln n is NP-hard. Greedy achieves about ln n too. Rounding matters because the same argument extends to problems where greedy has no analysis, such as covering with extra side constraints.
Implementation: round, repair, prune
The code solves the LP with SciPy's HiGHS backend, runs the repeated rounding, repairs any element still uncovered, and then prunes sets made redundant by later picks. The proof needs none of the last two steps. Production use needs both.
import math, random
import numpy as np
from scipy.optimize import linprog
def set_cover_lp(n, sets, costs):
A = np.zeros((n, len(sets)))
for j, S in enumerate(sets):
for e in S:
A[e, j] = 1.0
# sum_j A[e, j] x_j >= 1 is written as -A x <= -1 for linprog
res = linprog(costs, A_ub=-A, b_ub=-np.ones(n), bounds=(0, 1), method="highs")
if res.status != 0:
raise ValueError(res.message)
return res.x, res.fun
def round_cover(x, sets, costs, n, rng, rounds=None):
rounds = rounds or math.ceil(2 * math.log(max(n, 2)))
chosen = set()
for _ in range(rounds):
chosen |= {j for j, xj in enumerate(x) if rng.random() < xj}
covered = set().union(*(sets[j] for j in chosen))
for e in range(n): # repair: cheapest set for each uncovered element
if e not in covered:
j = min((j for j, S in enumerate(sets) if e in S), key=lambda j: costs[j])
chosen.add(j)
covered |= sets[j]
return chosen
def prune(chosen, sets, costs, n):
"""Drop sets, most expensive first, while the rest still covers everything."""
chosen = set(chosen)
for j in sorted(chosen, key=lambda j: -costs[j]):
rest = chosen - {j}
if len(set().union(*(sets[i] for i in rest))) == n:
chosen = rest
return chosen
Worked example: measured costs
A seeded random instance has 60 elements and 40 sets of 4 to 12 elements each, with costs from 1 to 5. The LP optimum is 26.0, and 14 variables are strictly fractional. SciPy's milp finds the true integer optimum, 27. The greedy algorithm finds 30. Each rounding strategy was then run 200 times:
| Strategy | Mean cost | Best |
|---|---|---|
| t = ceil(2 ln 60) = 9 rounds, raw | 44.47 | 35 |
| 9 rounds, then prune | 28.02 | 28 |
| 1 round + repair, raw | 31.52 | 27 |
| 1 round + repair, then prune | 28.80 | 27 |
| Greedy (deterministic) | 30 | 30 |
| Integer optimum (MILP) | 27 | 27 |
The worst-case schedule is wasteful. Nine rounds pick almost every set with positive weight, so the raw cost is 71% above the LP. Repairing what the proof treats as a rare failure, and then pruning, does much better. A single round with repair even found the optimum in some runs. Running the rounding several times and keeping the best is cheap, and the LP value certifies quality as you go: here the best result is within 27/26, about 4%, of a lower bound, without ever knowing OPT.
Sums of rounded variables: Chernoff bounds
Some problems care about a sum of rounded variables, not whether at least one is chosen. The original application was routing with minimum congestion. Each connection picks one path, and the goal is to minimise the maximum load on any edge. Solve the multicommodity-flow LP, decompose each connection's flow into paths, and pick one path per connection with probability equal to its flow. Each edge's load is then a sum of independent 0/1 variables with mean at most the LP congestion C. The Chernoff bound limits how far such a sum can drift above its mean:
P[ X >= (1 + d) * mu ] <= ( e^d / (1 + d)^(1 + d) ) ^ mu for independent 0/1 terms, mean mu
Choose d so the right-hand side is below 1 / (number of edges)^2; union bound over edges.
Result (Raghavan-Thompson): max load = O(log n / log log n) times C, with high probability, when C >= 1.The same pattern of independent choices, a Chernoff bound per constraint and a union bound over constraints handles packing problems, low-congestion embeddings and discrepancy problems. Independence is what makes it work. If your rounding makes correlated choices, for example by picking exactly one path per connection, check that the variables inside each constraint are still independent, or negatively correlated, before you quote the bound.
MAX-SAT and derandomization
Rounding also works for maximisation, and the analysis can depend on the shape of each constraint. In weighted MAX-SAT the LP has a variable y_i per boolean and z_j per clause, with z_j bounded by the sum of the literal values in clause j. Setting each variable true with probability y_i satisfies a clause of length k with probability at least (1 - (1 - 1/k)^k) z_j, which is at least (1 - 1/e) z_j, about 0.632 z_j. That is good for short clauses and weak for long ones. A plain coin flip satisfies a length-k clause with probability 1 - 2^-k, which is weak for short clauses and good for long ones. Goemans and Williamson showed in 1994 that the average of the two guarantees is at least 3/4 z_j for every k. So taking the better of the two assignments is a 3/4-approximation.
Derandomization. For MAX-SAT the expected satisfied weight under independent probabilities can be computed exactly. The method of conditional expectations then removes the randomness. Fix variables one at a time, setting each to whichever value keeps the conditional expectation higher. The expectation never decreases, so the final deterministic assignment is at least as good as the random guarantee:
def expected_weight(p, clauses, w):
"""Clauses are lists of signed 1-based ints; p[i] = P(variable i+1 is True)."""
total = 0.0
for clause, wj in zip(clauses, w):
miss = 1.0
for lit in clause:
miss *= (1 - p[lit - 1]) if lit > 0 else p[-lit - 1]
total += wj * (1 - miss)
return total
def derandomize(p, clauses, w):
p = list(p)
for i in range(len(p)):
p[i] = 1.0; hi = expected_weight(p, clauses, w)
p[i] = 0.0; lo = expected_weight(p, clauses, w)
p[i] = 1.0 if hi >= lo else 0.0
return pOn a seeded instance with 30 variables, 150 clauses of length 1 to 3 and total weight 373, the LP bound is 345.25. The expected weight is 284.38 under coin flips and 321.59 under LP rounding. After derandomization, the coin-flip start reaches 341 and the LP start reaches 338. Both are within about 2% of the LP bound, far above the 3/4 guarantee. The start with the better expectation did not give the better final answer, so run both and keep the better result. For set cover, exact conditional expectations are harder to compute, and pessimistic estimators (Raghavan, 1988) play the same role.
Failure modes
- Infeasible output. A proof that holds with probability 1 - 1/n still fails sometimes. Always verify the constraints and repair, or re-run, instead of trusting the analysis.
- Solver tolerance. LP values such as 0.9999999 or 1e-9 come from floating point. Clip x* to [0, 1] and treat values within a tolerance as exactly 0 or 1. Otherwise variables the LP set to 1 will occasionally be dropped.
- Large integrality gap. If the LP is weak, the rounding inherits the weakness. Add valid inequalities, or use a stronger formulation, before tuning the rounding.
- Hidden dependence. Chernoff bounds need independence, or at least negative correlation. Rounding that couples variables can break the analysis without any visible error.
- LP solve dominates runtime. On large instances, solving the LP costs far more than rounding it. Reuse warm starts, solve the dual, or use column generation.
- Unreproducible results. Seed the random generator and log the seed with each solution, so that an unusual plan can be reproduced and investigated.
Trade-offs
Randomized rounding versus deterministic threshold rounding. For vertex cover, rounding every x* >= 1/2 up gives a clean factor of 2 with no randomness. Threshold rounding fits when each constraint has few variables, and randomized rounding fits when constraints are long sums. Rounding versus exact MILP. For mid-sized instances a MILP solver finds the optimum, and rounding is better used as its primal heuristic. Rounding comes into its own at scale, or when you need a provable guarantee and a fixed time budget. Rounding versus greedy. Greedy is simpler and often competitive, but the LP gives a lower bound that greedy never gives you, so you know how far from optimal you might be. LP duality explains why that bound is trustworthy.
What to do next
- Write your problem as a 0/1 integer program, and check how weak the LP relaxation is on small instances by comparing LP* with MILP's OPT.
- Implement the round, repair, prune pipeline, and report each solution's cost as a ratio to LP*.
- Run the rounding several times with logged seeds and keep the best. Compare against greedy.
- If you need determinism, derandomize with conditional expectations where they can be computed exactly.
- For sum-type constraints, check independence before relying on Chernoff bounds.
- Read approximation algorithms for the wider toolbox.