Set cover: the fewest plants that reach every metro
Corrugated board is mostly air, so a box plant sells within about 250 miles of itself and not much further. Take the 69 largest US metros. A plant in any one of them reaches every metro inside that radius. What is the fewest plants that reach all 69? That question is set cover, the opening problem of Korte & Vygen's chapter on approximation algorithms, and it is NP-hard, which means the honest answers are a fast one that is provably not too bad and a slow one that is exact. Here are both, on the map.
The instance
The universe is the 69 metros. There are 69 candidate sets, one per metro: a plant at Detroit covers the eleven metros within 250 miles of Detroit, a plant at Boise covers two. The largest set has k = 11 elements. Minimize the number of sets whose union is everything.
The greedy algorithm
Pick the set that covers the most still-uncovered elements. Remove what it covers. Repeat until nothing is left. That is the whole algorithm, and the first pick is already a surprise: not New York, not Los Angeles, but Detroit, because the 250-mile disc around Detroit holds Chicago, Pittsburgh, Cincinnati, Columbus, Indianapolis, Cleveland, Buffalo, Grand Rapids, Toledo and Dayton. The densest circle of cities in the country is the Great Lakes, and the algorithm has no prior about which cities matter, only counts.
Watch the gains fall. The first five picks cover 40 metros. The last seven picks cover one metro each: New York for Boston, Houston for San Antonio, Washington for Raleigh, then Denver, New Orleans, Salt Lake City and Boise for themselves. This is the shape of every greedy set cover: big early wins, then a long tail of singletons, and it is exactly why the approximation ratio is logarithmic. Each pick covers at least a 1/OPT fraction of what remains (the optimum covers all of it with OPT sets, so some set covers at least that share), and shrinking the remainder by that factor each time takes about OPT · ln n picks to reach zero.
What the book proves, and what happened
Korte and Vygen give Chvátal's bound: greedy uses at most H(k) · OPT sets, where H(k) = 1 + ½ + ⅓ + … + 1/k and k is the largest set. Here H(11) ≈ 3.02, and the optimum turns out to be 18, so the guarantee is 54. Greedy used 20. The guarantee is a ceiling, not a forecast; on a geographic instance where sets are discs that overlap smoothly, greedy is within two of exact. The gap is real, though. The optimum is a different set of cities: it takes Riverside instead of Los Angeles, Orlando instead of Tampa, Indianapolis, Nashville and Raleigh instead of Knoxville, Philadelphia and Washington. Ten cities appear in both answers; eight do not.
Two more numbers came along for free. The LP relaxation, which lets you build 0.4 of a plant, gives a lower bound of exactly 18, so once CP-SAT found a cover of 18 it was done without having to prove anything further; the LP had already proved it. And the LP solution happened to land on an integer vertex (a different optimal cover from CP-SAT's, since there are several). That is luck, not a theorem; on the instance below it would not happen.
The instance that breaks greedy
The bound is tight, and the book's construction for showing it is worth seeing because it is so small. Two rows of 15 elements are the optimum: two sets cover everything. Four columns of sizes 2, 4, 8 and 16, each split evenly across the rows, tile the same 30 elements. Greedy counts: the biggest column has 16, each row has 15, so it takes the column. Now each row has 7 left and the next column has 8. It takes the column. And so on, four picks for an optimum of two, and with k columns the ratio is k/2 ≈ log₂ n. Greedy never looks back, and the rows are only ever the second-best choice.
Adding to a network that already exists
No one builds from scratch. The fictional producer in the network model has 14 plants, and they leave 16 metros outside any 250-mile disc: all of Florida, New England past New York, Northern California, New Orleans, Salt Lake, Albuquerque, El Paso, Boise, Little Rock, Reno, Buffalo. Fewest new plants to reach them is set cover again, on a smaller universe. Greedy and the optimum both say 9, though they disagree about where; switch the map to add to the 14 to see both.
And there is the number that an operations VP would notice. Fourteen existing plants plus nine new ones is 23 plants to cover a country that 18 can cover. Five plants, roughly a quarter of the network, is the price of having built it one acquisition at a time, each one locally sensible. That is the greedy algorithm again, running over decades instead of milliseconds, with the same ratio.
The Python
From setcover.py, which also builds the sets, runs the LP and CP-SAT, checks Chvátal's bound on every run, and verifies the bad instance for k = 2…5. The JavaScript on this page runs the same greedy on the same sets; the exact answers are the solver's, embedded.
def greedy(sets, universe):
"""Pick the set covering the most still-uncovered elements; repeat."""
uncovered, chosen, trace = set(universe), [], []
while uncovered:
best = max(range(len(sets)), key=lambda i: len(sets[i] & uncovered))
gain = len(sets[best] & uncovered)
if gain == 0:
raise ValueError("some element is in no set")
chosen.append(best)
uncovered -= sets[best]
trace.append((best, gain, len(uncovered)))
return chosen, trace
def bad_instance(k):
"""Two rows of 2^k - 1 elements (the optimum), tiled by columns of size 2, 4, ..., 2^k."""
half = 2**k - 1
rows = [frozenset(range(half)), frozenset(range(half, 2*half))]
cols, start = [], 0
for size in (2**j for j in range(1, k + 1)):
h = size // 2
cols.append(frozenset([*range(start, start + h), *range(half + start, half + start + h)]))
start += h
return list(range(2*half)), rows + cols
# the checks that run every time
assert math.ceil(lp_bound) <= len(optimum) <= len(greedy_cover)
assert len(greedy_cover) <= harmonic(k) * len(optimum)
for k in (2, 3, 4, 5):
U, S = bad_instance(k)
assert len(exact(S, U)) == 2 and len(greedy(S, U)[0]) == k