Fourier–Motzkin elimination

Fourier–Motzkin elimination is the analogue for inequalities of Gaussian elimination for equations^[1]^[2]. To remove a variable $x$ from a system of linear inequalities, write each inequality that mentions $x$ as a lower bound $L \leq x$ or an upper bound $x \leq U$ , where $L$ and $U$ are linear in the other variables. A value of $x$ exists if and only if every lower bound is at most every upper bound, so the inequalities in $x$ can be replaced by all the inequalities $L \leq U$ . The solution set of the new system is the projection of the old one along $x$ ^[4].

§ 01

The elimination step

Scale each inequality so that $x$ has coefficient $+ 1$ or $- 1$ :

x \leq U_{i} (i = 1.. p) L_{j} \leq x (j = 1.. q)

and replace them by the $pq$ inequalities $L_{j} \leq U_{i}$ , keeping the inequalities that do not mention $x$ . Each new inequality is a non-negative combination of two old ones, so every solution of the old system satisfies the new one; conversely, given a solution of the new system, any $x$ between $max_{j} L_{j}$ and $min_{i} U_{i}$ extends it. Once every variable is eliminated, the system is a list of comparisons between constants, and it is satisfiable if and only if they all hold^[4].

Fourier–Motzkin elimination over the rationals in OCaml, using Zarith.

(* A constraint a1 x1 + ... + an xn <= b over the rationals, with the
   coefficients in an array indexed by variable. *)
type c = { a : Q.t array; b : Q.t }
 
(* Eliminate variable k: every lower bound on xk is combined with every
   upper bound, and constraints without xk are kept. *)
let eliminate k cs =
  let pos = List.filter (fun c -> Q.sign c.a.(k) > 0) cs
  and neg = List.filter (fun c -> Q.sign c.a.(k) < 0) cs
  and zero = List.filter (fun c -> Q.sign c.a.(k) = 0) cs in
  let comb p n =
    (* scale so that xk has coefficient +1 in p and -1 in n, then add *)
    let sp = Q.inv p.a.(k) and sn = Q.inv (Q.neg n.a.(k)) in
    { a = Array.mapi (fun i x -> Q.add (Q.mul sp x) (Q.mul sn n.a.(i))) p.a; b = Q.add (Q.mul sp p.b) (Q.mul sn n.b) }
  in
  (* drop constraints 0 <= b that hold trivially *)
  List.filter (fun c -> Array.exists (fun x -> Q.sign x <> 0) c.a || Q.sign c.b < 0)
    (zero @ List.concat_map (fun p -> List.map (comb p) neg) pos)
 
(* Bounds on xk once x0 .. x(k-1) are fixed, from the constraints that
   mention no later variable. *)
let bounds k xs cs =
  List.fold_left (fun (lo, hi) c ->
    let rest = ref c.b in
    Array.iteri (fun i x -> if i < k then rest := Q.sub !rest (Q.mul c.a.(i) x)) xs;
    let s = Q.sign c.a.(k) in
    if s > 0 then (lo, Q.min hi (Q.div !rest c.a.(k)))
    else if s < 0 then (Q.max lo (Q.div !rest c.a.(k)), hi)
    else (lo, hi))
    (Q.minus_inf, Q.inf) cs
 
let pp names c =
  let term first (x, v) =
    let mag = Q.abs x in
    let coef = if Q.equal mag Q.one then "" else Q.to_string mag ^ " " in
    (if Q.sign x < 0 then (if first then "-" else " - ") else if first then "" else " + ") ^ coef ^ v
  in
  let terms = List.filter (fun (x, _) -> Q.sign x <> 0) (Array.to_list (Array.mapi (fun i x -> (x, names.(i))) c.a)) in
  (if terms = [] then "0" else String.concat "" (List.mapi (fun i t -> term (i = 0) t) terms))
  ^ " <= " ^ Q.to_string c.b

Eliminating z and y from a system in three variables, then back-substituting.

open Fm
 
let names = [| "x"; "y"; "z" |]
let c l b = { a = Array.of_list (List.map Q.of_int l); b = Q.of_int b }
 
(* Is there a point with x + y + z <= 8, x - y >= -1, y + 2z >= 7,
   x >= 1, z <= 3, and x + z >= 5 ? *)
let system = [ c [ 1; 1; 1 ] 8; c [ -1; 1; 0 ] 1; c [ 0; -1; -2 ] (-7); c [ -1; 0; 0 ] (-1); c [ 0; 0; 1 ] 3; c [ -1; 0; -1 ] (-5) ]
 
let show title cs = Printf.printf "%s (%d)\n" title (List.length cs); List.iter (fun c -> print_endline ("  " ^ pp names c)) cs

Running it.

system (6)
  x + y + z <= 8
  -x + y <= 1
  -y - 2 z <= -7
  -x <= -1
  z <= 3
  -x - z <= -5
z eliminated (6)
  -x + y <= 1
  -x <= -1
  x + 1/2 y <= 9/2
  y <= 3
  -1/2 y <= -1/2
  -x <= -2
y eliminated (4)
  -x <= -1
  -x <= -2
  -x <= 0
  2 x <= 8
x in [2, 4], take 3
y in [1, 3], take 2
z in [5/2, 3], take 11/4
all six constraints hold: true

Back-substitution goes in the opposite order: $x$ is chosen within the bounds of the final system, then $y$ within the bounds the previous system gives once $x$ is fixed, and so on. Each interval is non-empty because of the elimination step.

§ 02

Growth

One step can replace $p + q$ inequalities by $pq$ , up to $m^{2} /4$ from $m$ . After $k$ steps the bound is doubly exponential in $k$ , and most of the new inequalities are redundant^[4]^[6].

Twelve random inequalities in four variables, eliminating three of them.

open Fm
 
(* m random constraints in n variables; eliminate them one by one. *)
let () =
  Random.init 4;
  let n = 4 and m = 12 in
  let cs = ref (List.init m (fun _ ->
    { a = Array.init n (fun _ -> Q.of_int (Random.int 11 - 5)); b = Q.of_int (Random.int 20) })) in
  Printf.printf "start: %d constraints\n" (List.length !cs);
  for k = n - 1 downto 1 do
    cs := eliminate k !cs;
    Printf.printf "after eliminating x%d: %d constraints\n" k (List.length !cs)
  done

Running it.

start: 12 constraints
after eliminating x3: 31 constraints
after eliminating x2: 234 constraints
after eliminating x1: 13412 constraints

Removing redundant inequalities after each step, with a linear programming check for each, keeps the growth singly exponential, since a projection of a polyhedron has at most exponentially many facets. In practice, satisfiability of large systems is decided with the simplex method instead; Fourier–Motzkin remains useful for projection, for small systems, and for producing explanations, since each derived inequality records which originals it combines^[6].

§ 03

Integers

Over the integers the step is unsound for satisfiability: the projection of the integer points is not the set of integer points of the projection.

A system with rational but no integer solutions.

open Fm
 
(* 1 <= 2x <= 1 and y = x: satisfiable over the rationals, not the integers. *)
let () =
  let names = [| "x"; "y" |] in
  let c l b = { a = Array.of_list (List.map Q.of_int l); b = Q.of_int b } in
  let cs = [ c [ 2; 0 ] 1; c [ -2; 0 ] (-1); c [ 1; -1 ] 0; c [ -1; 1 ] 0 ] in
  let s = eliminate 1 cs in
  List.iter (fun c -> print_endline (pp names c)) s;
  let lo, hi = bounds 0 [| Q.zero; Q.zero |] s in
  Printf.printf "x in [%s, %s]\n" (Q.to_string lo) (Q.to_string hi)

Running it.

2 x <= 1
-2 x <= -1
x in [1/2, 1/2]

Elimination leaves $x = 1/2$ , which is not an integer. Pugh's Omega test extends Fourier–Motzkin to integers by computing an exact shadow when coefficients allow, and otherwise a dark shadow and a finite set of splinters^[5]. SMT solvers for linear integer arithmetic use related projections alongside branch and bound^[6].

§ 04

History

Joseph Fourier described the method in 1826 as part of his work on systems of inequalities^[1]. Lloyd Dines rediscovered it in 1919^[3], and Theodore Motzkin rediscovered it in his 1936 dissertation, from which the double name comes^[2]. William Pugh built the Omega test on it in 1991 for dependence analysis in compilers^[5].