Exact Linear Programming

symplex::linprog solves linear programs over ℚ: a two-phase dense simplex running on Ratio<BigInt>, with Bland’s rule after the first degenerate step so it cannot cycle. Optima, shadow prices and infeasibility certificates are exact — there are no tolerances and no “numerically infeasible” verdicts. That is what makes it useful as the engine behind certificate searches (Farkas lemmas, Positivstellensatz-style combinations, Carathéodory decompositions), where a floating-point solver can only say “probably”.

The module works on plain Vec<Q> data (Q = Ratio<BigInt>), with two small constructors: qi(n) for an integer and q(n, d) for n/d. LpSolution::x_ex(&ctx) converts an optimum back into Ex rationals, and linprog_matrix accepts Matrix data directly.

The builder

LpProblem::minimize(c) / maximize(c) start a program in c.len() variables; .le(row, rhs), .ge(row, rhs), .eq(row, rhs) add constraint rows; .bounds(j, Bounds::closed(lo, hi)) (or Bounds::at_least(lo), Bounds::at_most(hi)) and .free(j) change a variable’s bounds from the default 0 ≤ xⱼ < ∞; .solve() returns Result<LpSolution>.

use symplex::prelude::*;
use symplex::linprog::{LpProblem, Q, qi};

/// `Ratio<BigInt>` displays as `3/2`, but `{:?}` on a `Vec<Q>` is verbose — format by hand.
fn show(v: &[Q]) -> String {
    format!("({})", v.iter().map(|x| x.to_string()).collect::<Vec<_>>().join(", "))
}

fn main() {
    // max 5x + 4y  s.t.  6x + 4y ≤ 24,  x + 2y ≤ 6,  −x + y ≤ 1,  y ≤ 2,  x, y ≥ 0
    let sol = LpProblem::maximize(vec![qi(5), qi(4)])
        .le(vec![qi(6), qi(4)], qi(24))
        .le(vec![qi(1), qi(2)], qi(6))
        .le(vec![qi(-1), qi(1)], qi(1))
        .le(vec![qi(0), qi(1)], qi(2))
        .solve()
        .unwrap();
    println!("{:?}", sol.status);                     // Optimal
    println!("{}", show(&sol.x));                     // (3, 3/2)
    println!("{}", sol.objective.clone().unwrap());   // 21
    println!("{}", show(&sol.duals));                 // (3/4, 1/2, 0, 0)
    println!("{}", sol.is_optimal());                 // true
    let ctx = Context::new();
    println!("{:?}", sol.x_ex(&ctx));                 // [Ex(3), Ex(3/2)]
}

The remaining examples on this page reuse the show helper.

Statuses

solve() only returns Err for malformed input — no variables, a row of the wrong length, bounds on a variable that does not exist. The three mathematical outcomes are values of LpStatus:

use symplex::Bounds;
use symplex::linprog::{LpProblem, qi};

fn main() {
    let sol = LpProblem::maximize(vec![qi(1), qi(1)]).le(vec![qi(1), qi(-1)], qi(1)).solve().unwrap();
    println!("{:?} {:?} {:?}", sol.status, sol.objective, sol.farkas);   // Unbounded None None

    println!("{}", LpProblem::maximize(vec![qi(1), qi(1)]).le(vec![qi(1)], qi(1)).solve().unwrap_err());
    // linprog: invalid argument: constraint 0 has 1 coefficients but there are 2 variables
    println!("{}", LpProblem::maximize(vec![qi(1)]).bounds(3, Bounds::free()).solve().unwrap_err());
    // linprog: invalid argument: bounds were set for variable 3 but there are only 1 variables

    // Contradictory bounds: infeasible, but there is no constraint certificate to give.
    let bad = LpProblem::minimize(vec![qi(1)]).bounds(0, Bounds::closed(qi(3), qi(1))).solve().unwrap();
    println!("{:?} {:?}", bad.status, bad.farkas);                       // Infeasible None
}

Bounds and free variables

// min x − y  s.t.  x + y ≤ 3,  −2 ≤ x,  0 ≤ y ≤ 1
let sol = LpProblem::minimize(vec![qi(1), qi(-1)])
    .le(vec![qi(1), qi(1)], qi(3))
    .bounds(0, Bounds::at_least(qi(-2)))
    .bounds(1, Bounds::closed(qi(0), qi(1)))
    .solve()
    .unwrap();
println!("{:?} x* = {} objective {} duals {}", sol.status, show(&sol.x), sol.objective.unwrap(), show(&sol.duals));
// Optimal x* = (-2, 1) objective -3 duals (0)      — the constraint is slack, so its price is 0

// A free variable: min x  s.t.  2x ≥ −5
let free = LpProblem::minimize(vec![qi(1)]).free(0).ge(vec![qi(2)], qi(-5)).solve().unwrap();
println!("{:?} x* = {} objective {}", free.status, show(&free.x), free.objective.unwrap());
// Optimal x* = (-5/2) objective -5/2

Duals and complementary slackness

When the status is Optimal, duals holds one shadow price yᵢ per constraint, in insertion order: the rate of change of the optimal objective value of the problem as posed with respect to bᵢ. The sign conventions that follow from that definition are, quoting the module documentation:

for a minimisation yᵢ ≤ 0 on ≤ rows and yᵢ ≥ 0 on ≥ rows (the signs flip for a maximisation), yᵢ is free on = rows, and with the reduced costs r = c − Aᵀy:

• complementary slackness: yᵢ·(aᵢ·x* − bᵢ) = 0 for every row;
• rⱼ = 0 unless x*ⱼ sits at a finite bound (for a minimisation rⱼ ≥ 0 at a lower bound and rⱼ ≤ 0 at an upper bound; reversed for a maximisation);
• strong duality: cᵀx* = yᵀb + Σⱼ rⱼ x*ⱼ, which reduces to cᵀx* = yᵀb under the default bounds x ≥ 0.

All of these are identities you can check with exact arithmetic:

use symplex::linprog::{LpProblem, Q, qi};

fn show(v: &[Q]) -> String {
    format!("({})", v.iter().map(|x| x.to_string()).collect::<Vec<_>>().join(", "))
}

fn main() {
    let rows = [
        (vec![qi(6), qi(4)], qi(24)),
        (vec![qi(1), qi(2)], qi(6)),
        (vec![qi(-1), qi(1)], qi(1)),
        (vec![qi(0), qi(1)], qi(2)),
    ];
    let mut p = LpProblem::maximize(vec![qi(5), qi(4)]);
    for (row, rhs) in &rows {
        p = p.le(row.clone(), rhs.clone());
    }
    let sol = p.solve().unwrap();

    // Strong duality under x ≥ 0: cᵀx* = yᵀb.
    let ytb: Q = rows.iter().zip(&sol.duals).map(|((_, b), y)| b * y).sum();
    println!("cᵀx* = {}   yᵀb = {ytb}", sol.objective.clone().unwrap());   // cᵀx* = 21   yᵀb = 21

    // Complementary slackness: only binding rows have a non-zero price.
    for (i, (row, b)) in rows.iter().enumerate() {
        let slack: Q = row.iter().zip(&sol.x).map(|(a, x)| a * x).sum::<Q>() - b;
        println!("row {i}: a·x* − b = {slack:>4},  y = {:>3},  y·slack = {}",
            sol.duals[i], &slack * &sol.duals[i]);
    }
    // row 0: a·x* − b =    0,  y = 3/4,  y·slack = 0
    // row 1: a·x* − b =    0,  y = 1/2,  y·slack = 0
    // row 2: a·x* − b = -5/2,  y =   0,  y·slack = 0
    // row 3: a·x* − b = -1/2,  y =   0,  y·slack = 0

    // y₀ = 3/4 is ∂(optimum)/∂b₀: raising b₀ from 24 to 25 adds exactly 3/4.
    let sol2 = LpProblem::maximize(vec![qi(5), qi(4)])
        .le(vec![qi(6), qi(4)], qi(25))
        .le(vec![qi(1), qi(2)], qi(6))
        .le(vec![qi(-1), qi(1)], qi(1))
        .le(vec![qi(0), qi(1)], qi(2))
        .solve()
        .unwrap();
    println!("b₀ = 25 → objective {}", sol2.objective.unwrap());          // 87/4

    // Minimisation with ≥ rows: prices are ≥ 0.
    let m = LpProblem::minimize(vec![qi(1), qi(1)])
        .ge(vec![qi(1), qi(2)], qi(1))
        .ge(vec![qi(3), qi(1)], qi(1))
        .solve()
        .unwrap();
    println!("x* = {}, objective {}, duals {}", show(&m.x), m.objective.unwrap(), show(&m.duals));
    // x* = (1/5, 2/5), objective 3/5, duals (2/5, 1/5)
}

Farkas certificates

When the status is Infeasible, farkas is a vector y with one entry per constraint that proves infeasibility. From the module documentation:

yᵢ ≥ 0 on ≤ rows, yᵢ ≤ 0 on ≥ rows, free on = rows, such that, with g = Aᵀy,

Σⱼ  inf { gⱼ·xⱼ : lⱼ ≤ xⱼ ≤ uⱼ }   >   yᵀb

where every infimum is finite (gⱼ > 0 ⇒ lⱼ finite, gⱼ < 0 ⇒ uⱼ finite, gⱼ = 0 contributes 0). Any feasible x would satisfy (Aᵀy)·x ≤ yᵀb, so the inequality proves that none exists. With no finite bounds this is the textbook form Aᵀy = 0, yᵀb < 0.

Under the default bounds x ≥ 0 the infima are all 0, so the certificate reads Aᵀy ≥ 0 and yᵀb < 0:

use symplex::linprog::{LpProblem, LpStatus, qi};
use num_traits::Signed;

fn main() {
    // x + y ≤ 1  and  x + y ≥ 2  cannot both hold.
    let sol = LpProblem::minimize(vec![qi(1), qi(1)])
        .le(vec![qi(1), qi(1)], qi(1))
        .ge(vec![qi(1), qi(1)], qi(2))
        .solve()
        .unwrap();
    assert_eq!(sol.status, LpStatus::Infeasible);
    let y = sol.farkas.clone().unwrap();
    println!("y = ({}, {})", y[0], y[1]);                    // y = (1, -1):  y₀ ≥ 0 on the ≤ row, y₁ ≤ 0 on the ≥ row
    let g = &y[0] + &y[1];                                   // both columns of A are (1, 1)
    let ytb = &y[0] + &(&y[1] * qi(2));
    println!("Aᵀy = ({g}, {g}),  yᵀb = {ytb}");              // Aᵀy = (0, 0),  yᵀb = -1
    assert!(!g.is_negative() && ytb.is_negative());
    println!("{} {}", sol.x.len(), sol.duals.is_empty());    // 0 true   (no point, no prices)
}

In words: adding the first row to −1 times the second gives 0 ≤ −1. The cookbook shows a Farkas vector being read as a linear functional that separates a polynomial from a cone.

feasible_nonneg: is b a non-negative combination?

feasible_nonneg(&a_eq, &b_eq) answers “is there an x ≥ 0 with A·x = b?” — Ok(Some(x)) with a witness, or Ok(None). It is the query behind most certificate searches and is exact even when the data have denominators like 1/3 and 1/7.

use symplex::linprog::{Q, feasible_nonneg, q, qi};

fn show(v: &[Q]) -> String {
    format!("({})", v.iter().map(|x| x.to_string()).collect::<Vec<_>>().join(", "))
}

fn main() {
    // μ ≥ 0 with  μ₁/3 + μ₂/7 + 2μ₃/5 = 1  and  μ₁ + μ₂ + μ₃ = 4
    let a = vec![vec![q(1, 3), q(1, 7), q(2, 5)], vec![qi(1), qi(1), qi(1)]];
    let b = vec![qi(1), qi(4)];
    match feasible_nonneg(&a, &b).unwrap() {
        Some(mu) => {
            println!("μ = {}", show(&mu));                               // μ = (0, 7/3, 5/3)
            for (row, rhs) in a.iter().zip(&b) {
                let lhs: Q = row.iter().zip(&mu).map(|(c, m)| c * m).sum();
                assert_eq!(&lhs, rhs);
            }
            println!("A·μ = b exactly");
        }
        None => println!("no non-negative combination"),
    }

    // x + y = −1 has no non-negative solution.
    println!("{:?}", feasible_nonneg(&[vec![qi(1), qi(1)]], &[qi(-1)]).unwrap());   // None

    // "Is (2, 3, 3) in the cone spanned by (1,0,1), (0,1,1), (1,1,0)?"
    // Columns are the generators, so build the rows by transposing.
    let cols = [[qi(1), qi(0), qi(1)], [qi(0), qi(1), qi(1)], [qi(1), qi(1), qi(0)]];
    let target = [qi(2), qi(3), qi(3)];
    let rows: Vec<Vec<Q>> = (0..3).map(|i| cols.iter().map(|c| c[i].clone()).collect()).collect();
    println!("{:?}", feasible_nonneg(&rows, &target).unwrap().map(|v| show(&v)));   // Some("(1, 2, 1)")
}

Note the orientation: feasible_nonneg takes rows of A. When your generators are naturally columns — a list of vectors, or polynomials laid out by Poly::coefficient_matrix — use nonneg_combination(&vectors, &target) instead, which asks the cone-membership question directly. Both have a certified form: feasible_nonneg_certified and nonneg_combination return a Feasibility, whose Infeasible { farkas } variant carries the separating vector y (y·vⱼ ≥ 0 for every generator, y·target < 0), so there is no need to re-pose the system as an LpProblem to obtain the proof of impossibility.

use symplex::linprog::{Feasibility, Q, nonneg_combination, qi};

fn show(v: &[Q]) -> String {
    format!("({})", v.iter().map(|x| x.to_string()).collect::<Vec<_>>().join(", "))
}

fn main() {
    let cone = [vec![qi(1), qi(0), qi(1)], vec![qi(0), qi(1), qi(1)], vec![qi(1), qi(1), qi(0)]];
    match nonneg_combination(&cone, &[qi(2), qi(3), qi(3)]).unwrap() {
        Feasibility::Feasible(lambda) => println!("λ = {}", show(&lambda)),      // λ = (1, 2, 1)
        Feasibility::Infeasible { .. } => println!("outside the cone"),
    }
    // (1, 0, 0) is outside: it would need λ₁ + λ₃ = 1, λ₂ + λ₃ = 0, λ₁ + λ₂ = 0.
    match nonneg_combination(&cone, &[qi(1), qi(0), qi(0)]).unwrap() {
        Feasibility::Infeasible { farkas: Some(y) } => println!("separating y = {}", show(&y)),   // separating y = (-1, 1, 1)
        other => println!("{other:?}"),
    }
    // y·(1,0,1) = 0, y·(0,1,1) = 2, y·(1,1,0) = 0 are all ≥ 0, while y·(1,0,0) = −1 < 0.
}

SciPy-shaped linprog

linprog(c, a_ub, b_ub, a_eq, b_eq, bounds) minimises cᵀx subject to A_ub·x ≤ b_ub, A_eq·x = b_eq and per-variable bounds (a &[Bounds<Q>], one per variable; empty means x ≥ 0). Constraints are numbered ≤ rows first, then = rows — that is the order of duals and farkas.

// min −x − y   s.t.  x + 2y ≤ 4,  3x + y ≤ 6,  x, y ≥ 0
let sol = linprog(
    &[qi(-1), qi(-1)],
    &[vec![qi(1), qi(2)], vec![qi(3), qi(1)]],
    &[qi(4), qi(6)],
    &[],
    &[],
    &[],
)
.unwrap();
println!("{:?} x* = {} objective {}", sol.status, show(&sol.x), sol.objective.unwrap());
// Optimal x* = (8/5, 6/5) objective -14/5

Matrix input: linprog_matrix

linprog_matrix(objective, &c, a_ub, b_ub, a_eq, b_eq) takes Matrix data. Entries are constant-folded with eval() first, so 1 + 2 or 1/2 + 1/3 are fine; a symbol or π is rejected with a clear error rather than approximated. Bounds are the default x ≥ 0.

use symplex::prelude::*;
use symplex::linprog::{Objective, linprog_matrix};

fn main() {
    let ctx = Context::new();
    let c = matrix![ctx, [3], [2]];
    let a = Matrix::new(vec![
        vec![ctx.int(1), ctx.int(1)],
        vec![ctx.int(1), &ctx.int(1) + &ctx.int(2)],     // folded to 3
    ])
    .unwrap();
    let b = matrix![ctx, [4], [6]];
    let sol = linprog_matrix(Objective::Maximize, &c, Some(&a), Some(&b), None, None).unwrap();
    println!("x* = {:?}, objective {}", sol.x_ex(&ctx), sol.objective.unwrap());
    // x* = [Ex(4), Ex(0)], objective 12

    let x = ctx.symbol("x");
    let bad = Matrix::new(vec![vec![x, ctx.int(1)]]).unwrap();
    println!("{}", linprog_matrix(Objective::Minimize, &c, Some(&bad), Some(&matrix![ctx, [1]]), None, None).unwrap_err());
    // linprog_matrix: invalid argument: A_ub must contain only numeric literals; found `x`
}

Polytopes from half-spaces

symplex::polytope::Polytope (0.4) is a convex polyhedron {x ∈ ℚⁿ : aᵢ·x + bᵢ ≥ 0} with exact geometry built on the LP and on QMatrix: is_empty / any_point / bounding_box / is_bounded are LP calls; vertices solves every n × n sub-system exactly and keeps the points inside; volume (any dimension) is an exact facet decomposition around the vertex centroid; irredundant drops half-spaces that touch no vertex; split cuts by a hyperplane; from_exprs / to_exprs translate to and from affine Ex hypotheses, so a cell can go straight into prove_nonnegative_on_polyhedron.

use symplex::prelude::*;
use symplex::polytope::Polytope;
use symplex::linprog::{q, qi};

fn main() {
    let ctx = Context::new();
    let (r, t) = (ctx.symbol("r"), ctx.symbol("t"));
    // The unit box cut by t ≥ r and r + t ≤ 3/2.
    let cell = Polytope::from_exprs(
        &[r.clone(), 1 - &r, t.clone(), 1 - &t, &t - &r, ctx.rational(3, 2) - &r - &t],
        &[r.clone(), t.clone()],
    )
    .unwrap();
    let v: Vec<String> = cell.vertices().unwrap().iter().map(|p| format!("({}, {})", p[0], p[1])).collect();
    println!("{}", v.join(", "));                                 // (0, 0), (0, 1), (1/2, 1), (3/4, 3/4)
    println!("{}", cell.volume().unwrap());                       // 7/16
    println!("{}", cell.irredundant().unwrap().num_halfspaces()); // 5  (1 - r is implied)
    let halves = cell.split(&[qi(-1), qi(0)], q(1, 2));           // cut at r = 1/2: `1/2 - r ≥ 0` is the left piece
    println!("{} {}", halves.nonnegative.volume().unwrap(), halves.nonpositive.volume().unwrap());   // 3/8 1/16
    println!("{}", cell.contains(&[q(1, 4), q(1, 2)]));            // true
}

Everything is exact and every answer is a rational; the enumeration is O(C(m, n)) linear solves and the volume recursion visits every face, which is the right trade for the handful of cells a decision tree produces (dimension ≤ 5) and the wrong one for large polyhedra. Since 0.6.1 the enumeration runs in integer arithmetic (half-spaces scaled once, distinct hyperplanes only, containment as the sign of a·X + b·D), the vertex list is cached on the polytope, and volume hands each facet its own vertices instead of re-enumerating — about 10× on the vertex work. Ask is_full_dimensional() (one LP) rather than volume() > 0 when that is the question, and interior_point() for a point with positive slack everywhere; HalfSpace::normalized() is the key that identifies a candidate cut with its flip and rescalings. When the cell’s facets depend on a parameter, ParametricPolytope::new(&hyps, &vars, &j) holds the family and polytope_at / vertices_at / volume_at(&j_value) instantiate it exactly with a per-sample cache (0.5).

Cut scoring — hundreds of candidate planes per tree node, each needing the vertex sets of both pieces — should not re-enumerate: cell.clip(&h) derives them from the cached vertices in one pass. The vertex cache carries each vertex’s tight set (vertices_with_tight(): the indices of the half-spaces through it), and two vertices are joined by an edge exactly when the normals of their common tight half-spaces have rank n − 1 — for simple vertices that is just “share n − 1 indices”, for degenerate ones (a pyramid’s apex, a cube cut through a vertex, a face given by h ≥ 0 and h ≤ 0) an exact rank decides, so no crossing is invented or lost. Each edge with endpoints on opposite sides contributes the exact crossing vᵢ + h(vᵢ)/(h(vᵢ) − h(vⱼ))·(vⱼ − vᵢ); the result Clip { pos, neg, on } equals, as sets, with_halfspace(h).vertices() and its flipped counterpart. pos_is_full_dimensional() / neg_is_full_dimensional() answer the dimension question by a rank on the vertices (no LP, exact for bounded cells — also is_full_dimensional_from_vertices() on any polytope), and pos_polytope(&cell, &h) / neg_polytope return the pieces with their vertex cache already filled, so a tie-breaking volume() on them, or a further clip, enumerates nothing. In the downstream tree builder this replaced split + two vertices() per candidate and the same scores fell out byte for byte.

Performance

Since 0.3.5 the tableau uses integer pivoting: each constraint row is scaled once to clear its denominators, and every pivot then follows Bareiss’s fraction-free rule, so all entries stay integers sharing one common denominator (the current pivot, ±det B). Nothing in the inner loop computes a gcd; the ratio test and every sign test are integer comparisons. Results are identical to the rational tableau — same pivots, same optimum, same duals — because the same Dantzig/Bland choices are made on the same rational values, only represented differently.

Each pivot is still O(m·n) big-integer operations, but the constants are much smaller: in a release build a 40-row × 100-variable program went from 1.1 s to 41 ms, a 60 × 160 one from 2.8 s to 85 ms, and the degree-5 three-variable Handelman search in the certificates cookbook from 3.6 s to 0.9 s. Beyond a few hundred rows, or when the data are floating-point measurements to begin with, an exact solver is still the wrong tool — the numerical routines in Numerical Optimisation or an external LP library are. The pivot count is capped at 10 000 + 50·(m + n); exceeding it is reported as ComputationFailed, though Bland’s rule makes that a theoretical rather than a practical concern.

See cargo run --example exact_lp for the complete program.