Cutting planes

Deep dive

How it works

Parts:

• LP relaxation — current optimal basis x*, re-solved by dual simplex after each round
• separators — one routine per family (GMI, MIR, cover, flow cover, clique, zero-half) that returns violated inequalities
• cut pool and filter — keeps the cuts with best efficacy, parallelism and numerics
• branch-and-cut driver — decides: another round, or branch

Operations (steps below): solve, separate, filter, add, re-solve; branch when the bound stalls.

LP solve ─► x* integral? ─yes─► done
   ▲            │no
   │            ▼
 add best ◄─ filter ◄─ separate (violated cuts?) ─none/stalled─► branch

1. Solve the relaxation; take optimal x* and its basis. If integral, stop.
2. Separate: find an inequality a·x <= b valid for all integer points with a·x* > b.
3. Gomory fractional cut (pure integer, all variables and slacks integer): a tableau row x_i + sum_{j nonbasic} t_j x_j = r with f0 = frac(r) > 0 gives sum_j frac(t_j) x_j >= f0. Why: sum_j frac(t_j) x_j = r - x_i - sum_j floor(t_j) x_j is f0 plus an integer and is >= 0, hence >= f0.
4. GMI (mixed): with f_j = frac(t_j), nonbasic variables at 0: sum_{int, f_j<=f0} f_j/f0 x_j + sum_{int, f_j>f0} (1-f_j)/(1-f0) x_j + sum_{cont, t_j>0} t_j/f0 x_j + sum_{cont, t_j<0} -t_j/(1-f0) x_j >= 1. Dominates the fractional cut.
5. Cover for sum a_j x_j <= b, x binary: a set C with sum_C a_j > b gives sum_C x_j <= |C|-1. Separation is a small knapsack: minimize sum (1-x*_j) z_j s.t. sum a_j z_j > b; a cover cut is violated iff the minimum is < 1.
6. MIR for sum a_j x_j - s <= b, x integer, s >= 0, f = frac(b) > 0: sum (floor(a_j) + max(0, f_j - f)/(1-f)) x_j - s/(1-f) <= floor(b). Solvers aggregate rows and complement variables to find violated ones.
7. Filter (efficacy, parallelism, density, numerical range), add the few best, re-solve with the dual simplex, loop. Stop after a round limit or when the bound stalls (tailing off), then branch.

Worked example (Wolsey): max 4x1 - x2 s.t. 7x1 - 2x2 <= 14, x2 <= 3, 2x1 - 2x2 <= 3. LP optimum (20/7, 3), value 8.43. The x1 row yields s1 + 2 s2 >= 6, i.e. x1 <= 2; re-solve gives (2, 0.5); the next cut is x1 - x2 <= 1; re-solve gives (2, 1), value 7, the integer optimum.

Code

Dense-tableau Gomory fractional cuts plus dual-simplex re-solve; for small integer-data problems.

const EPS: f64 = 1e-9;

/// Dense tableau for  min d·x  s.t. T x = rhs, x >= 0. Row i has basic column basis[i];
/// each row is [coeffs for all columns | rhs]. Initially: max c·x, Ax <= b, b >= 0, integer data.
pub struct Lp { t: Vec<Vec<f64>>, d: Vec<f64>, z: f64, basis: Vec<usize>, n_orig: usize }

impl Lp {
    pub fn new(a: &[Vec<f64>], b: &[f64], c: &[f64]) -> Lp {
        let (m, n) = (a.len(), c.len());
        let t = (0..m).map(|i| {
            let mut r = a[i].clone();
            r.extend((0..m).map(|k| (k == i) as u8 as f64)); // slack i
            r.push(b[i]);
            r
        }).collect();
        let mut d: Vec<f64> = c.iter().map(|v| -v).collect();
        d.resize(n + m, 0.0);
        Lp { t, d, z: 0.0, basis: (n..n + m).collect(), n_orig: n }
    }

    fn pivot(&mut self, r: usize, k: usize) {
        let p = self.t[r][k];
        for v in self.t[r].iter_mut() { *v /= p; }
        let pr = self.t[r].clone();
        for i in 0..self.t.len() {
            let f = self.t[i][k];
            if i != r && f != 0.0 { for (v, w) in self.t[i].iter_mut().zip(&pr) { *v -= f * w; } }
        }
        let f = self.d[k];
        for (v, w) in self.d.iter_mut().zip(&pr) { *v -= f * w; }
        self.z -= f * pr[pr.len() - 1];
        self.basis[r] = k;
    }

    /// Primal simplex, Bland's rule (slow but cycle-free). Needs rhs >= 0.
    fn primal(&mut self) {
        while let Some(k) = (0..self.d.len()).find(|&j| self.d[j] < -EPS) {
            let w = self.t[0].len() - 1;
            let r = (0..self.t.len()).filter(|&i| self.t[i][k] > EPS)
                .min_by(|&i, &j| (self.t[i][w] / self.t[i][k]).partial_cmp(&(self.t[j][w] / self.t[j][k])).unwrap())
                .expect("unbounded");
            self.pivot(r, k);
        }
    }

    /// Dual simplex: restores primal feasibility after a cut row with negative rhs was added.
    fn dual(&mut self) -> bool {
        loop {
            let w = self.t[0].len() - 1;
            let Some(r) = (0..self.t.len()).filter(|&i| self.t[i][w] < -EPS)
                .min_by(|&i, &j| self.t[i][w].partial_cmp(&self.t[j][w]).unwrap()) else { return true };
            let Some(k) = (0..w).filter(|&j| self.t[r][j] < -EPS)
                .min_by(|&i, &j| (self.d[i] / -self.t[r][i]).partial_cmp(&(self.d[j] / -self.t[r][j])).unwrap())
                else { return false }; // row proves infeasibility
            self.pivot(r, k);
        }
    }

    /// Gomory fractional cut from row r:  sum_j frac(t_rj) x_j >= frac(rhs)  over nonbasic j.
    /// Valid only if all variables, slacks included, are integer. Added as -sum f_j x_j + s = -f0.
    fn add_gomory_cut(&mut self, r: usize) {
        let w = self.t[0].len() - 1;
        let frac = |v: f64| { let r = v.round(); if (v - r).abs() < 1e-9 { 0.0 } else { v - v.floor() } };
        let mut row: Vec<f64> = (0..w).map(|j| if self.basis.contains(&j) { 0.0 } else { -frac(self.t[r][j]) }).collect();
        row.push(1.0);
        row.push(-frac(self.t[r][w]));
        for i in 0..self.t.len() { self.t[i].insert(w, 0.0); }
        self.d.push(0.0);
        self.basis.push(w);
        self.t.push(row);
    }

    /// Cutting-plane loop. Returns (objective, x, rounds).
    pub fn solve_integer(&mut self, max_rounds: usize) -> Option<(f64, Vec<f64>, usize)> {
        self.primal();
        for round in 0..=max_rounds {
            let w = self.t[0].len() - 1;
            // most fractional basic variable among the original ones (ties: first)
            let best = (0..self.t.len()).filter(|&i| self.basis[i] < self.n_orig)
                .map(|i| (i, (self.t[i][w] - self.t[i][w].round()).abs()))
                .filter(|&(_, f)| f > 1e-6).max_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
            let Some((r, _)) = best else {
                let mut x = vec![0.0; self.n_orig];
                for (i, &j) in self.basis.iter().enumerate() { if j < self.n_orig { x[j] = self.t[i][w].round(); } }
                return Some((self.z, x, round));
            };
            if round == max_rounds { return None; }
            self.add_gomory_cut(r);
            if !self.dual() { return None; }
        }
        None
    }
}

fn main() {
    // Wolsey: max 4x1 - x2  s.t. 7x1 - 2x2 <= 14, x2 <= 3, 2x1 - 2x2 <= 3, x integer >= 0
    let a = vec![vec![7.0, -2.0], vec![0.0, 1.0], vec![2.0, -2.0]];
    let mut lp = Lp::new(&a, &[14.0, 3.0, 3.0], &[4.0, -1.0]);
    let (z, x, rounds) = lp.solve_integer(50).unwrap();
    assert!((z - 7.0).abs() < 1e-9 && x == [2.0, 1.0] && rounds == 2);
}

Compiled and run: LP bound 8.4286, two cuts, z = 7, x = (2, 1); also matched brute force on 300 random small instances.

Complexity and constants

• One GMI cut needs one row of B^-1 A: a sparse solve with the basis factors plus a pass over the row, order of O(nnz); the dense code above pays O(m n) per pivot.
• Each round costs one warm-started dual-simplex re-solve (typically few pivots), far cheaper than a cold start.
• Each cut adds a row, so the basis grows, density grows, and later LPs slow down; cut pools exist to limit that.
• A cut derived from rows that contain earlier cuts has higher rank; high rank means large denominators and unstable numbers.

Efficiency

Bound: LP re-solve time, and LP size after cuts, not separation itself.

• Few, strong cuts per round beat many weak ones: SCIP's Gomory separator defaults to at most 200 cuts per root round and 10 root rounds (5 rounds, 50 cuts per node), and only tries rows whose basic variable is at least 0.01 from integral.
• Root focus: separating at every node multiplies LP cost by the tree size; SCIP runs its Gomory separator at depth frequency 10.
• Drop cuts that stay non-binding so the LP does not keep growing (SCIP: dynamiccuts).
• Sparsity filter: reject near-parallel and dense cuts; each non-zero is paid in every later LP.
• Cheaper: only structural cuts (cover, clique) when the model has them.

Implementation notes

• Compute the tableau row from the factored basis (e_i^T B^-1 A), never form B^-1. Skip rows whose r is within ~1e-6 of an integer.
• Snap tiny fractional parts to 0 before computing frac, as frac(0.99999999) creates a bogus coefficient near 1. Reject cuts whose largest/smallest coefficient ratio is huge (order of 1e6 and above is risky).
• A cut derived using branching bounds is only valid in that subtree (local); never put it in the global pool.

Pitfalls

• Floating-point noise producing an invalid cut that removes the optimum: the solver reports a worse "optimal" value with no error.
• Tailing off: later rounds improve the bound by tiny amounts; stop and branch.
• Pure-integer formula on rows with continuous variables or non-integer slacks: invalid; use GMI.
• Adding every violated cut: the LP gets denser and slower than branching would have been.

Variants and alternatives

In production

Gurobi, CPLEX, Xpress, SCIP, HiGHS, CBC (via Cgl), and the LP cuts of OR-Tools CP-SAT. Branch-and-cut: Crowder-Johnson-Padberg (0-1 programs), Padberg-Rinaldi (TSP). The Concorde TSP solver (branch-and-cut) is applied to logistics routing (parcel pickup trucks, warehouse stacker cranes), manufacturing (drilling holes in circuit boards) and VLSI chip-design instances (applications, VLSI data sets).

Implementations

Sources

• Ralph E. Gomory, Outline of an algorithm for integer solutions to linear programs, Bull. AMS 64(5), 1958. doi
• Harlan Crowder, Ellis L. Johnson, Manfred Padberg, Solving large-scale zero-one linear programming problems, Oper. Res. 31(5), 1983. doi
• George L. Nemhauser, Laurence A. Wolsey, A recursive procedure to generate all cuts for 0-1 mixed integer programs, Math. Program. 46, 1990. doi
• Manfred Padberg, Giovanni Rinaldi, A branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems, SIAM Review 33(1), 1991. doi
• Egon Balas, Sebastian Ceria, Gerard Cornuejols, N. Natraj, Gomory cuts revisited, Oper. Res. Lett. 19, 1996. doi
• Hugues Marchand, Laurence A. Wolsey, Aggregation and mixed integer rounding to solve MIPs, Oper. Res. 49(3), 2001. doi
• Gerard Cornuejols, Valid inequalities for mixed integer linear programs, Math. Program. 112, 2008. doi
• Giacomo Zanette, Matteo Fischetti, Egon Balas, Lexicography and degeneracy: can a pure cutting plane algorithm work?, Math. Program. 130, 2011. doi
• Michele Conforti, Gerard Cornuejols, Giacomo Zambelli, Integer Programming, Springer 2014. doi
• SCIP source, src/scip/sepa_gomory.c (parameter defaults quoted above). code