Simplex method

The simplex method solves linear programs: it maximizes a linear function of real variables subject to linear inequalities. The feasible set is a convex polyhedron, and if an optimum exists it is attained at a vertex; the method starts at a vertex and repeatedly moves along an edge to a neighbouring vertex with a better objective value, until no neighbour is better^[2]. Each move is a pivot on a tableau of the constraints. The final tableau also contains a solution of the dual linear program, whose value equals the optimum and certifies it^[14]. The method is fast in practice, but with Dantzig's original pivoting rule it takes an exponential number of steps on the Klee–Minty cubes^[5], and on degenerate problems a pivoting rule such as Bland's is needed to guarantee termination^[6]. Dantzig published it in 1947–1951^[1]; it remains, with interior-point methods, one of the two standard ways of solving linear programs, and it is the core of the linear-arithmetic solvers in SMT solvers^[13].

§ 01

Linear programs

A linear program in standard inequality form is

maximize c^{T} x subject to A x \leq b, x \geq 0

with $A$ an $m \times n$ matrix. Adding a slack variable $s_{i} \geq 0$ to each inequality turns it into an equation $a_{i}^{T} x + s_{i} = b_{i}$ . A basic solution chooses $m$ of the $n + m$ variables as basic, sets the others to zero, and solves for the basic ones; it is feasible if they come out non-negative. The basic feasible solutions are exactly the vertices of the feasible polyhedron, and if the objective is bounded on it, some vertex is optimal^[4].

The feasible region of maximize 3x + 2y subject to x + y ≤ 4, x + 3y ≤ 9, x ≤ 3, x, y ≥ 0. The simplex method starts at the origin and follows two edges to the optimal vertex (3, 1), each step increasing the objective.

§ 02

The algorithm

When $b \geq 0$ , setting $x = 0$ and $s = b$ is a basic feasible solution with the slacks as basis. The method then repeats three steps^[4]. Pricing: the objective row of the tableau gives, for each non-basic variable $j$ , its reduced cost $c ˉ_{j}$ , the rate at which the objective changes if $x_{j}$ is increased from zero; if every $c ˉ_{j} \leq 0$ the current vertex is optimal. Ratio test: otherwise choose an entering variable with $c ˉ_{e} > 0$ and increase it until the first basic variable reaches zero:

r = ar g i : t_{i e} > 0 min \frac{β _{i}}{t _{i e}}

where $t_{i e}$ is the entering column and $β_{i}$ the current value of the $i$ th basic variable. If no $t_{i e}$ is positive, the entering variable can grow without bound and the program is unbounded. Pivot: the variable of row $r$ leaves the basis, and Gaussian elimination on the column $e$ rewrites the tableau in terms of the new basis. Geometrically, the pivot moves along an edge of the polyhedron to the next vertex.

Exact rational numbers.

(* Exact rationals on native integers, always in lowest terms with a
   positive denominator. *)
type t = { n : int; d : int }
let rec gcd a b = if b = 0 then abs a else gcd b (a mod b)
let make n d = if d = 0 then raise Division_by_zero else
  let g = gcd n d in let s = if d < 0 then -1 else 1 in { n = s * n / g; d = s * d / g }
let of_int n = { n; d = 1 }
let zero = of_int 0 and one = of_int 1
let add a b = make ((a.n * b.d) + (b.n * a.d)) (a.d * b.d)
let sub a b = make ((a.n * b.d) - (b.n * a.d)) (a.d * b.d)
let mul a b = make (a.n * b.n) (a.d * b.d)
let div a b = make (a.n * b.d) (a.d * b.n)
let compare a b = Stdlib.compare (a.n * b.d) (b.n * a.d)
let sign a = Stdlib.compare a.n 0
let to_string a = if a.d = 1 then string_of_int a.n else Printf.sprintf "%d/%d" a.n a.d

A tableau simplex method with Dantzig's and Bland's pivoting rules, which returns the primal and dual solutions.

(* The simplex method on a dense tableau, for
     maximize c.x  subject to  A x <= b,  x >= 0,  with b >= 0,
   so that the slack variables form a feasible starting basis. Variables
   0..n-1 are the original ones and n..n+m-1 the slacks. *)
type rule = Dantzig | Bland
 
type result = Optimal of Q.t * Q.t array * Q.t array | Unbounded | Cycled
 
let solve ?(rule = Bland) ?(trace = fun _ _ _ _ -> ()) ?(max_pivots = 100_000) a b c =
  let m = Array.length a and n = Array.length c in
  let t = Array.init m (fun i -> Array.init (n + m + 1) (fun j ->
    if j < n then a.(i).(j) else if j < n + m then (if j - n = i then Q.one else Q.zero) else b.(i))) in
  (* The objective row holds the reduced costs, and minus the objective
     value in the last column. *)
  let z = Array.init (n + m + 1) (fun j -> if j < n then c.(j) else Q.zero) in
  let basis = Array.init m (fun i -> n + i) in
  let seen = Hashtbl.create 64 in
  let pivot r e =
    let pv = t.(r).(e) in
    t.(r) <- Array.map (fun x -> Q.div x pv) t.(r);
    let elim row = let k = row.(e) in if Q.sign k <> 0 then Array.iteri (fun j x -> row.(j) <- Q.sub row.(j) (Q.mul k x)) t.(r) in
    Array.iteri (fun i row -> if i <> r then elim row) t;
    elim z;
    basis.(r) <- e
  in
  let rec loop k =
    (* The objective never decreases, so returning to an earlier basis
       means the method is cycling through degenerate pivots. *)
    let key = Array.to_list (Array.map string_of_int basis) |> String.concat "," in
    if Hashtbl.mem seen key || k >= max_pivots then Cycled
    else begin
      Hashtbl.replace seen key ();
      let candidates = List.filter (fun j -> Q.sign z.(j) > 0) (List.init (n + m) Fun.id) in
      match candidates with
      | [] ->
          let x = Array.make n Q.zero in
          Array.iteri (fun i v -> if v < n then x.(v) <- t.(i).(n + m)) basis;
          let y = Array.init m (fun i -> Q.sub Q.zero z.(n + i)) in
          Optimal (Q.sub Q.zero z.(n + m), x, y)
      | first :: _ ->
          let e = match rule with
            | Bland -> first
            | Dantzig -> List.fold_left (fun best j -> if Q.compare z.(j) z.(best) > 0 then j else best) first candidates in
          (* Ratio test; ties broken by the smallest basic variable (Bland)
             or the first row (Dantzig). *)
          let best = ref None in
          Array.iteri (fun i row ->
            if Q.sign row.(e) > 0 then begin
              let ratio = Q.div row.(n + m) row.(e) in
              match !best with
              | None -> best := Some (i, ratio)
              | Some (i', r') ->
                  let cmp = Q.compare ratio r' in
                  if cmp < 0 || (cmp = 0 && rule = Bland && basis.(i) < basis.(i')) then best := Some (i, ratio)
            end) t;
          match !best with
          | None -> Unbounded
          | Some (r, _) ->
              let leaving = basis.(r) in
              pivot r e;
              trace (k + 1) e leaving (Q.sub Q.zero z.(n + m));
              loop (k + 1)
    end
  in
  loop 0

Chvátal's first example.

(* Chvatal's example: maximize 5x1 + 4x2 + 3x3 subject to
     2x1 + 3x2 +  x3 <= 5
     4x1 +  x2 + 2x3 <= 11
     3x1 + 4x2 + 2x3 <= 8,   x >= 0. *)
let () =
  let q = Q.of_int in
  let a = [| [| q 2; q 3; q 1 |]; [| q 4; q 1; q 2 |]; [| q 3; q 4; q 2 |] |] and b = [| q 5; q 11; q 8 |] and c = [| q 5; q 4; q 3 |] in
  let name j = if j < 3 then Printf.sprintf "x%d" (j + 1) else Printf.sprintf "s%d" (j - 2) in
  let trace k e l v = Printf.printf "pivot %d: %s enters, %s leaves, objective %s\n" k (name e) (name l) (Q.to_string v) in
  match Simplex.solve ~rule:Simplex.Dantzig ~trace a b c with
  | Simplex.Optimal (v, x, y) ->
      Printf.printf "optimum %s at x = (%s)\n" (Q.to_string v) (String.concat ", " (Array.to_list (Array.map Q.to_string x)));
      Printf.printf "dual solution y = (%s), b.y = %s\n" (String.concat ", " (Array.to_list (Array.map Q.to_string y)))
        (Q.to_string (Array.fold_left Q.add Q.zero (Array.map2 Q.mul b y)))
  | _ -> print_endline "no optimum"

Running it.

pivot 1: x1 enters, s1 leaves, objective 25/2
pivot 2: x3 enters, s3 leaves, objective 13
optimum 13 at x = (2, 0, 1)
dual solution y = (1, 0, 1), b.y = 13

The method takes two pivots to reach the optimum 13 at $x = (2, 0, 1)$ , the solution given by Chvátal^[4]. The arithmetic is exact: floating-point implementations have to cope with rounding in the ratio test and in the elimination, and production solvers use careful tolerances and refactorization for this.

§ 03

Duality

Every linear program has a dual. For the program above it is^[14]:

minimize b^{T} y subject to A^{T} y \geq c, y \geq 0

Weak duality, $c^{T} x \leq y^{T} A x \leq b^{T} y$ for any feasible $x$ and $y$ , shows that any dual solution bounds the primal optimum. Strong duality says the two optima are equal when either exists. The simplex method proves it constructively: at the optimal tableau, the negated reduced costs of the slack variables form a feasible dual solution with the same value. In the example above $y = (1, 0, 1)$ and $b^{T} y = 13$ , a certificate that no feasible $x$ does better, checkable without repeating the computation. Complementary slackness says that at the optimum, $y_{i} > 0$ only for constraints that are tight: the second constraint has slack, and its dual variable is 0.

§ 04

Degeneracy and cycling

A vertex is degenerate when more than $n$ constraints are tight at it, so that some basic variable is zero. A pivot at such a vertex can change the basis without moving, with a ratio of zero and no improvement in the objective. The method can then return to a basis it has already visited, and cycle forever. Beale gave a small example in 1955^[7]:

Beale's example under Dantzig's rule and Bland's rule.

(* Beale's example, degenerate at the origin:
     maximize 3/4 x1 - 20 x2 + 1/2 x3 - 6 x4
     subject to 1/4 x1 -  8 x2 -     x3 + 9 x4 <= 0
                1/2 x1 - 12 x2 - 1/2 x3 + 3 x4 <= 0
                                     x3        <= 1 *)
let () =
  let q = Q.of_int and f a b = Q.make a b in
  let a = [| [| f 1 4; q (-8); q (-1); q 9 |]; [| f 1 2; q (-12); f (-1) 2; q 3 |]; [| q 0; q 0; q 1; q 0 |] |] in
  let b = [| q 0; q 0; q 1 |] and c = [| f 3 4; q (-20); f 1 2; q (-6) |] in
  let name j = if j < 4 then Printf.sprintf "x%d" (j + 1) else Printf.sprintf "s%d" (j - 3) in
  List.iter
    (fun (label, rule) ->
      Printf.printf "%s:\n" label;
      let trace k e l v = Printf.printf "  pivot %d: %s enters, %s leaves, objective %s\n" k (name e) (name l) (Q.to_string v) in
      match Simplex.solve ~rule ~trace a b c with
      | Simplex.Optimal (v, _, _) -> Printf.printf "  optimum %s\n" (Q.to_string v)
      | Simplex.Cycled -> print_endline "  back at an earlier basis: cycling"
      | Simplex.Unbounded -> print_endline "  unbounded")
    [ ("Dantzig's rule", Simplex.Dantzig); ("Bland's rule", Simplex.Bland) ]

Running it.

Dantzig's rule:
  pivot 1: x1 enters, s1 leaves, objective 0
  pivot 2: x2 enters, s2 leaves, objective 0
  pivot 3: x3 enters, x1 leaves, objective 0
  pivot 4: x4 enters, x2 leaves, objective 0
  pivot 5: s1 enters, x3 leaves, objective 0
  pivot 6: s2 enters, x4 leaves, objective 0
  back at an earlier basis: cycling
Bland's rule:
  pivot 1: x1 enters, s1 leaves, objective 0
  pivot 2: x2 enters, s2 leaves, objective 0
  pivot 3: x3 enters, x1 leaves, objective 0
  pivot 4: x4 enters, x2 leaves, objective 0
  pivot 5: x1 enters, s3 leaves, objective 1/5
  pivot 6: s1 enters, x4 leaves, objective 5/4
  optimum 5/4

With the largest reduced cost entering and ties in the ratio test broken by row, the method makes six pivots at the origin without improving the objective and returns to the starting basis. Bland's rule, which chooses the entering variable with the smallest index among those with positive reduced cost and breaks ties in the ratio test by the smallest index, provably never cycles^[6], and here reaches the optimum 5/4. Lexicographic tie-breaking and perturbation of $b$ are the other standard remedies.

§ 05

Complexity

Klee and Minty constructed a deformed cube in $n$ dimensions on which Dantzig's rule visits all $2^{n}$ vertices^[5]:

maximize j = 1 \sum n 2^{n - j} x_{j} subject to j = 1 \sum i - 1 2^{i - j + 1} x_{j} + x_{i} \leq 5^{i} (i = 1, \dots, n), x \geq 0

Klee–Minty cubes of dimension 2 to 8.

(* The Klee-Minty cube in n dimensions:
     maximize  sum_j 2^(n-j) x_j
     subject to  sum_(j<i) 2^(i-j+1) x_j + x_i <= 5^i   (i = 1..n),  x >= 0.
   With Dantzig's largest-coefficient rule the simplex method visits all
   2^n vertices. *)
let klee_minty n =
  let q = Q.of_int in
  let pow b e = let r = ref 1 in for _ = 1 to e do r := !r * b done; !r in
  let a = Array.init n (fun i -> Array.init n (fun j -> if j < i then q (pow 2 (i - j + 1)) else if j = i then Q.one else Q.zero)) in
  let b = Array.init n (fun i -> q (pow 5 (i + 1))) and c = Array.init n (fun j -> q (pow 2 (n - j - 1))) in
  (a, b, c)

Running it.

n   pivots (Dantzig) pivots (Bland)    optimum
2                  3              3         25
3                  7              5        125
4                 15              9        625
5                 31             15       3125
6                 63             25      15625
7                127             41      78125
8                255             67     390625

Dantzig's rule takes $2^{n} - 1$ pivots, exactly. Bland's rule does better on this family but is exponential on others, and exponential examples are known for most deterministic pivoting rules. Whether some pivoting rule is polynomial is open, and is related to the polynomial Hirsch conjecture on the diameter of polytopes; Santos disproved the original Hirsch conjecture in 2012^[12]. In practice the number of pivots is usually a small multiple of the number of constraints. Borgwardt proved polynomial average-case behaviour for a probabilistic model in 1982^[11], and Spielman and Teng's smoothed analysis showed that the expected number of pivots is polynomial when the input is slightly perturbed at random, which explains why bad cases are not met in practice^[10].

Linear programming itself is solvable in polynomial time. Khachiyan's ellipsoid method proved it in 1979^[8], though it is slow in practice, and Karmarkar's interior-point method of 1984 was both polynomial and practical^[9]. Modern solvers offer both simplex and interior-point methods; the simplex method remains preferred when a sequence of related problems is solved, because it can restart from the previous optimal basis.

§ 06

Variants

The revised simplex method keeps a factorization of the basis matrix instead of the whole tableau, which is far cheaper for large sparse problems. The dual simplex method maintains dual feasibility and works toward primal feasibility, and is the method of choice after adding a constraint, as in branch-and-bound for integer programming. When $x = 0$ is infeasible, a first phase solves an auxiliary program to find a feasible basis, the two-phase method^[4].

In SMT solvers, linear real arithmetic is decided by a general simplex method in the form of Dutertre and de Moura, which works with bounds on variables instead of a fixed objective, checks feasibility incrementally as the SAT solver asserts and retracts bounds, and produces a conflict explanation from the row that cannot be satisfied^[13]; see theory solver and DPLL(T). Difference logic, a special case, has a faster graph-based decision procedure.

§ 07

History

Kantorovich formulated linear programming problems of production planning in 1939^[3]. Dantzig devised the simplex method in 1947 while working on planning problems for the United States Air Force, and published it in 1951^[1]; his 1963 book became the standard reference^[2]. The duality theorem was developed with von Neumann and proved by Gale, Kuhn and Tucker in 1951^[14]. Beale showed that cycling can occur in 1955^[7], Klee and Minty that the method can take exponential time in 1972^[5], and Bland gave his anti-cycling rule in 1977^[6]. Khachiyan and Karmarkar gave polynomial algorithms in 1979 and 1984^[8]^[9], and Spielman and Teng explained the method's practical efficiency in 2004^[10].