TensorBernsteinCertificate.lean

1 import Mathlib.RingTheory.Polynomial.Bernstein
2 import Mathlib.Algebra.MvPolynomial.Monad
3 import Mathlib.Analysis.SpecialFunctions.Bernstein
4 import Mathlib.Tactic
5 
6 /-! Finite exact rational tensor-Bernstein witnesses and complete binary subdivision trees. -/
7 set_option autoImplicit false
8 namespace Komorebi.LSM.TensorBernsteinCertificate
9 open scoped BigOperators
10 open MvPolynomial
11 noncomputable section
12 abbrev Polynomial₂ := MvPolynomial (Fin 2) ℚ
13 
14 noncomputable def value (p : Polynomial₂) (x : Fin 2 → ℝ) : ℝ :=
15   MvPolynomial.eval₂ (algebraMap ℚ ℝ) x p
16 
17 /-- All endpoints, widths, and submitted coefficients are rational. -/
18 structure Rectangle where
19   lo : Fin 2 → ℚ
20   hi : Fin 2 → ℚ
21   proper : ∀ i,lo i<hi i
22 
23 def Rectangle.width (R : Rectangle) (i : Fin 2) : ℚ := R.hi i-R.lo i
24 
25 def Rectangle.domain (R : Rectangle) : Set (Fin 2 → ℝ) :=
26   {x | ∀ i,x i∈Set.Icc (R.lo i : ℝ) (R.hi i : ℝ)}
27 
28 def unitSquare : Set (Fin 2 → ℝ) := {x | ∀ i,x i∈Set.Icc (0 : ℝ) 1}
29 
30 def Rectangle.localPolynomial (R : Rectangle) (p : Polynomial₂) : Polynomial₂ :=
31   MvPolynomial.bind₁ (fun i => C (R.lo i)+C (R.width i)*X i) p
32 
33 noncomputable def Rectangle.coordinate (R : Rectangle) (x : Fin 2 → ℝ) (i : Fin 2) : ℝ :=
34   (x i-R.lo i)/(R.width i : ℝ)
35 
36 theorem Rectangle.width_pos (R : Rectangle) (i : Fin 2) : 0<(R.width i : ℝ) := by
37   exact_mod_cast sub_pos.mpr (R.proper i)
38 
39 theorem Rectangle.coordinate_mem (R : Rectangle) {x : Fin 2 → ℝ} (hx : x∈R.domain) :
40     R.coordinate x∈unitSquare := by
41   intro i
42   have hi := hx i
43   have hw := R.width_pos i
44   constructor
45   · exact div_nonneg (sub_nonneg.mpr hi.1) hw.le
46   · apply (div_le_one hw).mpr
47     simpa [Rectangle.width] using sub_le_sub_right hi.2 (R.lo i : ℝ)
48 
49 theorem Rectangle.value_local (R : Rectangle) (p : Polynomial₂) (x : Fin 2 → ℝ) :
50     value (R.localPolynomial p) (R.coordinate x)=value p x := by
51   change MvPolynomial.eval₂Hom (algebraMap ℚ ℝ) (R.coordinate x)
52     (MvPolynomial.bind₁ (fun i => C (R.lo i)+C (R.width i)*X i) p)=_
53   rw [MvPolynomial.eval₂Hom_bind₁]
54   have he : (fun i => (MvPolynomial.eval₂Hom (algebraMap ℚ ℝ) (R.coordinate x))
55       (C (R.lo i)+C (R.width i)*X i))=x := by
56     funext i
57     simp [Rectangle.coordinate]
58     field_simp [(R.width_pos i).ne']
59     ring
60   rw [he]
61   rfl
62 
63 /-- The literal univariate Bernstein polynomial inserted in one tensor coordinate. -/
64 def basisPolynomial (n k : ℕ) (i : Fin 2) : Polynomial₂ :=
65   C (n.choose k : ℚ)*X i^k*(1-X i)^(n-k)
66 
67 noncomputable def basis (n k : ℕ) (x : ℝ) : ℝ :=
68   (n.choose k : ℝ)*x^k*(1-x)^(n-k)
69 
70 theorem basis_nonneg (n k : ℕ) {x : ℝ} (hx : x∈Set.Icc (0 : ℝ) 1) : 0≤basis n k x := by
71   unfold basis
72   have h0 := hx.1
73   have h1 : 0≤1-x := sub_nonneg.mpr hx.2
74   positivity
75 
76 theorem basis_sum (n : ℕ) (x : ℝ) : ∑ k : Fin (n+1),basis n k x=1 := by
77   have h := congrArg (Polynomial.eval x) (bernsteinPolynomial.sum ℝ n)
78   rw [Fin.sum_univ_eq_sum_range (fun k => basis n k x)]
79   simpa [basis,bernsteinPolynomial,Polynomial.eval_finsetSum] using h
80 
81 theorem value_basisPolynomial (n k : ℕ) (i : Fin 2) (x : Fin 2 → ℝ) :
82     value (basisPolynomial n k i) x=basis n k (x i) := by
83   simp [value,basisPolynomial,basis]
84 
85 def tensorPolynomial {m n : ℕ} (c : Fin (m+1) → Fin (n+1) → ℚ) : Polynomial₂ :=
86   ∑ i,∑ j,C (c i j)*basisPolynomial m i 0*basisPolynomial n j 1
87 
88 theorem value_tensorPolynomial {m n : ℕ} (c : Fin (m+1) → Fin (n+1) → ℚ)
89     (x : Fin 2 → ℝ) : value (tensorPolynomial c) x=
90       ∑ i,∑ j,(c i j : ℝ)*basis m i (x 0)*basis n j (x 1) := by
91   simp [tensorPolynomial,value,MvPolynomial.eval₂_sum,basisPolynomial,basis]
92 
93 /-- The identity is an exact equality of finitely supported rational polynomials. -/
94 structure LocalCertificate (p : Polynomial₂) (R : Rectangle) where
95   degreeX : ℕ
96   degreeY : ℕ
97   coefficient : Fin (degreeX+1) → Fin (degreeY+1) → ℚ
98   nonneg : ∀ i j,0≤coefficient i j
99   identity : R.localPolynomial p=tensorPolynomial coefficient
100 
101 theorem LocalCertificate.sound {p : Polynomial₂} {R : Rectangle} (cert : LocalCertificate p R)
102     {x : Fin 2 → ℝ} (hx : x∈R.domain) : 0≤value p x := by
103   rw [← R.value_local p x,cert.identity,value_tensorPolynomial]
104   have hu := R.coordinate_mem hx
105   apply Finset.sum_nonneg
106   intro i _
107   apply Finset.sum_nonneg
108   intro j _
109   exact mul_nonneg (mul_nonneg (by exact_mod_cast cert.nonneg i j)
110     (basis_nonneg _ _ (hu 0))) (basis_nonneg _ _ (hu 1))
111 
112 /-- A split is always a genuine midpoint split of the specified coordinate. -/
113 def Rectangle.lower (R : Rectangle) (axis : Fin 2) : Rectangle where
114   lo := R.lo
115   hi i := if i=axis then (R.lo i+R.hi i)/2 else R.hi i
116   proper i := by
117     split_ifs <;> linarith [R.proper i]
118 
119 def Rectangle.upper (R : Rectangle) (axis : Fin 2) : Rectangle where
120   lo i := if i=axis then (R.lo i+R.hi i)/2 else R.lo i
121   hi := R.hi
122   proper i := by
123     split_ifs <;> linarith [R.proper i]
124 
125 theorem Rectangle.split_cover (R : Rectangle) (axis : Fin 2) {x : Fin 2 → ℝ}
126     (hx : x∈R.domain) : x∈(R.lower axis).domain ∨ x∈(R.upper axis).domain := by
127   by_cases h : x axis≤((R.lo axis+R.hi axis)/2 : ℚ)
128   · left
129     intro i
130     constructor
131     · exact (hx i).1
132     · by_cases he : i=axis
133       · subst i
134         simpa [Rectangle.lower] using h
135       · simpa [Rectangle.lower,he] using (hx i).2
136   · right
137     intro i
138     constructor
139     · by_cases he : i=axis
140       · subst i
141         simpa [Rectangle.upper] using (le_of_lt (lt_of_not_ge h))
142       · simpa [Rectangle.upper,he] using (hx i).1
143     · exact (hx i).2
144 
145 /-- Every internal node contains both children: an arbitrary rectangle list is not accepted. -/
146 inductive SubdivisionWitness (p : Polynomial₂) : Rectangle → Type
147   | leaf {R} (certificate : LocalCertificate p R) : SubdivisionWitness p R
148   | split {R} (axis : Fin 2) (lower : SubdivisionWitness p (R.lower axis))
149       (upper : SubdivisionWitness p (R.upper axis)) : SubdivisionWitness p R
150 
151 theorem SubdivisionWitness.sound {p : Polynomial₂} {R : Rectangle} (w : SubdivisionWitness p R)
152     {x : Fin 2 → ℝ} (hx : x∈R.domain) : 0≤value p x := by
153   induction w with
154   | leaf cert => exact cert.sound hx
155   | @split R axis lower upper ihl ihu =>
156     rcases R.split_cover axis hx with hl|hu
157     · exact ihl hl
158     · exact ihu hu
159 
160 def square : Rectangle where
161   lo _ := -1
162   hi _ := 1
163   proper _ := by norm_num
164 
165 theorem square_domain (x : Fin 2 → ℝ) :
166     x∈square.domain ↔ ∀ i,x i∈Set.Icc (-1 : ℝ) 1 := by simp [Rectangle.domain,square]
167 
168 end
169 end Komorebi.LSM.TensorBernsteinCertificate

Back to source controls