Cutting planes: how a solver rounds without rounding
Every scheduling, cubing and network model I've written ends in an integer program, and every one of them is solved by a machine that starts from the wrong answer: the LP relaxation, where you can commit 7.57 truckloads. Chapter 14 of Papadimitriou & Steiglitz is Gomory's 1958 answer to the question of how to get from there to an integer without guessing. You read a row off the simplex tableau, write down an inequality every integer point must satisfy and the current fraction doesn't, add it, re-solve, repeat. Here it is on an instance small enough to draw, with every row shown.
The week
A plant has 19 corrugator hours and 24 converting hours uncommitted next week and two customers who will take as many whole truckloads as it can promise. A load of printed RSC takes 2 corrugator hours and 3 on the folder-gluer and contributes $900. A load of double-wall die-cut takes 3 corrugator hours (double wall runs slow) and 1 on the die cutter, and contributes $1,100. In hundreds of dollars:
maximize 9·x₁ + 11·x₂ x₁ = RSC loads, x₂ = double-wall loads
subject to 2·x₁ + 3·x₂ ≤ 19 corrugator hours
3·x₁ + 1·x₂ ≤ 24 converting hours
x₁, x₂ ≥ 0 and integer half a truck doesn't shipDrop the word integer and the simplex method says: promise 7.57 RSC loads and 1.29 double-wall loads for $8,229. Round it to (8, 1) and you've oversold the converting floor by an hour. Round down to (7, 1) and you make $7,400. The true integer optimum is (5, 3), $7,800: fewer of the product the LP loved, more of the other, and no rounding of the LP point gets you anywhere near it.
Watch the cuts
The shaded region is the LP's feasible set; the dots inside it are the plans that actually ship. The dashed outline is the integer hull, the smallest convex set containing those dots. If you could hand the simplex method that hull instead of the LP polygon, its answer would be integer every time, because every corner of the hull is a dot. Cutting planes approach the hull from outside, one inequality at a time, and they only have to get it right near the optimum.
Where a cut comes from
At the LP optimum the tableau has a row for each basic variable. The row for x₁ reads
x₁ − (1/7)·s₁ + (3/7)·s₂ = 53/7where s₁ and s₂ are the slacks of the two constraints (unused corrugator and converting hours). Now the chapter's one idea. In any integer solution, x₁, s₁ and s₂ are all non-negative integers (slack hours are whole if everything else is). Split every coefficient into its integer part and its fractional part, where the fractional part is always in [0, 1): −1/7 = −1 + 6/7, 3/7 = 0 + 3/7, 53/7 = 7 + 4/7. Move the integer parts to one side:
(6/7)·s₁ + (3/7)·s₂ − 4/7 = 7 − x₁ − (−1)·s₁ − 0·s₂ = an integerThe left side is an integer, and because s₁, s₂ ≥ 0 it is at least −4/7, so it is at least 0:
(6/7)·s₁ + (3/7)·s₂ ≥ 4/7 the Gomory cutEvery feasible integer point satisfies this. The current LP optimum doesn't: s₁ = s₂ = 0 there, and 0 ≥ 4/7 is false. Substitute s₁ = 19 − 2x₁ − 3x₂ and s₂ = 24 − 3x₁ − x₂ and multiply through by 7, and the cut in the original variables is 3x₁ + 3x₂ ≤ 26, which is the first red line in the figure. Add it with its own slack t₁, re-optimize (the dual simplex does it in one pivot, because the old basis is still dual-feasible), and the new LP optimum is (7, 5/3). Still fractional. The row for x₂ now gives (1/3)·t₁ ≥ 2/3, which in the original variables is x₁ + x₂ ≤ 8, and the next LP optimum is (5, 3). Integer. Done, and provably optimal, because it is the optimum of a relaxation that contains every integer point.
What the chapter proves, and what practice found
Gomory proved the method terminates: with a lexicographic rule for choosing pivots and source rows, a finite number of cuts reaches an integer optimum. The proof is in the chapter and it is the kind that tells you nothing about how many. On this instance, two. On a neighboring instance I tried first (9 and 14 for the objective, 32 and 29 hours), six cuts to move the optimum from (3.5, 5) to (3, 5), each one shaving a sliver. Cuts creep. For forty years that was the verdict on pure cutting planes: elegant, finite, and slow, while branch and bound did the real work.
What changed in the 1990s was not the cut but where you put it. Add a few rounds of Gomory's mixed-integer cuts at each node of a branch-and-bound tree and the LP bounds tighten enough to prune most of the tree. That hybrid, branch and cut, is what sits inside CBC, SCIP, Gurobi and CPLEX, and it is why the truck-cubing MILP on sixty units solves in seconds instead of hours. When I say a solver "found the optimum," this is what it did: relaxed, cut, branched, cut again, and proved the bound met the incumbent.
The Python
gomory.py is the whole method in exact rational arithmetic (Python's Fraction), so the tableau above is printed from the code, not typed. The test solves 290 random small instances, checks the answer against brute force, and checks that no cut ever removed a feasible integer point.
def add_gomory_cut(self):
# source row: a basic variable with a fractional value
r = next(i for i in range(len(self.rows)) if frac(self.rhs[i]) != 0)
coeffs = [frac(a) if j not in self.basis else F(0) for j, a in enumerate(self.rows[r])]
f0 = frac(self.rhs[r])
# cut: sum_j frac(a_rj) x_j >= f0 -> with a new slack t >= 0:
# -sum_j frac(a_rj) x_j + t = -f0 (infeasible for the current point; dual simplex fixes it)
for row in self.rows:
row.append(F(0))
self.rows.append([-v for v in coeffs] + [F(1)])
self.rhs.append(-f0)
self.z.append(F(0))
self.basis.append(len(self.z) - 1)
def gomory(c, A, b):
T = Tableau(c, A, b)
T.primal_simplex()
while not T.is_integer():
T.add_gomory_cut()
T.dual_simplex()
return T