More work towards proving farkas lemma

PolyhedralCombinatorics.lean

@@ -5,12 +5,14 @@ import Utils.Finset
 
 import PolyhedralCombinatorics.LinearSystem.Basic
 import PolyhedralCombinatorics.LinearSystem.LinearConstraints
+import PolyhedralCombinatorics.LinearSystem.Notation
 
 import PolyhedralCombinatorics.Polyhedron.Basic
-import PolyhedralCombinatorics.Polyhedron.Duality
 import PolyhedralCombinatorics.Polyhedron.Notation
 
 import PolyhedralCombinatorics.Projection.Basic
 import PolyhedralCombinatorics.Projection.SemiSpace
 import PolyhedralCombinatorics.Projection.Computable
 import PolyhedralCombinatorics.Projection.FourierMotzkin
+
+import PolyhedralCombinatorics.Duality.Farkas

PolyhedralCombinatorics/Polyhedron/Duality.lean → PolyhedralCombinatorics/Duality/Farkas.lean

@@ -1,20 +1,19 @@
 import PolyhedralCombinatorics.Projection.FourierMotzkin
+
 import Mathlib.LinearAlgebra.Matrix.DotProduct
 
-namespace LinearSystem.Duality
-open Polyhedron Matrix
+namespace LinearSystem
+open Matrix
 
 variable {𝔽 : Type*} [LinearOrderedField 𝔽] {m n : ℕ}
 
-theorem farkas_lemma {A : Matrix (Fin m) (Fin n) 𝔽} {b : Fin m → 𝔽}
-  : (∃ x, x ∈ P(A, b)) ↔ ∀ y, y ≥ 0 → y ᵥ* A = 0 → y ⬝ᵥ b ≥ 0 := by
+theorem farkas_lemma {S : LinearSystem 𝔽 n}
+  : (∃ x, x ∈ S.solutions) ↔ ∀ y ≥ 0, y ᵥ* S.mat = 0 → y ⬝ᵥ S.vec ≥ 0 := by
   constructor <;> intro h
   . intro y y_nonneg hy
     obtain ⟨x, hx⟩ := h
-    rw [mem_ofLinearSystem_of] at hx
     have := dotProduct_le_dotProduct_of_nonneg_left hx y_nonneg
     rw [dotProduct_mulVec, hy, zero_dotProduct] at this
     assumption
   . by_contra hc
-    simp_rw [not_exists, ←Polyhedron.eq_empty_iff] at hc
     sorry

PolyhedralCombinatorics/LinearSystem/LinearConstraints.lean

@@ -92,7 +92,7 @@ theorem mem_ofConstraints_nil_solutions : x ∈ (@ofConstraints 𝔽 _ n []).sol
     cons y b (cons (-y) (-b) $ ofConstraints cs) := by
   rfl
 
-@[simp] theorem mem_solutions_ofConstraints (x : Fin n → 𝔽)
+@[simp] theorem mem_ofConstraints (x : Fin n → 𝔽)
   : x ∈ (ofConstraints cs).solutions ↔ ∀ c ∈ cs, c.valid x := by
   induction cs
   case nil =>
@@ -104,7 +104,7 @@ theorem mem_ofConstraints_nil_solutions : x ∈ (@ofConstraints 𝔽 _ n []).sol
   all constraints. -/
 @[simp] theorem solutions_ofConstraints
   : (ofConstraints cs).solutions = { x : Fin n → 𝔽 | ∀ c ∈ cs, c.valid x } :=
-  Set.ext_iff.mpr mem_solutions_ofConstraints
+  Set.ext_iff.mpr mem_ofConstraints
 
 end toSet_ofConstraints
 

PolyhedralCombinatorics/LinearSystem/Notation.lean

@@ -0,0 +1,40 @@
+import PolyhedralCombinatorics.LinearSystem.LinearConstraints
+
+namespace LinearSystem
+open Matrix
+
+variable {𝔽} [LinearOrderedField 𝔽] {n : ℕ}
+
+open Lean.Parser Lean.Elab.Term
+
+/-- Syntax for addressing variables in linear constraints.
+
+  `x_[n]` is shorthand for `Pi.single n 1`. -/
+notation "x_[" n "]" => Pi.single (n : Fin _) 1
+
+/-- Syntax category for linear constraints used to define polyhedra. -/
+declare_syntax_cat linConstraint
+syntax:60 term:61 " ≤ " term : linConstraint
+syntax:60 term:61 " <= " term : linConstraint
+syntax:60 term:61 " = " term : linConstraint
+syntax:60 term:61 " ≥ " term : linConstraint
+syntax:60 term:61 " >= " term : linConstraint
+
+/-- Syntax for declaring polyhedra as systems of linear constraints. -/
+syntax:max (name := linSystem)
+  "!L" ("[" term:90 "^" term "]")? "{" linConstraint,*,? "}" : term
+
+macro_rules
+  | `(!L[$t^$n]{$[$constraints],*}) => `((!L{$constraints,*} : LinearSystem $t $n))
+  | `(!L{$[$constraints],*}) => do
+    let constraints ← constraints.mapM (fun
+      | `(linConstraint| $x:term ≤ $y:term) => `(LinearConstraint.le $x $y)
+      | `(linConstraint| $x:term <= $y:term) => `(LinearConstraint.le $x $y)
+      | `(linConstraint| $x:term = $y:term) => `(LinearConstraint.eq $x $y)
+      | `(linConstraint| $x:term ≥ $y:term) => `(LinearConstraint.ge $x $y)
+      | `(linConstraint| $x:term >= $y:term) => `(LinearConstraint.ge $x $y)
+      | _ => Lean.Macro.throwUnsupported)
+    `(ofConstraints [$constraints,*])
+
+example := !L[ℚ^2]{2 • x_[1] ≤ 1}
+example : LinearSystem 𝔽 2 := !L{2 • x_[1] ≤ 1, x_[0] + x_[1] = 0}

PolyhedralCombinatorics/Polyhedron/Basic.lean

@@ -131,7 +131,7 @@ theorem toSet_convex : Convex 𝔽 p.toSet := Quotient.recOn p solutions_convex
 
 @[simp] theorem mem_ofConstraints {constraints : List (LinearConstraint 𝔽 n)}
   : ∀ x, x ∈ ofLinearSystem (ofConstraints constraints) ↔ ∀ c ∈ constraints, c.valid x :=
-  mem_solutions_ofConstraints
+  LinearSystem.mem_ofConstraints
 
 /-- The set of points of a polyhedron defined by a given list of `constraints` is the set of
   points that satisfy those constraints. -/

PolyhedralCombinatorics/Polyhedron/Notation.lean

@@ -1,4 +1,5 @@
 import PolyhedralCombinatorics.Polyhedron.Basic
+import PolyhedralCombinatorics.LinearSystem.Notation
 
 namespace Polyhedron
 open Matrix LinearSystem
@@ -10,34 +11,13 @@ open Lean.Parser Lean.Elab.Term
 @[inherit_doc ofLinearSystem]
 scoped notation:max (name := matVecPolyhedron) "P(" A " , " b ")" => ofLinearSystem $ of A b
 
-/-- Syntax for addressing variables in linear constraints.
-
-  `x_[n]` is shorthand for `Pi.single n 1`. -/
-notation "x_[" n "]" => Pi.single n 1
-
-/-- Syntax category for linear constraints used to define polyhedra. -/
-declare_syntax_cat linConstraint
-syntax:60 term:61 " ≤ " term : linConstraint
-syntax:60 term:61 " <= " term : linConstraint
-syntax:60 term:61 " = " term : linConstraint
-syntax:60 term:61 " ≥ " term : linConstraint
-syntax:60 term:61 " >= " term : linConstraint
-
 /-- Syntax for declaring polyhedra as systems of linear constraints. -/
 syntax:max (name := systemPolyhedron)
   "!P" ("[" term:90 "^" term "]")? "{" linConstraint,*,? "}" : term
 
 macro_rules
-  | `(!P[$t^$n]{$[$constraints],*}) => `((!P{$constraints,*} : Polyhedron $t $n))
-  | `(!P{$[$constraints],*}) => do
-    let constraints ← constraints.mapM (fun
-      | `(linConstraint| $x:term ≤ $y:term) => `(LinearConstraint.le $x $y)
-      | `(linConstraint| $x:term <= $y:term) => `(LinearConstraint.le $x $y)
-      | `(linConstraint| $x:term = $y:term) => `(LinearConstraint.eq $x $y)
-      | `(linConstraint| $x:term ≥ $y:term) => `(LinearConstraint.ge $x $y)
-      | `(linConstraint| $x:term >= $y:term) => `(LinearConstraint.ge $x $y)
-      | _ => Lean.Macro.throwUnsupported)
-    `(ofLinearSystem $ ofConstraints [$constraints,*])
+  | `(!P[$t^$n]{$[$constraints],*}) => `(ofLinearSystem !L[$t^$n]{$constraints,*})
+  | `(!P{$[$constraints],*}) => `(ofLinearSystem $ !L{$constraints,*})
 
 example := !P[ℚ^2]{2 • x_[1] ≤ 1}
 example : Polyhedron 𝔽 2 := !P{2 • x_[1] ≤ 1, x_[0] + x_[1] = 0}

PolyhedralCombinatorics/Projection/Computable.lean

@@ -1,5 +1,7 @@
 import PolyhedralCombinatorics.Projection.SemiSpace
+import PolyhedralCombinatorics.LinearSystem.Notation
 
+import Mathlib.Data.Matrix.Reflection
 import Mathlib.Data.Finset.Sort
 import Mathlib.Data.Sum.Order
 
@@ -75,6 +77,58 @@ def computeProjection (S : LinearSystem 𝔽 n) (c : Fin n → 𝔽) : LinearSys
   simp_rw [Set.ext_iff, mem_projection]
   apply mem_computeProjection
 
+theorem computeProjection_mat_ortho {S : LinearSystem 𝔽 n} {c : Fin n → 𝔽}
+  : (computeProjection S c).mat *ᵥ c = 0 := by
+  unfold computeProjection
+  lift_lets
+  extract_lets _ Z _ _ r p D _
+  funext i
+  simp_rw [mulVec, Pi.zero_apply, D, Matrix.of_apply]
+  eta_reduce
+  split
+  . rename_i s _
+    have := mem_filter_univ.mp s.property
+    assumption
+  . simp only [sub_dotProduct, smul_dotProduct, smul_eq_mul]
+    rw [mul_comm, sub_self]
+
+theorem computeProjection_mat_conic {S : LinearSystem 𝔽 n} {c : Fin n → 𝔽}
+  : ∃ U : Matrix _ _ 𝔽,
+    (∀ i, U i ≥ 0)
+    ∧ U * S.mat = (computeProjection S c).mat
+    ∧ U *ᵥ S.vec = (computeProjection S c).vec := by
+  unfold computeProjection
+  lift_lets
+  extract_lets _ _ _ _ r p D d
+  let U : Matrix (Fin r) (Fin S.m) 𝔽 :=
+    Matrix.of fun i ↦
+      match p i with
+      | .inl s => Pi.single s 1
+      | .inr (s, t) => Pi.single ↑s (S.mat t ⬝ᵥ c) - Pi.single ↑t (S.mat s ⬝ᵥ c)
+  exists U
+  constructor
+  . simp_rw [U, Pi.le_def, of_apply, Pi.zero_apply]
+    intro i j
+    rcases p i with s | ⟨s, t⟩ <;> simp only
+    . rw [Pi.single_apply]
+      split <;> simp only [zero_le_one, le_refl]
+    . simp_rw [Pi.sub_apply, Pi.single_apply]
+      have hs := (mem_filter_univ.mp s.prop).le
+      have ht := (mem_filter_univ.mp t.prop).le
+      split <;> split <;> simp_all
+  constructor
+  funext i j
+  simp_rw [mul_apply', U, D, of_apply]
+  rotate_left
+  funext i
+  simp_rw [U, d, mulVec, of_apply]
+  all_goals (
+    rcases p i with s | ⟨s, t⟩ <;> simp only
+    . simp only [single_dotProduct, one_mul]
+    . simp_rw [sub_dotProduct, single_dotProduct]
+      rfl
+  )
+
 end LinearSystem
 
 namespace Polyhedron

PolyhedralCombinatorics/Projection/FourierMotzkin.lean

@@ -25,40 +25,23 @@ theorem mem_elim0 {S : LinearSystem 𝔽 (n + 1)} {x : Fin n → 𝔽}
   intro i
   simp_rw [mulVec, Function.comp_apply, dotProduct_cons, mul_zero, zero_add]
 
-end LinearSystem
-
-namespace Polyhedron
-open Matrix
-
-def fourierMotzkin (P : Polyhedron 𝔽 n) (j : Fin n) :=
-  !P{x_[j] = 0} ∩ P.computeProjection x_[j]
+def fourierMotzkin (S : LinearSystem 𝔽 n) (j : Fin n) :=
+  !L{x_[j] = 0} ++ S.computeProjection x_[j]
 
-theorem mem_fourierMotzkin (P : Polyhedron 𝔽 n) (j : Fin n) :
-  x ∈ P.fourierMotzkin j ↔ x j = 0 ∧ ∃ (α : 𝔽), x + Pi.single j α ∈ P := by
+theorem mem_fourierMotzkin  {S : LinearSystem 𝔽 n} {j : Fin n} :
+  x ∈ (S.fourierMotzkin j).solutions ↔ x j = 0 ∧ ∃ (α : 𝔽), x + Pi.single j α ∈ S.solutions := by
   simp_rw [
-    fourierMotzkin, mem_inter, mem_computeProjection, mem_ofConstraints,
+    fourierMotzkin, mem_concat, mem_computeProjection, mem_ofConstraints,
     List.mem_singleton, forall_eq, LinearConstraint.valid,
-    single_dotProduct, one_mul, and_congr_right_iff,
+    single_dotProduct, one_mul, and_congr_right_iff, mem_projection,
     ←Pi.single_smul, smul_eq_mul, mul_one, implies_true
   ]
 
-def elim0 (P : Polyhedron 𝔽 (n + 1)) : Polyhedron 𝔽 n :=
-  Quotient.recOn P (Polyhedron.ofLinearSystem ∘ LinearSystem.elim0)
-    fun S S' heq ↦ by
-      ext x
-      simp_rw [LinearSystem.equiv_def, Set.ext_iff] at heq
-      simp only [eq_rec_constant, Function.comp_apply, mem_ofLinearSystem, LinearSystem.mem_elim0]
-      apply heq
+def fourierMotzkinElim0 (S : LinearSystem 𝔽 (n + 1)) : LinearSystem 𝔽 n :=
+  (S.fourierMotzkin 0).elim0
 
-theorem mem_elim0 {P : Polyhedron 𝔽 (n + 1)} {x : Fin n → 𝔽}
-  : x ∈ P.elim0 ↔ vecCons 0 x ∈ P :=
-  Quotient.recOn P (fun S ↦ S.mem_elim0) (fun _ _ _ ↦ rfl)
-
-def fourierMotzkinElim0 (P : Polyhedron 𝔽 (n + 1)) : Polyhedron 𝔽 n :=
-  (P.fourierMotzkin 0).elim0
-
-theorem mem_fourierMotzkinElim0 {P : Polyhedron 𝔽 (n + 1)} {x : Fin n → 𝔽}
-  : x ∈ P.fourierMotzkinElim0 ↔ ∃ x₀ : 𝔽, vecCons x₀ x ∈ P := by
+theorem mem_fourierMotzkinElim0 {S : LinearSystem 𝔽 (n + 1)} {x : Fin n → 𝔽}
+  : x ∈ S.fourierMotzkinElim0.solutions ↔ ∃ x₀ : 𝔽, vecCons x₀ x ∈ S.solutions := by
   unfold fourierMotzkinElim0
   suffices ∀ α, (Pi.single 0 α : Fin (n + 1) → 𝔽) = vecCons α 0 by
     simp_rw [mem_elim0, mem_fourierMotzkin, cons_val_zero, true_and, cons_add,
@@ -70,25 +53,27 @@ theorem mem_fourierMotzkinElim0 {P : Polyhedron 𝔽 (n + 1)} {x : Fin n → 
   . simp_rw [Pi.single_eq_of_ne (f := fun _ ↦ 𝔽) (Fin.succ_ne_zero _) _, cons_val_succ,
       Pi.zero_apply]
 
-theorem fourierMotzkinElim0_eq_empty_iff {P : Polyhedron 𝔽 (n + 1)}
-  : P.fourierMotzkinElim0 = ∅ ↔ P = ∅ := by
-  simp_rw [eq_empty_iff, mem_fourierMotzkinElim0, not_exists]
+theorem fourierMotzkinElim0_eq_empty_iff {S : LinearSystem 𝔽 (n + 1)}
+  : S.fourierMotzkinElim0.solutions = ∅ ↔ S.solutions = ∅ := by
+  simp_rw [Set.eq_empty_iff_forall_not_mem, mem_fourierMotzkinElim0, not_exists]
   constructor <;> intro h x
   . rw [←Matrix.cons_head_tail x]
     apply h
   . intros
     apply h
 
-def fourierMotzkinElimRec {n : ℕ} (P : Polyhedron 𝔽 n) : Polyhedron 𝔽 0 :=
+def fourierMotzkinElimRec {n : ℕ} (S : LinearSystem 𝔽 n) : LinearSystem 𝔽 0 :=
   match n with
-  | 0 => P
-  | _ + 1 => fourierMotzkinElimRec P.fourierMotzkinElim0
+  | 0 => S
+  | _ + 1 => fourierMotzkinElimRec S.fourierMotzkinElim0
 
-@[simp] theorem recElimDimensions_eq_empty_iff (P : Polyhedron 𝔽 n)
-  : P.fourierMotzkinElimRec = ∅ ↔ P = ∅ := by
+@[simp] theorem recElimDimensions_eq_empty_iff (S : LinearSystem 𝔽 n)
+  : S.fourierMotzkinElimRec.solutions = ∅ ↔ S.solutions = ∅ := by
   induction n
   case zero => rfl
   case succ n ih =>
     transitivity
     . apply ih
-    . exact P.fourierMotzkinElim0_eq_empty_iff
+    . exact S.fourierMotzkinElim0_eq_empty_iff
+
+end LinearSystem