Interior-point method (primal-dual, Mehrotra)

Deep dive

How it works

Parts:

• Iterate (x, y, s) — primal, multipliers, dual slacks; x, s > 0, may be infeasible.
• Residuals — rp = b − Ax, rd = c − Aᵀy − s, complementarity xᵢsᵢ = μ.
• Normal matrix A D Aᵀ, D = diag(xᵢ/sᵢ) — SPD, factored by Cholesky, sparse via direct solver and one fill-reducing ordering.
• Step control — fraction-to-boundary, Mehrotra's rule for μ.

1. Newton on Ax = b, Aᵀy + s = c, xᵢsᵢ = μ: A dx = rp, Aᵀdy + ds = rd, S dx + X ds = rc. Eliminate: A D Aᵀ dy = rp − A S⁻¹(rc − X rd), ds = rd − Aᵀdy, dx = S⁻¹(rc − X ds).
2. Predictor: rc = −XSe (σ = 0); μ_aff after its step.
3. Centering: σ = (μ_aff/μ)³.
4. Corrector: rc = σμe − XSe − dxₐdsₐ; same factorization, one more solve.
5. Step length: 0.995 × the largest t keeping x, s > 0, capped at 1, separate for x and (y, s).
6. Stop when residuals and xᵀs are small.

 start x=s=e -> predictor -> corrector -> ... mu -> 0 (middle of optimal face)

Infeasible start: x = s = e, y = 0, no phase I. Infeasibility is certified by the homogeneous self-dual embedding (Ye, Todd, Mizuno 1994; Xu, Hung, Ye 1996): an always-feasible LP whose solution is an optimum or a certificate. Crossover pivots to a vertex basis for simplex or branch and bound. QP adds Q; prefer the augmented system [Q + X⁻¹S, Aᵀ; A, 0] with LDLᵀ.

Code

Dense LP, Mehrotra, dense Cholesky; compiled and run. Tests: Wyndor LP (optimum −36) against vertex enumeration, and a seeded random 20 × 50 LP built around a known primal-dual interior point.

type Mat = Vec<Vec<f64>>;
fn chol_solve(mut m: Mat, r: &[f64]) -> Vec<f64> {
    let n = r.len();
    for j in 0..n {
        let mut d = m[j][j];
        for k in 0..j { d -= m[j][k] * m[j][k]; }
        let d = if d > 1e-14 { d.sqrt() } else { 1e7 };
        m[j][j] = d;
        for i in j + 1..n {
            let mut s = m[i][j];
            for k in 0..j { s -= m[i][k] * m[j][k]; }
            m[i][j] = s / d;
        }
    }
    let mut z = r.to_vec();
    for i in 0..n { for k in 0..i { z[i] -= m[i][k] * z[k]; } z[i] /= m[i][i]; }
    for i in (0..n).rev() { for k in i + 1..n { z[i] -= m[k][i] * z[k]; } z[i] /= m[i][i]; }
    z
}
fn dot(a: &[f64], b: &[f64]) -> f64 { a.iter().zip(b).map(|(p, q)| p * q).sum() }
fn at_y(a: &Mat, y: &[f64], n: usize) -> Vec<f64> { (0..n).map(|j| (0..a.len()).map(|i| a[i][j] * y[i]).sum()).collect() }
fn a_x(a: &Mat, x: &[f64]) -> Vec<f64> { a.iter().map(|row| dot(row, x)).collect() }
fn norm(v: &[f64]) -> f64 { dot(v, v).sqrt() }
fn step(v: &[f64], dv: &[f64]) -> f64 {
    (0.995 * v.iter().zip(dv).filter(|(_, &d)| d < 0.0).map(|(&v, &d)| -v / d).fold(f64::INFINITY, f64::min)).min(1.0)
}
/// min cᵀx, Ax = b, x >= 0 -> (x, y, s, iterations)
fn ipm(a: &Mat, b: &[f64], c: &[f64], tol: f64) -> (Vec<f64>, Vec<f64>, Vec<f64>, usize) {
    let (m, n) = (a.len(), c.len());
    let (mut x, mut y, mut s) = (vec![1.0; n], vec![0.0; m], vec![1.0; n]);
    for it in 0..100 {
        let ax = a_x(a, &x);
        let rp: Vec<f64> = (0..m).map(|i| b[i] - ax[i]).collect();        
        let aty = at_y(a, &y, n);
        let rd: Vec<f64> = (0..n).map(|j| c[j] - aty[j] - s[j]).collect();
        let mu = dot(&x, &s) / n as f64;
        if norm(&rp) <= tol * (1.0 + norm(b)) && norm(&rd) <= tol * (1.0 + norm(c)) && n as f64 * mu <= tol {
            return (x, y, s, it);
        }
        let d: Vec<f64> = (0..n).map(|j| x[j] / s[j]).collect();
        let mut nm = vec![vec![0.0; m]; m];
        for i in 0..m { for k in 0..=i {
            let v: f64 = (0..n).map(|j| a[i][j] * d[j] * a[k][j]).sum();
            nm[i][k] = v; nm[k][i] = v;
        } }
        let solve = |rc: &[f64]| {
            let t: Vec<f64> = (0..n).map(|j| (rc[j] - x[j] * rd[j]) / s[j]).collect(); // S⁻¹(rc - X rd)
            let at = a_x(a, &t);
            let r: Vec<f64> = (0..m).map(|i| rp[i] - at[i]).collect();
            let dy = chol_solve(nm.clone(), &r);
            let aty = at_y(a, &dy, n);
            let ds: Vec<f64> = (0..n).map(|j| rd[j] - aty[j]).collect();
            let dx: Vec<f64> = (0..n).map(|j| (rc[j] - x[j] * ds[j]) / s[j]).collect();
            (dx, dy, ds)
        };
        let rc: Vec<f64> = (0..n).map(|j| -x[j] * s[j]).collect();
        let (dxa, _, dsa) = solve(&rc); // predictor
        let (ap, ad) = (step(&x, &dxa), step(&s, &dsa));
        let mu_aff = (0..n).map(|j| (x[j] + ap * dxa[j]) * (s[j] + ad * dsa[j])).sum::<f64>() / n as f64;
        let sigma = (mu_aff / mu).powi(3);
        // centering + second-order corrector
        let rc: Vec<f64> = (0..n).map(|j| sigma * mu - x[j] * s[j] - dxa[j] * dsa[j]).collect();
        let (dx, dy, ds) = solve(&rc);
        let (ap, ad) = (step(&x, &dx), step(&s, &ds));
        for j in 0..n { x[j] += ap * dx[j]; s[j] += ad * ds[j]; }
        for i in 0..m { y[i] += ad * dy[i]; }
    }
    panic!("no convergence");
}
fn main() {
    let a: Mat = vec![vec![1., 0., 1., 0., 0.], vec![0., 2., 0., 1., 0.], vec![3., 2., 0., 0., 1.]];
    let (b, c) = (vec![4., 12., 18.], vec![-3., -5., 0., 0., 0.]);
    let (x, _, _, it) = ipm(&a, &b, &c, 1e-10);
    let obj = dot(&c, &x); println!("tiny LP: obj {obj:.9} x {:.6?} in {it} iterations", &x[..2]);
    assert!((obj + 36.0).abs() < 1e-6);
    assert!((brute(&a, &b, &c) - obj).abs() < 1e-6);
    let mut seed = 12345u64;
    let mut rnd = || { seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); (seed >> 11) as f64 / (1u64 << 53) as f64 };
    let (m, n) = (20, 50); let a: Mat = (0..m).map(|_| (0..n).map(|_| 2.0 * rnd() - 1.0).collect()).collect();
    let (x0, y0, s0): (Vec<f64>, Vec<f64>, Vec<f64>) = ((0..n).map(|_| rnd() + 0.1).collect(), (0..m).map(|_| 2.0 * rnd() - 1.0).collect(), (0..n).map(|_| rnd() + 0.1).collect());
    let b = a_x(&a, &x0);
    let c: Vec<f64> = at_y(&a, &y0, n).iter().zip(&s0).map(|(p, q)| p + q).collect();
    let (x, y, s, it) = ipm(&a, &b, &c, 1e-10);
    let rp = norm(&a_x(&a, &x).iter().zip(&b).map(|(p, q)| p - q).collect::<Vec<_>>());
    let rd = norm(&at_y(&a, &y, n).iter().zip(&s).zip(&c).map(|((p, q), r)| r - p - q).collect::<Vec<_>>());
    let gap = dot(&x, &s);
    println!("random LP {m}x{n}: obj {:.9} primal res {rp:.1e} dual res {rd:.1e} gap {gap:.1e} in {it} iterations", dot(&c, &x));
    assert!(rp < 1e-8 && rd < 1e-8 && gap < 1e-8);
}
/// Best objective over all feasible basic solutions.
fn brute(a: &Mat, b: &[f64], c: &[f64]) -> f64 {
    let (m, n) = (a.len(), c.len());
    let mut best = f64::INFINITY;
    for mask in 0u32..(1 << n) {
        if mask.count_ones() as usize != m { continue; }
        let cols: Vec<usize> = (0..n).filter(|j| mask >> j & 1 == 1).collect();
        let mut t: Mat = (0..m).map(|i| { let mut r: Vec<f64> = cols.iter().map(|&j| a[i][j]).collect(); r.push(b[i]); r }).collect();
        let mut ok = true;
        for k in 0..m {
            let p = (k..m).max_by(|&i, &j| t[i][k].abs().partial_cmp(&t[j][k].abs()).unwrap()).unwrap();
            if t[p][k].abs() < 1e-12 { ok = false; break; }
            t.swap(k, p);
            for i in 0..m { if i != k { let f = t[i][k] / t[k][k]; for j in k..=m { let v = t[k][j]; t[i][j] -= f * v; } } }
        }
        if !ok { continue; }
        let xb: Vec<f64> = (0..m).map(|i| t[i][m] / t[i][i]).collect();
        if xb.iter().all(|&v| v >= -1e-9) { best = best.min(cols.iter().zip(&xb).map(|(&j, &v)| c[j] * v).sum()); }
    }
    best
}

Measured: Wyndor 8 iterations; random LP 10 (primal residual 1.2e-12, dual 3.7e-15, xᵀs 7.1e-12), cold started.

Complexity and constants

• Per iteration: forming A D Aᵀ O(m²n) dense; Cholesky m³/3 flops; each solve 2m².
• Iterations: 8 and 10 measured here; "tens, size-independent" is general knowledge.

Efficiency

• Cost: one sparse Cholesky per iteration.
• Levers:
  • One ordering and symbolic analysis — nnz(L) far below natural order — always.
  • Predictor and corrector share the factor — second solve O(nnz(L)) (Mehrotra 1992) — always.
  • Presolve — fewer rows, fewer dense columns — before every solve.
  • Preconditioned CG instead of factoring — memory — L does not fit (HiGHS IPX).
• Naive → standard → state of the art: dense normal equations · O(m³) per iteration; sparse Cholesky + Mehrotra · few dozen factorizations; plus presolve and parallel supernodal factors.

Implementation notes

• Late D spans many orders of magnitude: pin tiny pivots (as the code does) or use the augmented form.

Pitfalls

• Check both residuals and the gap, not the objective alone.
• Splitting free variables x = x⁺ − x⁻: both parts drift off to infinity.
• Expecting a vertex: branching needs crossover.

Variants and alternatives

In production

• Logistics, telecom, energy (general knowledge, unchecked): large network-design, capacity and dispatch LPs use barrier plus crossover.
• Finance, machine learning: portfolio and SVM training are convex QPs, the target of MOSEK and Clarabel.
• Robotics: model predictive control is a QP per step (unverified).
• HiGHS ships IPX and HiPO; Ipopt does nonlinear.

Implementations

Sources

• Karmarkar, A new polynomial-time algorithm for linear programming, Combinatorica 4, 1984. paper
• Mehrotra, On the implementation of a primal-dual interior point method, SIAM J. Optimization 2(4), 1992. paper
• Ye, Todd, Mizuno, An O(√nL)-iteration homogeneous and self-dual linear programming algorithm, MOR 19(1), 1994. paper
• Xu, Hung, Ye, A simplified homogeneous and self-dual linear programming algorithm and its implementation, Ann. Oper. Res. 62, 1996. paper
• Goulart, Chen, Clarabel, 2024. arXiv · code