EllipsoidPolynomial.lean

1 import KomorebiLSM.OrientedEllipsoid
2 import KomorebiLSM.ContinuousRecords
3 
4 /-! Finite polynomial/SOS witnesses for the oriented quadratic certificate. -/
5 set_option autoImplicit false
6 namespace Komorebi.LSM.OrientedEllipsoid
7 open scoped BigOperators Matrix Polynomial
8 open MeasureTheory
9 variable {D : Type*} [Fintype D] [DecidableEq D]
10 
11 /-- The literal matrix action as a continuous linear map, with no extra positivity field. -/
12 noncomputable def matrixOperator : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ →ₗ[ℝ]
13     (State D →L[ℝ] State D) :=
14   LinearMap.toContinuousLinearMap.toLinearMap.comp Matrix.toLin'.toLinearMap
15 
16 @[simp] theorem matrixOperator_apply (B : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ) (z : State D) :
17     matrixOperator B z = B *ᵥ z := rfl
18 
19 noncomputable def matrixPolynomial {n : ℕ}
20     (B : Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ) :
21     Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X] :=
22   fun i j => Polynomial.ofFn n (fun k => B k i j)
23 
24 omit [Fintype D] [DecidableEq D] in
25 theorem matrixPolynomial_eval {n : ℕ} (B : Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ)
26     (u : ℝ) : PolynomialGram.evalMatrix (matrixPolynomial B) u = ∑ k : Fin n, u^(k:ℕ) • B k := by
27   ext i j
28   simp [matrixPolynomial,PolynomialGram.evalMatrix,UniformIntervalMoments.eval_ofFn,
29     Matrix.sum_apply,Matrix.smul_apply,mul_comm]
30 
31 noncomputable def polynomialOperator {n : ℕ}
32     (B : Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ) : ℝ → State D →L[ℝ] State D :=
33   PolynomialGram.operatorPolynomial (fun k => matrixOperator (B k))
34 
35 theorem polynomialOperator_matrix {n : ℕ}
36     (B : Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ) (u : ℝ) :
37     polynomialOperator B u = matrixOperator (PolynomialGram.evalMatrix (matrixPolynomial B) u) := by
38   rw [matrixPolynomial_eval,map_sum]
39   simp only [map_smul]
40   rfl
41 
42 /-- A finite check of scalar coefficients, including the coefficients above neither degree. -/
43 theorem scalar_coefficient_identity (p q : ℝ[X]) (bound : ℕ)
44     (hp : p.natDegree≤bound) (hq : q.natDegree≤bound)
45     (hc : ∀ k : Fin (bound+1),p.coeff k=q.coeff k) : p=q := by
46   ext k
47   by_cases hk : k<bound+1
48   · exact hc ⟨k,hk⟩
49   · rw [Polynomial.coeff_eq_zero_of_natDegree_lt (lt_of_le_of_lt hp (by omega)),
50       Polynomial.coeff_eq_zero_of_natDegree_lt (lt_of_le_of_lt hq (by omega))]
51 
52 noncomputable def scalarSOS {n m : ℕ} (f : Fin n → ℝ[X]) (g : Fin m → ℝ[X]) : ℝ[X] :=
53   (∑ i,(f i)^2) + (1-Polynomial.X^2)*(∑ j,(g j)^2)
54 
55 theorem scalarSOS_nonneg {n m : ℕ} (f : Fin n → ℝ[X]) (g : Fin m → ℝ[X])
56     (u : ℝ) (hu : u∈Set.Icc (-1 : ℝ) 1) : 0≤(scalarSOS f g).eval u := by
57   have hfactor : 0≤1-u^2 := by nlinarith [hu.1,hu.2]
58   simp only [scalarSOS,Polynomial.eval_add,Polynomial.eval_finsetSum,Polynomial.eval_pow,
59     Polynomial.eval_mul,Polynomial.eval_sub,Polynomial.eval_one,Polynomial.eval_X]
60   exact add_nonneg (Finset.sum_nonneg (fun i _ => sq_nonneg _))
61     (mul_nonneg hfactor (Finset.sum_nonneg (fun j _ => sq_nonneg _)))
62 
63 /-- Scalar interval SOS identities are checked by finitely many actual coefficients. -/
64 structure ScalarCertificate (p : ℝ[X]) where
65   terms0 : ℕ
66   terms1 : ℕ
67   square0 : Fin terms0 → ℝ[X]
68   square1 : Fin terms1 → ℝ[X]
69   bound : ℕ
70   degree_left : p.natDegree≤bound
71   degree_right : (scalarSOS square0 square1).natDegree≤bound
72   coefficients : ∀ k : Fin (bound+1),p.coeff k=(scalarSOS square0 square1).coeff k
73 
74 theorem ScalarCertificate.nonneg {p : ℝ[X]} (cert : ScalarCertificate p)
75     (u : ℝ) (hu : u∈Set.Icc (-1 : ℝ) 1) : 0≤p.eval u := by
76   rw [scalar_coefficient_identity p _ cert.bound cert.degree_left cert.degree_right cert.coefficients]
77   exact scalarSOS_nonneg cert.square0 cert.square1 u hu
78 
79 noncomputable def betaPolynomial (c : D → ℝ)
80     (P : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) : ℝ[X] :=
81   P (Sum.inl ()) (Sum.inl ()) + ∑ i,Polynomial.C (c i)*P (Sum.inl ()) (Sum.inr i)
82 
83 noncomputable def dualPolynomial (H : Matrix D D ℝ)
84     (P : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) : ℝ[X] :=
85   ∑ i,∑ j,Polynomial.C (H⁻¹ i j)*P (Sum.inl ()) (Sum.inr i)*P (Sum.inl ()) (Sum.inr j)
86 
87 omit [DecidableEq D] in
88 theorem betaPolynomial_eval (c : D → ℝ)
89     (P : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) (u : ℝ) :
90     (betaPolynomial c P).eval u = beta c (PolynomialGram.evalMatrix P u) := by
91   simp [Polynomial.eval_finsetSum,betaPolynomial,beta,massConstant,massSpatial,PolynomialGram.evalMatrix,dotProduct,mul_comm]
92 
93 theorem dualPolynomial_eval (H : Matrix D D ℝ)
94     (P : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) (u : ℝ) :
95     (dualPolynomial H P).eval u = massSpatial (PolynomialGram.evalMatrix P u) ⬝ᵥ
96       (H⁻¹ *ᵥ massSpatial (PolynomialGram.evalMatrix P u)) := by
97   simp only [dualPolynomial,Polynomial.eval_finsetSum,Polynomial.eval_mul,Polynomial.eval_C,
98     dotProduct,Matrix.mulVec,massSpatial,PolynomialGram.evalMatrix,Finset.mul_sum]
99   apply Finset.sum_congr rfl
100   intro i _
101   apply Finset.sum_congr rfl
102   intro j _
103   ring
104 
105 noncomputable def quadraticPolynomial (c : D → ℝ) (H : Matrix D D ℝ)
106     (P : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) (tau : ℝ[X]) :
107     Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X] :=
108   P.transpose * PolynomialGram.constantMatrix (Q c H) * P -
109     tau • PolynomialGram.constantMatrix (Q c H)
110 
111 omit [Fintype D] [DecidableEq D] in
112 theorem evalMatrix_sub (P T : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) (u : ℝ) :
113     PolynomialGram.evalMatrix (P-T) u = PolynomialGram.evalMatrix P u - PolynomialGram.evalMatrix T u := by
114   ext i j
115   change (P i j-T i j).eval u = (P i j).eval u - (T i j).eval u
116   exact Polynomial.eval_sub (P i j) (T i j) u
117 
118 theorem quadraticPolynomial_eval (c : D → ℝ) (H : Matrix D D ℝ)
119     (P : Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ[X]) (tau : ℝ[X]) (u : ℝ) :
120     PolynomialGram.evalMatrix (quadraticPolynomial c H P tau) u =
121       (PolynomialGram.evalMatrix P u).transpose * Q c H * PolynomialGram.evalMatrix P u -
122         tau.eval u • Q c H := by
123   rw [quadraticPolynomial,evalMatrix_sub,PolynomialGram.evalMatrix_mul,
124     PolynomialGram.evalMatrix_mul,PolynomialGram.evalMatrix_transpose,
125     PolynomialGram.evalMatrix_constant,PolynomialGram.evalMatrix_scalar,
126     PolynomialGram.evalMatrix_constant]
127 
128 /-- Finite scalar SOS and matrix-Gram coefficient checks certify the complete interval. -/
129 structure PolynomialCertificate (c : D → ℝ) (H : Matrix D D ℝ) {A : Type*} {n : ℕ}
130     (B : A → Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ) (q bound : ℕ) where
131   tau : A → ℝ[X]
132   tauSOS : ∀ a,ScalarCertificate (tau a)
133   betaSOS : ∀ a,ScalarCertificate (betaPolynomial c (matrixPolynomial (B a)))
134   orientationSOS : ∀ a,ScalarCertificate
135     ((betaPolynomial c (matrixPolynomial (B a)))^2-dualPolynomial H (matrixPolynomial (B a)))
136   gram0 : A → Matrix (Fin (q+1) × (Unit ⊕ D)) (Fin (q+1) × (Unit ⊕ D)) ℝ
137   gram1 : A → Matrix (Fin q × (Unit ⊕ D)) (Fin q × (Unit ⊕ D)) ℝ
138   gram0_positive : ∀ a,(gram0 a).PosSemidef
139   gram1_positive : ∀ a,(gram1 a).PosSemidef
140   degree_left : ∀ a i j,(quadraticPolynomial c H (matrixPolynomial (B a)) (tau a) i j).natDegree≤bound
141   degree_right : ∀ a i j,(PolynomialGram.intervalGramPolynomial q (gram0 a) (gram1 a) i j).natDegree≤bound
142   coefficients : ∀ a i j (k : Fin (bound+1)),
143     (quadraticPolynomial c H (matrixPolynomial (B a)) (tau a) i j).coeff k =
144       (PolynomialGram.intervalGramPolynomial q (gram0 a) (gram1 a) i j).coeff k
145   mass_moments : ∀ a i,mass ((∑ k : Fin n,UniformIntervalMoments.moment k • matrixOperator (B a k))
146     (Pi.basisFun ℝ (Unit ⊕ D) i))=mass (Pi.basisFun ℝ (Unit ⊕ D) i)
147 
148 /-- A full polynomial-family quadratic witness constructs the literal P52 instrument. -/
149 noncomputable def polynomialDensityInstrument (c : D → ℝ) (H : Matrix D D ℝ)
150     (hH : H.PosDef) {A : Type*} {n : ℕ}
151     (B : A → Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ)
152     (branches : ∀ a u,u∈Set.Icc (-1 : ℝ) 1 →
153       BranchCertificate c H (PolynomialGram.evalMatrix (matrixPolynomial (B a)) u))
154     (hm : ∀ a z,mass ((∑ k : Fin n,UniformIntervalMoments.moment k • matrixOperator (B a k)) z)=mass z) :
155     ContinuousInstrument.DensityInstrument (A := A) (probabilityCone c H hH)
156       UniformIntervalMoments.measure where
157   op a := polynomialOperator (B a)
158   measurable_op a := (PolynomialGram.operatorPolynomial_continuous
159     (fun k => matrixOperator (B a k))).measurable
160   integrable_op a := PolynomialGram.operatorPolynomial_integrable (fun k => matrixOperator (B a k))
161   positive a := by
162     filter_upwards [PolynomialGram.uniform_ae_interval] with u hu
163     intro z hz
164     change polynomialOperator (B a) u z∈cone c H
165     rw [polynomialOperator_matrix,matrixOperator_apply]
166     exact BranchCertificate.preserves hH (branches a u hu) hz
167   mass_preserving a z := by
168     change mass ((∫ u,PolynomialGram.operatorPolynomial (fun k => matrixOperator (B a k)) u
169       ∂UniformIntervalMoments.measure) z)=mass z
170     rw [PolynomialGram.operatorPolynomial_integral]
171     exact hm a z
172 
173 @[simp] theorem polynomialDensityInstrument_op (c : D → ℝ) (H : Matrix D D ℝ)
174     (hH : H.PosDef) {A : Type*} {n : ℕ}
175     (B : A → Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ)
176     (branches : ∀ a u,u∈Set.Icc (-1 : ℝ) 1 →
177       BranchCertificate c H (PolynomialGram.evalMatrix (matrixPolynomial (B a)) u))
178     (hm : ∀ a z,mass ((∑ k : Fin n,UniformIntervalMoments.moment k • matrixOperator (B a k)) z)=mass z)
179     (a : A) (u : ℝ) :
180     (polynomialDensityInstrument c H hH B branches hm).op a u =
181       matrixOperator (PolynomialGram.evalMatrix (matrixPolynomial (B a)) u) :=
182   polynomialOperator_matrix (B a) u
183 
184 namespace PolynomialCertificate
185 variable {c : D → ℝ} {H : Matrix D D ℝ} {A : Type*} {n q bound : ℕ}
186   {B : A → Fin n → Matrix (Unit ⊕ D) (Unit ⊕ D) ℝ}
187   (cert : PolynomialCertificate c H B q bound)
188 include cert
189 
190 /-- Finite SOS/coefficient data imply every hypothesis at every interval point. -/
191 noncomputable def branch (a : A) (u : ℝ) (hu : u∈Set.Icc (-1 : ℝ) 1) :
192     BranchCertificate c H (PolynomialGram.evalMatrix (matrixPolynomial (B a)) u) where
193   tau := (cert.tau a).eval u
194   remainder := PolynomialGram.intervalGram q (cert.gram0 a) (cert.gram1 a) u
195   tau_nonneg := (cert.tauSOS a).nonneg u hu
196   remainder_positive := PolynomialGram.interval_gram_positive q (cert.gram0 a) (cert.gram1 a)
197     (cert.gram0_positive a) (cert.gram1_positive a) u hu
198   identity := by
199     have he := (PolynomialGram.polynomial_identity_finite
200       (quadraticPolynomial c H (matrixPolynomial (B a)) (cert.tau a))
201       (PolynomialGram.intervalGramPolynomial q (cert.gram0 a) (cert.gram1 a)) bound
202       (cert.degree_left a) (cert.degree_right a)).mpr (cert.coefficients a) u
203     rwa [quadraticPolynomial_eval,PolynomialGram.intervalGramPolynomial_eval] at he
204   beta_nonneg := by
205     simpa only [betaPolynomial_eval] using (cert.betaSOS a).nonneg u hu
206   orientation := by
207     simpa only [Polynomial.eval_sub,Polynomial.eval_pow,betaPolynomial_eval,dualPolynomial_eval]
208       using (cert.orientationSOS a).nonneg u hu
209 
210 /-- The finite basis moment checks imply the exact integrated mass functional identity. -/
211 theorem mass_preserving (a : A) (z : State D) :
212     mass ((∑ k : Fin n,UniformIntervalMoments.moment k • matrixOperator (B a k)) z)=mass z :=
213   (PolynomialGram.moment_mass_finite (Pi.basisFun ℝ (Unit ⊕ D)) mass
214     (fun k => matrixOperator (B a k))).mpr (cert.mass_moments a) z
215 
216 noncomputable def instrument (hH : H.PosDef) :
217     ContinuousInstrument.DensityInstrument (A := A) (probabilityCone c H hH)
218       UniformIntervalMoments.measure :=
219   polynomialDensityInstrument c H hH B cert.branch cert.mass_preserving
220 
221 @[simp] theorem instrument_op (hH : H.PosDef) (a : A) (u : ℝ) :
222     (cert.instrument hH).op a u = matrixOperator
223       (PolynomialGram.evalMatrix (matrixPolynomial (B a)) u) :=
224   polynomialDensityInstrument_op c H hH B cert.branch cert.mass_preserving a u
225 
226 /-- Positivity is pointwise on the entire closed interval, including its endpoints. -/
227 theorem preserves (hH : H.PosDef) (a : A) (u : ℝ) (hu : u∈Set.Icc (-1 : ℝ) 1)
228     {z : State D} (hz : z∈cone c H) : polynomialOperator (B a) u z∈cone c H := by
229   rw [polynomialOperator_matrix,matrixOperator_apply]
230   exact BranchCertificate.preserves hH (cert.branch a u hu) hz
231 
232 /-- Every literal finite product of interval observation operators stays in the actual cone. -/
233 theorem word_cone (hH : H.PosDef) (word : List (A × ℝ))
234     (hword : ∀ s∈word,s.2∈Set.Icc (-1 : ℝ) 1) {z : State D} (hz : z∈cone c H) :
235     FiniteInstruments.wordAction (fun s : A × ℝ => (polynomialOperator (B s.1) s.2).toLinearMap)
236       word z∈cone c H := by
237   induction word generalizing z with
238   | nil => exact hz
239   | cons s word ih =>
240     apply ih (fun t ht => hword t (List.mem_cons_of_mem s ht))
241     exact cert.preserves hH s.1 s.2 (hword s (List.mem_cons_self)) hz
242 
243 /-- Finite-word densities are nonnegative, without incorrectly bounding point densities by one. -/
244 theorem word_mass_nonneg (hH : H.PosDef) (word : List (A × ℝ))
245     (hword : ∀ s∈word,s.2∈Set.Icc (-1 : ℝ) 1) {z : State D} (hz : z∈cone c H) :
246     0≤mass (FiniteInstruments.wordAction
247       (fun s : A × ℝ => (polynomialOperator (B s.1) s.2).toLinearMap) word z) :=
248   (probabilityCone c H hH).mass_nonneg (cert.word_cone hH word hword hz)
249 
250 /-- The complete finite observation measures have total mass one at every horizon. -/
251 theorem observation_law_probability (hH : H.PosDef) (actions : ℕ → A)
252     (z : State D) (hz : z∈cone c H) (hm : mass z=1) (horizon : ℕ) :
253     IsProbabilityMeasure (ContinuousRecords.observationLaw (cert.instrument hH) actions z horizon) :=
254   ContinuousRecords.observationLaw_probability (cert.instrument hH) actions ⟨z,hz,hm⟩ horizon
255 
256 /-- Arbitrary common measurable randomized history policies give genuine normalized joint laws. -/
257 theorem record_law_probability [Fintype A] [MeasurableSpace A] [MeasurableSingletonClass A]
258     (hH : H.PosDef) (policy : ContinuousRecords.Policy (A := A) (Y := ℝ))
259     (z : State D) (hz : z∈cone c H) (hm : mass z=1) (horizon : ℕ) :
260     IsProbabilityMeasure (ContinuousRecords.recordLaw (cert.instrument hH) policy z horizon) :=
261   ContinuousRecords.recordLaw_probability (cert.instrument hH) policy ⟨z,hz,hm⟩ horizon
262 
263 /-- The normalized laws are the literal prefix marginals of every longer finite law. -/
264 theorem record_law_prefix [Fintype A] [MeasurableSpace A] [MeasurableSingletonClass A]
265     (hH : H.PosDef) (policy : ContinuousRecords.Policy (A := A) (Y := ℝ))
266     (z : State D) (hz : z∈cone c H) (horizon extra : ℕ) :
267     (ContinuousRecords.recordLaw (cert.instrument hH) policy z (horizon+extra)).map
268       (ContinuousRecords.truncate horizon extra)=
269       ContinuousRecords.recordLaw (cert.instrument hH) policy z horizon :=
270   ContinuousRecords.recordLaw_prefix (cert.instrument hH) policy hz horizon extra
271 end PolynomialCertificate
272 
273 end Komorebi.LSM.OrientedEllipsoid

Back to source controls