RationalFeasibility.lean

1 import KomorebiLSM.PositiveLift
2 import Mathlib.Data.Rat.Cast.Order
3 import Mathlib.Data.Finset.Max
4 
5 /-! Exact rational feasibility, including boundary faces. Fourier–Motzkin
6 elimination transfers a finite weak linear system from the reals to rationals. -/
7 set_option autoImplicit false
8 namespace Komorebi.LSM.RationalSearch
9 open scoped BigOperators
10 
11 /-- A finite collection of rational lower and upper bounds that contains a real
12 point contains a rational point. Endpoints need not have positive separation. -/
13 theorem rational_interval {I : Type*} [Fintype I]
14     (lo hi : I → ℚ) (L U : I → Prop) [DecidablePred L] [DecidablePred U]
15     (t : ℝ) (hl : ∀ i, L i → (lo i : ℝ) ≤ t)
16     (hu : ∀ i, U i → t ≤ (hi i : ℝ)) :
17     ∃ q : ℚ, (∀ i, L i → lo i ≤ q) ∧ (∀ i, U i → q ≤ hi i) := by
18   classical
19   let ls := Finset.univ.filter L
20   let us := Finset.univ.filter U
21   by_cases hn : ls.Nonempty
22   · obtain ⟨i, hiL, hmax⟩ := ls.exists_max_image lo hn
23     have hLi : L i := (Finset.mem_filter.mp hiL).2
24     refine ⟨lo i, ?_, ?_⟩
25     · intro j hj
26       exact hmax j (Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩)
27     · intro j hj
28       exact_mod_cast (hl i hLi).trans (hu j hj)
29   · by_cases hm : us.Nonempty
30     · obtain ⟨i, hiU, hmin⟩ := us.exists_min_image hi hm
31       refine ⟨hi i, ?_, ?_⟩
32       · intro j hj
33         exact False.elim (hn ⟨j,Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩⟩)
34       · intro j hj
35         exact hmin j (Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩)
36     · refine ⟨0, ?_, ?_⟩
37       · intro j hj
38         exact False.elim (hn ⟨j,Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩⟩)
39       · intro j hj
40         exact False.elim (hm ⟨j,Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩⟩)
41 
42 /-- Every feasible finite weak linear system with rational coefficients has a
43 rational solution. Equality and nonnegativity are encoded as weak inequalities. -/
44 theorem rational_feasible (n : ℕ) : ∀ {I : Type*} [Fintype I]
45     (a : I → Fin n → ℚ) (b : I → ℚ),
46     (∃ x : Fin n → ℝ, ∀ i, ∑ j, (a i j : ℝ) * x j ≤ (b i : ℝ)) →
47     ∃ x : Fin n → ℚ, ∀ i, ∑ j, a i j * x j ≤ b i := by
48   induction n with
49   | zero =>
50     intro I _ a b ⟨x,hx⟩
51     refine ⟨fun _ => 0, ?_⟩
52     intro i
53     simpa using (show (0 : ℚ) ≤ b i by exact_mod_cast (by simpa using hx i : (0 : ℝ) ≤ b i))
54   | succ n ih =>
55     intro I _ a b ⟨x,hx⟩
56     classical
57     let c : I → ℚ := fun i => a i 0
58     let d : I → Fin n → ℚ := fun i j => a i j.succ
59     let aa : I ⊕ (I × I) → Fin n → ℚ
60       | .inl i => if c i = 0 then d i else 0
61       | .inr (i,k) => if c i < 0 ∧ 0 < c k then
62           fun j => d i j / (-c i) + d k j / c k else 0
63     let bb : I ⊕ (I × I) → ℚ
64       | .inl i => if c i = 0 then b i else 0
65       | .inr (i,k) => if c i < 0 ∧ 0 < c k then b i / (-c i) + b k / c k else 0
66     have hx' (i : I) : (c i : ℝ) * x 0 + ∑ j, (d i j : ℝ) * x j.succ ≤ b i := by
67       simpa only [c,d,Fin.sum_univ_succ] using hx i
68     have hred : ∃ y : Fin n → ℝ, ∀ i, ∑ j, (aa i j : ℝ) * y j ≤ bb i := by
69       refine ⟨fun j => x j.succ, ?_⟩
70       intro row
71       cases row with
72       | inl i =>
73         by_cases h : c i = 0
74         · simpa [aa,bb,h] using hx' i
75         · simp [aa,bb,h]
76       | inr p =>
77         obtain ⟨i,k⟩ := p
78         by_cases h : c i < 0 ∧ 0 < c k
79         · have hi : (0 : ℝ) < -(c i : ℝ) := by exact_mod_cast neg_pos.mpr h.1
80           have hk : (0 : ℝ) < (c k : ℝ) := by exact_mod_cast h.2
81           have hh := add_le_add (div_le_div_of_nonneg_right (hx' i) hi.le)
82             (div_le_div_of_nonneg_right (hx' k) hk.le)
83           simp only [aa,bb,ite_eq_left h,Rat.cast_add,Rat.cast_div,Rat.cast_neg]
84           simp_rw [add_mul,div_mul_eq_mul_div,Finset.sum_add_distrib,← Finset.sum_div]
85           have ci : (c i : ℝ) ≠ 0 := ne_of_lt (by linarith)
86           have ck : (c k : ℝ) ≠ 0 := ne_of_gt hk
87           convert hh using 1 <;> first | rfl | (field_simp; ring)
88         · simp [aa,bb,h]
89     obtain ⟨y,hy⟩ := ih aa bb hred
90     let r : I → ℚ := fun i => b i - ∑ j, d i j * y j
91     have hz (i : I) (h : c i = 0) : 0 ≤ r i := by
92       have := hy (.inl i)
93       simp only [aa,bb,ite_eq_left h] at this
94       exact sub_nonneg.mpr this
95     have hp (i k : I) (hi : c i < 0) (hk : 0 < c k) : r i / c i ≤ r k / c k := by
96       have hh := hy (.inr (i,k))
97       simp only [aa,bb,ite_eq_left (show c i < 0 ∧ 0 < c k from ⟨hi,hk⟩)] at hh
98       simp_rw [add_mul,div_mul_eq_mul_div,Finset.sum_add_distrib,← Finset.sum_div] at hh
99       dsimp [r]
100       simp only [sub_div,div_neg] at hh ⊢
101       linarith
102     -- The eliminated rational system makes every lower bound ≤ every upper bound.
103     have hpoint : ∃ t : ℝ, (∀ i, c i < 0 → (r i / c i : ℚ) ≤ t) ∧
104         (∀ i, 0 < c i → t ≤ (r i / c i : ℚ)) := by
105       let ls := Finset.univ.filter (fun i => c i < 0)
106       by_cases hn : ls.Nonempty
107       · obtain ⟨i,hi,hmax⟩ := ls.exists_max_image (fun i => r i / c i) hn
108         refine ⟨(r i / c i : ℚ),?_,?_⟩
109         · intro j hj
110           exact_mod_cast hmax j (Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩)
111         · intro j hj
112           exact_mod_cast hp i j (Finset.mem_filter.mp hi).2 hj
113       · let us := Finset.univ.filter (fun i => 0 < c i)
114         by_cases hm : us.Nonempty
115         · obtain ⟨i,hi,hmin⟩ := us.exists_min_image (fun i => r i / c i) hm
116           refine ⟨(r i / c i : ℚ),?_,?_⟩
117           · intro j hj
118             exact False.elim (hn ⟨j,Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩⟩)
119           · intro j hj
120             exact_mod_cast hmin j (Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩)
121         · refine ⟨0,?_,?_⟩
122           · intro j hj
123             exact False.elim (hn ⟨j,Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩⟩)
124           · intro j hj
125             exact False.elim (hm ⟨j,Finset.mem_filter.mpr ⟨Finset.mem_univ _,hj⟩⟩)
126     obtain ⟨t,htl,htu⟩ := hpoint
127     obtain ⟨q,hql,hqu⟩ := rational_interval (fun i => r i / c i) (fun i => r i / c i)
128       (fun i => c i < 0) (fun i => 0 < c i) t htl htu
129     refine ⟨Fin.cons q y, ?_⟩
130     intro i
131     rw [Fin.sum_univ_succ]
132     change c i * q + ∑ j, d i j * y j ≤ b i
133     have hr : c i * q ≤ r i := by
134       rcases lt_trichotomy (c i) 0 with hi | hi | hi
135       · exact (div_le_iff_of_neg' hi).mp (hql i hi)
136       · simpa [hi] using hz i hi
137       · simpa only [mul_comm] using (le_div_iff₀ hi).mp (hqu i hi)
138     dsimp [r] at hr
139     linarith
140 /-- Nonnegative real solutions of rational equations have nonnegative rational
141 solutions. This includes solutions on faces with zero coordinates. -/
142 theorem rational_nonnegative_solution {I : Type*} [Fintype I] (m : ℕ)
143     (V : I → Fin m → ℚ) (z : I → ℚ)
144     (x : Fin m → ℝ) (hx : ∀ j, 0 ≤ x j)
145     (he : ∀ i, ∑ j, (V i j : ℝ) * x j = (z i : ℝ)) :
146     ∃ q : Fin m → ℚ, (∀ j, 0 ≤ q j) ∧ (∀ i, ∑ j, V i j * q j = z i) := by
147   classical
148   let a : (I ⊕ I) ⊕ Fin m → Fin m → ℚ
149     | .inl (.inl i) => V i
150     | .inl (.inr i) => -V i
151     | .inr j => fun k => if k = j then -1 else 0
152   let b : (I ⊕ I) ⊕ Fin m → ℚ
153     | .inl (.inl i) => z i
154     | .inl (.inr i) => -z i
155     | .inr _ => 0
156   have hr : ∃ y : Fin m → ℝ, ∀ i, ∑ j, (a i j : ℝ) * y j ≤ (b i : ℝ) := by
157     refine ⟨x, ?_⟩
158     intro i
159     rcases i with (i | i) | j
160     · exact (he i).le
161     · simpa [a,b,Finset.sum_neg_distrib] using neg_le_neg (he i).ge
162     · simpa [a,b,apply_ite,ite_mul] using neg_nonpos.mpr (hx j)
163   obtain ⟨q,hq⟩ := rational_feasible m a b hr
164   refine ⟨q,?_,?_⟩
165   · intro j
166     have := hq (.inr j)
167     simpa [a,b] using this
168   · intro i
169     have h₁ := hq (.inl (.inl i))
170     have h₂ := hq (.inl (.inr i))
171     simp only [a,b,Pi.neg_apply,neg_mul,Finset.sum_neg_distrib,neg_le_neg_iff] at h₁ h₂
172     exact le_antisymm h₁ h₂
173 
174 /-- Rational equations with a real solution have a rational solution; no signs
175 are imposed in this variant, used for the right inverse. -/
176 theorem rational_solution {I : Type*} [Fintype I] (m : ℕ)
177     (V : I → Fin m → ℚ) (z : I → ℚ) (x : Fin m → ℝ)
178     (he : ∀ i, ∑ j, (V i j : ℝ) * x j = (z i : ℝ)) :
179     ∃ q : Fin m → ℚ, ∀ i, ∑ j, V i j * q j = z i := by
180   let a : I ⊕ I → Fin m → ℚ | .inl i => V i | .inr i => -V i
181   let b : I ⊕ I → ℚ | .inl i => z i | .inr i => -z i
182   have hr : ∃ y : Fin m → ℝ, ∀ i, ∑ j, (a i j : ℝ) * y j ≤ (b i : ℝ) := by
183     refine ⟨x, ?_⟩
184     intro i
185     cases i with
186     | inl i => exact (he i).le
187     | inr i => simpa [a,b,Finset.sum_neg_distrib] using neg_le_neg (he i).ge
188   obtain ⟨q,hq⟩ := rational_feasible m a b hr
189   refine ⟨q,?_⟩
190   intro i
191   have h₁ := hq (.inl i)
192   have h₂ := hq (.inr i)
193   simp only [a,b,Pi.neg_apply,neg_mul,Finset.sum_neg_distrib,neg_le_neg_iff] at h₁ h₂
194   exact le_antisymm h₁ h₂
195 
196 end Komorebi.LSM.RationalSearch

Back to source controls