Simplex method

Deep dive

How it works

In: c, A, b (and optionally a starting basis) · Out: optimal vertex x, duals, final basis, or an infeasible/unbounded certificate · State: the basis (and its factorization), which makes warm starts possible · Cost: O(m(n+m)) per dense pivot; sparse LU in production

Tableau form: add a slack per row, Ax + s = b. Keep a basis (initially the slacks) and the tableau B⁻¹[A I | b] plus a reduced-cost row z − c·x.

1. Price: pick an entering column j with negative reduced cost. Dantzig's rule takes the most negative.
2. Ratio test: among rows i with a_ij > 0, pick the leaving row with the smallest b_i / a_ij. No such row means the LP is unbounded.
3. Pivot: scale the row, eliminate column j from all other rows (including the objective row).
4. Repeat until no reduced cost is negative: optimal.

• Degeneracy: ties in the ratio test give pivots that do not move; cycling is possible. Bland's rule (lowest index in both choices) provably terminates but is slow, so use it only after a stall.
• Phase I when some b_i < 0: add an artificial x0 with column −1, maximise −x0; one pivot of x0 into the most negative row yields a feasible basis. Final objective < 0 means infeasible.
• Revised simplex stores only B⁻¹ (as an LU factorisation) and computes per iteration: y = c_B B⁻¹ (BTRAN), reduced costs c_j − y·a_j, and the entering column B⁻¹ a_j (FTRAN).
• Dual simplex runs the same machinery on the dual: pick the leaving row first (most negative b_i), then the entering column by a ratio test on reduced costs. It stays dual feasible, so it restarts cheaply after adding a row or tightening a bound.

Code

Dense two-phase tableau, max c·x, Ax ≤ b, x ≥ 0, b of any sign. Compiled and tested against Beale's cycling example and brute-force vertex enumeration.

#[derive(Debug, PartialEq)]
pub enum Status { Optimal { x: Vec<f64>, value: f64 }, Infeasible, Unbounded }

const EPS: f64 = 1e-9;

struct Tab { t: Vec<Vec<f64>>, basis: Vec<usize>, cols: usize } // last row = objective z - c·x, last col = rhs

impl Tab {
    fn pivot(&mut self, r: usize, c: usize) {
        let p = self.t[r][c];
        for v in self.t[r].iter_mut() { *v /= p; }
        let pr = self.t[r].clone();
        for (i, row) in self.t.iter_mut().enumerate() {
            if i != r && row[c].abs() > 0.0 {
                let f = row[c];
                for (v, w) in row.iter_mut().zip(&pr) { *v -= f * w; }
            }
        }
        self.basis[r] = c;
    }
    /// Maximize; only columns < `allowed` may enter. Dantzig rule, Bland's rule once stalling.
    fn run(&mut self, allowed: usize) -> bool {
        let (m, rhs) = (self.basis.len(), self.cols);
        let mut stall = 0;
        loop {
            let obj = &self.t[m];
            let enter = if stall > 50 {
                (0..allowed).find(|&j| obj[j] < -EPS) // Bland: lowest index
            } else {
                (0..allowed).filter(|&j| obj[j] < -EPS)
                    .min_by(|&a, &b| obj[a].partial_cmp(&obj[b]).unwrap())
            };
            let Some(c) = enter else { return true };
            let mut leave: Option<(usize, f64)> = None; // ratio test, ties -> lowest basis index
            for i in 0..m {
                if self.t[i][c] > EPS {
                    let ratio = self.t[i][rhs] / self.t[i][c];
                    let better = match leave {
                        None => true,
                        Some((l, r)) => ratio < r - EPS || (ratio < r + EPS && self.basis[i] < self.basis[l]),
                    };
                    if better { leave = Some((i, ratio)); }
                }
            }
            let Some((r, ratio)) = leave else { return false }; // unbounded
            stall = if ratio < EPS { stall + 1 } else { 0 };
            self.pivot(r, c);
        }
    }
}

/// maximize c·x subject to A x <= b, x >= 0 (b may be negative).
pub fn simplex(a: &[Vec<f64>], b: &[f64], c: &[f64]) -> Status {
    let (m, n) = (a.len(), c.len());
    let (x0, rhs) = (n + m, n + m + 1); // columns: x, slacks, artificial x0, rhs
    let mut t = vec![vec![0.0; rhs + 1]; m + 1];
    for i in 0..m {
        t[i][..n].copy_from_slice(&a[i]);
        t[i][n + i] = 1.0;
        t[i][x0] = -1.0;
        t[i][rhs] = b[i];
    }
    let mut tab = Tab { t, basis: (n..n + m).collect(), cols: rhs };
    if let Some((k, _)) = b.iter().enumerate().filter(|(_, &v)| v < 0.0)
        .min_by(|x, y| x.1.partial_cmp(y.1).unwrap())
    {
        // Phase I: maximize -x0; one pivot of x0 into the most negative row makes the basis feasible.
        tab.t[m][x0] = 1.0;
        tab.pivot(k, x0);
        tab.run(rhs);
        if tab.t[m][rhs] < -1e-7 { return Status::Infeasible; }
        if let Some(r) = tab.basis.iter().position(|&v| v == x0) { // degenerate: drive x0 out
            if let Some(j) = (0..x0).find(|&j| tab.t[r][j].abs() > EPS) { tab.pivot(r, j); }
        }
        tab.t[m].iter_mut().for_each(|v| *v = 0.0);
    }
    // Phase II objective: z - c·x = 0, then eliminate basic columns from the objective row.
    for j in 0..n { tab.t[m][j] = -c[j]; }
    for r in 0..m {
        let f = tab.t[m][tab.basis[r]];
        if f != 0.0 { let row = tab.t[r].clone(); for (v, w) in tab.t[m].iter_mut().zip(&row) { *v -= f * w; } }
    }
    if !tab.run(x0) { return Status::Unbounded; }
    let mut x = vec![0.0; n];
    for (r, &v) in tab.basis.iter().enumerate() { if v < n { x[v] = tab.t[r][rhs]; } }
    Status::Optimal { x, value: tab.t[m][rhs] }
}

Not for real work: no scaling, no sparse algebra, dense O(m·n) memory.

Complexity and constants

• Dense tableau: 8·(m+1)(n+m+2) bytes; m=n=1,000 is about 16 MB, m=n=10,000 about 1.6 GB. One pivot touches all of it: about 2(m+1)(n+m+2) flops, streamed from memory.
• Measured by me (code above, random dense LP with positive data, release build, Apple Silicon): m=n=600 took 370 pivots, 60 ms, about 0.16 ms per pivot (about 9 Gflop/s). Pivots were 0.27·m, 0.49·m and 0.62·m for m=100, 300, 600. Real LPs vary widely.
• Revised simplex per iteration: one BTRAN, one FTRAN, pricing over nonzero columns; cost is in the nonzeros of the LU factors, not m·n.
• Worst case: Klee-Minty cube in n dimensions has 2ⁿ vertices and Dantzig's rule visits all of them.

Efficiency

Bound: memory bandwidth on the tableau; sparse index chasing and LU solves in the revised form.

• Revised form with sparse LU: never build the tableau; touch only nonzeros. Update the factors per pivot (Forrest-Tomlin) and refactor periodically to limit fill and error.
• Pricing: steepest-edge or Devex weights choose better entering columns than Dantzig and cut iteration counts (Forrest-Goldfarb 1992; Harris 1973); partial pricing scans a subset of columns when m ≪ n.
• Warm start: reuse the previous basis after a change; usually a handful of pivots instead of a full solve.
• Presolve: removes rows, columns and fixed variables before any pivot; the cheapest pivot is the one never made.
• Scaling to similar magnitudes keeps pivots stable.
• Cheaper alternative: if 1e-4 relative accuracy suffices on a huge sparse LP, first-order methods (PDHG) use only matrix-vector products and trade accuracy for much less work per iteration.

Implementation notes

• Layout: store the tableau row-major so the pivot row streams and each row update is an axpy that auto-vectorises; skip rows with a_ic = 0 (as the code does).
• Use tolerances, not zero: pivot element ≳ 1e-9 relative, reduced cost −1e-9, feasibility 1e-7. Harris' two-pass ratio test picks the largest pivot among near-ties for stability.
• Revised simplex: keep A column-wise (CSC) for pricing, rows for dual pricing; LU with Markowitz ordering (LU); refactor every order of 50-100 updates, or on numerical trouble.
• Bounded variables (l ≤ x ≤ u) are handled natively: nonbasic variables sit at a bound and may flip bound without a basis change. Do not add rows for bounds.
• Perturb costs or bounds slightly only if stalling shows up.

Pitfalls

• Absolute zero comparisons: tiny pivots amplify error and silently corrupt the tableau.
• Cycling under Dantzig on degenerate LPs (assignment, network structures): needs Bland, perturbation or a lexicographic rule.
• Badly scaled data (coefficients spanning 1e-6 to 1e6) makes optimality and feasibility tolerances meaningless.
• Reading a Phase I result as exact: "infeasible" means residual above tolerance.
• Using a dense tableau on a sparse model: fill-in turns a 1 MB problem into gigabytes.

Variants and alternatives

In production

HiGHS (primal and dual revised simplex; the default linprog method in SciPy), Clp and SoPlex (the LP engines of CBC and SCIP respectively), GLPK; commercial Gurobi, CPLEX, Xpress and Mosek offer primal and dual simplex as well. Uses: diet and blending (Stigler's diet problem was an early test of the method, Dantzig); power-system dispatch and capacity expansion (PyPSA); scientific computing (SciPy linprog); and the LP relaxations inside every MIP solver, hence logistics, scheduling and packing models.

Implementations

Sources

• George B. Dantzig, The Diet Problem, Interfaces 20(4), 1990. paper
• George B. Dantzig, Linear Programming and Extensions, Princeton University Press 1963.
• Victor Klee, George J. Minty, How good is the simplex algorithm?, in Inequalities III, Academic Press 1972, pp. 159-175.
• Robert G. Bland, New finite pivoting rules for the simplex method, Mathematics of Operations Research 2(2), 1977. paper
• Daniel A. Spielman, Shang-Hua Teng, Smoothed analysis of algorithms: why the simplex algorithm usually takes polynomial time, JACM 51(3), 2004. arXiv
• John J. Forrest, Donald Goldfarb, Steepest-edge simplex algorithms for linear programming, Mathematical Programming 57, 1992. paper
• Paula M. J. Harris, Pivot selection methods of the Devex LP code, Mathematical Programming 5, 1973. paper
• J. J. H. Forrest, J. A. Tomlin, Updated triangular factors of the basis to maintain sparsity in the product form simplex method, Mathematical Programming 2, 1972. paper
• Qi Huangfu, J. A. Julian Hall, Parallelizing the dual revised simplex method, Mathematical Programming Computation 10, 2018. paper
• Robert J. Vanderbei, Linear Programming: Foundations and Extensions, Springer, 4th ed. 2014 (auxiliary-problem Phase I). book