Presolve reductions

Deep dive

How it works

Parts:

• row activity bounds — minact/maxact per row, updated as bounds change
• reduction rules — redundant/forcing rows, bound tightening, dual fixing, coefficient tightening, aggregation, probing
• work queue — rows and columns touched by the last change, revisited until a fixpoint
• postsolve stack — append-only undo records that map the reduced solution back

Operations (steps below): compute activities, apply rules, repeat to a fixpoint, solve, undo in reverse.

model ─► [activities ─► rules ─► changed rows] ─repeat─► reduced model ─► solver
                                      │ push undo records                    │
                                      └──────────► postsolve stack ◄── solution (reverse scan)

Model: min c·x, rows lo ≤ a·x ≤ hi, bounds lb ≤ x ≤ ub, some x_j integer. Row minact = Σ a_j·(a_j>0 ? lb_j : ub_j); maxact mirrors it. Infinite terms are counted, not summed.

1. Infeasible / redundant row: minact > hi or maxact < lo proves infeasibility; lo ≤ minact and maxact ≤ hi makes the row redundant, delete it. minact == hi is a forcing row: fix every variable at its minimising bound.
2. Bound tightening: for a_j > 0, finite hi: x_j ≤ (hi − (minact − a_j·lb_j)) / a_j; three symmetric cases. Round for integers. A singleton row has empty residual and just becomes a bound.
3. Fixed variable (lb == ub): substitute out, shift row sides, push on the postsolve stack.
4. Dual fixing: a down-lock is a row that decreasing x_j may violate (a>0 with finite lo, a<0 with finite hi). If c_j ≥ 0 and x_j has none, set x_j = lb_j; mirror with up-locks.
5. Coefficient tightening (Savelsbergh): in Σ a x ≤ b, binary x_j, a_j > 0: if maxact − a_j < b, let d = b − (maxact − a_j), then a_j -= d; b -= d. Same integer points, stronger LP: 5x + y ≤ 4 becomes 2x + y ≤ 1, cutting off x = 0.8, y = 0.
6. Aggregation: from a doubleton equation a x + b y = c substitute x = (c − b y)/a; watch fill-in.
7. Duplicates and dominance: parallel rows (hash the normalised row, keep tighter sides), parallel columns (merge), dominated columns (fix at a bound).
8. Probing: fix a binary to 0 then 1 and propagate. One branch infeasible ⇒ fix the other; same fixing in both ⇒ fix it; implications feed clique tables.
9. Drive 1-8 from a worklist of changed rows/columns until a fixpoint or round limit.

Postsolve: every reduction pushes a record (fixed value, x = (c − b y)/a, …); after solving, pop in reverse to rebuild the full x.

Code

Steps 1-4 with integer rounding; postsolve here is just fixed[j].

const INF: f64 = f64::INFINITY;

struct Row { a: Vec<(usize, f64)>, lo: f64, hi: f64 } // lo <= a.x <= hi
struct Prob { lb: Vec<f64>, ub: Vec<f64>, int: Vec<bool>, c: Vec<f64>, rows: Vec<Row> } // min c.x

#[derive(Debug, PartialEq)]
enum Status { Infeasible, Reduced(Vec<Option<f64>>) } // Reduced(fixed): fixed[j] is the postsolve value

/// (finite part, number of infinite terms) of the min or max activity of a row.
fn activity(r: &Row, p: &Prob, max: bool) -> (f64, usize) {
    let (mut s, mut inf) = (0.0, 0);
    for &(j, a) in &r.a {
        let b = if (a > 0.0) == max { p.ub[j] } else { p.lb[j] };
        if b.is_infinite() { inf += 1 } else { s += a * b }
    }
    (s, inf)
}

fn presolve(p: &mut Prob) -> Status {
    let n = p.lb.len();
    let mut fixed = vec![None; n];
    loop {
        let mut changed = false;
        // 1. fixed variables: substitute out, shift row sides
        for j in 0..n {
            if fixed[j].is_none() && p.ub[j] - p.lb[j] <= 1e-9 { fixed[j] = Some(p.lb[j]); changed = true; }
        }
        for r in p.rows.iter_mut() {
            let mut shift = 0.0;
            r.a.retain(|&(j, a)| match fixed[j] { Some(v) => { shift += a * v; false } None => true });
            r.lo -= shift; r.hi -= shift;
        }
        // 2. per row: infeasible, redundant, bound tightening (a singleton row just becomes a bound)
        for k in 0..p.rows.len() {
            if p.rows[k].a.is_empty() { continue; }
            let (mn, mni) = activity(&p.rows[k], p, false);
            let (mx, mxi) = activity(&p.rows[k], p, true);
            let (lo, hi) = (p.rows[k].lo, p.rows[k].hi);
            if (mni == 0 && mn > hi + 1e-7) || (mxi == 0 && mx < lo - 1e-7) { return Status::Infeasible; }
            if mni + mxi == 0 && mn >= lo - 1e-9 && mx <= hi + 1e-9 { p.rows[k].a.clear(); changed = true; continue; }
            for idx in 0..p.rows[k].a.len() {
                let (j, a) = p.rows[k].a[idx];
                // activity of the row without term j, if bounded
                let res = |max: bool, tot: (f64, usize)| {
                    let b = if (a > 0.0) == max { p.ub[j] } else { p.lb[j] };
                    if b.is_infinite() { (tot.1 == 1).then_some(tot.0) } else { (tot.1 == 0).then(|| tot.0 - a * b) }
                };
                let (mut nl, mut nu) = (-INF, INF);
                if let (true, Some(r)) = (hi.is_finite(), res(false, (mn, mni))) { let v = (hi - r) / a; if a > 0.0 { nu = v } else { nl = v } }
                if let (true, Some(r)) = (lo.is_finite(), res(true, (mx, mxi))) { let v = (lo - r) / a; if a > 0.0 { nl = nl.max(v) } else { nu = nu.min(v) } }
                if p.int[j] { nl = (nl - 1e-6).ceil(); nu = (nu + 1e-6).floor(); }
                // continuous bounds only move if the gain is real, else we chase 1e-12 steps
                let big = |old: f64, new: f64| p.int[j] || (old - new).abs() > 1e-3 * (1.0 + new.abs());
                if nl > p.lb[j] && big(p.lb[j], nl) { p.lb[j] = nl; changed = true; }
                if nu < p.ub[j] && big(p.ub[j], nu) { p.ub[j] = nu; changed = true; }
                if p.lb[j] > p.ub[j] + 1e-7 { return Status::Infeasible; }
            }
        }
        // 3. dual fixing via locks: no row blocks the objective-improving direction
        let (mut dn, mut up) = (vec![0; n], vec![0; n]);
        for r in &p.rows { for &(j, a) in &r.a {
            if (a > 0.0 && r.lo.is_finite()) || (a < 0.0 && r.hi.is_finite()) { dn[j] += 1 } // decreasing x_j may violate
            if (a > 0.0 && r.hi.is_finite()) || (a < 0.0 && r.lo.is_finite()) { up[j] += 1 }
        } }
        for j in 0..n {
            if fixed[j].is_some() { continue; }
            if p.c[j] >= 0.0 && dn[j] == 0 && p.lb[j].is_finite() { p.ub[j] = p.lb[j]; changed = true; }
            else if p.c[j] <= 0.0 && up[j] == 0 && p.ub[j].is_finite() { p.lb[j] = p.ub[j]; changed = true; }
        }
        if !changed { break; }
    }
    p.rows.retain(|r| !r.a.is_empty());
    Status::Reduced(fixed)
}

fn main() {
    let row = |a: &[(usize, f64)], lo, hi| Row { a: a.to_vec(), lo, hi };
    // min -x0 - x1 + x3; x0,x3 integer, x1 binary, x2 continuous
    let mut p = Prob { lb: vec![0.0; 4], ub: vec![10.0, 1.0, INF, 5.0], int: vec![true, true, false, true], c: vec![-1.0, -1.0, 0.0, 1.0],
        rows: vec![row(&[(0, 2.0)], -INF, 7.0),            // singleton: x0 <= 3.5, rounded to 3
                   row(&[(0, 1.0), (1, 1.0)], -INF, 10.0), // redundant once x0 <= 3
                   row(&[(1, 1.0), (2, 1.0)], 1.5, INF),
                   row(&[(2, 1.0), (3, 1.0)], -INF, 4.0)] };
    // fixings solve it outright: x0 = 3, x1 = 1, x3 = 0, x2 = 0.5
    assert_eq!(presolve(&mut p), Status::Reduced(vec![Some(3.0), Some(1.0), Some(0.5), Some(0.0)]));
    assert!(p.rows.is_empty());
}

Compiled and run with the asserts shown. Real codes add the record stack, a worklist instead of full sweeps and tolerance control.

Complexity and constants

• One sweep is O(nnz): activities, then a residual per term. Rounds to fixpoint are few; codes cap rounds and stop when a round removes little.
• Memory: row-wise and column-wise copies of the matrix (order of 12 bytes per nonzero each: u32 index + f64) plus the reduction stack.
• Cost drivers: nonzeros touched, one hash per row/column for duplicates, and probing (one propagation per binary and value; budget it).
• The payoff is downstream: every simplex iteration, LU refactorisation and branch-and-bound node runs on the smaller model.

Efficiency

• Bound: memory bandwidth over the matrix; almost no floating-point work.
• Work avoided: each removed row or column shrinks every later LP and node. The Gurobi presolve paper reports presolve is often the difference between an intractable and a solvable model.
• Worklist, not sweeps: re-examine only rows and columns touched by the last change.
• Operation-count budgets for probing, dominance and sparsification (not wall-clock), so runs stay cheap and deterministic.
• Parallel: PaPILO runs independent reductions concurrently and applies them in a fixed order, so results stay reproducible.
• Cheaper alternative: for tiny models, or repeated solves where only right-hand sides change, keep the reduced model and rerun only bound propagation.

Implementation notes

• Keep the original model untouched; the stack is an append-only Vec of small enum records, so postsolve is one reverse scan.
• Use tolerances to detect redundancy (order of 1e-6) but apply exact arithmetic; for continuous variables skip tightenings under ~1e-3 relative, or rounds never end.
• Round integer bounds with slack, ceil(v − 1e-6) / floor(v + 1e-6), so float noise cannot cut off valid integers.
• Count infinite contributions instead of adding inf; inf − inf never arises.
• Reject aggregations with small pivots relative to the row and cap fill-in.
• Re-run after restarts and root cuts: new fixings are cheap to find because most of the model is unchanged.

Pitfalls

• Dual reductions break solution counting, enumeration and original-row duals.
• Basis and dual postsolve (LP) needs extra records per reduction (Andersen & Andersen); primal-only MIP postsolve does not.
• Tolerance chasing and mutually undoing rules loop forever; cap rounds.
• Big-M: presolve tightens M only as far as the bounds allow; loose variable bounds leave loose M. Give tight bounds.
• A buggy presolve yields "infeasible" with no row to blame; keep an off switch for debugging.
• Integrality-based reductions are invalid after relaxing variables to continuous.

Variants and alternatives

In production

Gurobi, CPLEX, Xpress, SCIP, HiGHS and Cbc presolve by default; the reduced LP then goes to simplex or an interior-point method. Pricing masters in column generation are presolved too. It is on by default in SciPy linprog (presolve=True, docs) and so runs for every field that hands generated models to a solver: logistics, energy dispatch, scheduling and packing; the MIPLIB 2017 benchmark gathers instances from many application areas (Gleixner et al.).

Implementations

Sources

• Erling D. Andersen, Knud D. Andersen, Presolving in linear programming, Mathematical Programming 71, 1995. paper
• Martin W. P. Savelsbergh, Preprocessing and probing techniques for mixed integer programming problems, ORSA Journal on Computing 6(4), 1994. paper
• Tobias Achterberg, Robert E. Bixby, Zonghao Gu, Edward Rothberg, Dieter Weninger, Presolve reductions in mixed integer programming, INFORMS Journal on Computing 32(2), 2020. paper
• Gerald Gamrath, Thorsten Koch, Alexander Martin, Matthias Miltenberger, Dieter Weninger, Progress in presolving for mixed integer programming, Mathematical Programming Computation 7, 2015. paper
• Ambros Gleixner et al., MIPLIB 2017: data-driven compilation of the 6th mixed-integer programming library, Math. Program. Comput. 13, 2021. paper
• Ambros Gleixner, Leona Gottwald, Alexander Hoen, PaPILO: A parallel presolving library for integer and linear optimization with multiprecision support, INFORMS Journal on Computing, 2023. paper · code