FiniteMomentIdentification.lean

1 import KomorebiLSM.UniformIntervalMoments
2 import KomorebiLSM.ResetCertificate
3 
4 /-! P54: conditional reconstruction from actual uniform moments, a full square
5 population core, and bounded signed reset experiments. -/
6 set_option autoImplicit false
7 open MeasureTheory Matrix
8 open scoped BigOperators
9 
10 namespace Komorebi.LSM.FiniteMomentIdentification
11 
12 section Mixing
13 variable {n : ℕ} {E : Type*} [AddCommGroup E] [Module ℝ E]
14 
15 /-- A scalar matrix acts on a finite family in any real vector space. -/
16 def mix (W : Matrix (Fin n) (Fin n) ℝ) (X : Fin n → E) (k : Fin n) : E :=
17   ∑ j, W k j • X j
18 
19 @[simp] theorem mix_one (X : Fin n → E) : mix (1 : Matrix (Fin n) (Fin n) ℝ) X = X := by
20   classical
21   funext k
22   simp [mix, Matrix.one_apply]
23 
24 theorem mix_mul (W V : Matrix (Fin n) (Fin n) ℝ) (X : Fin n → E) :
25     mix W (mix V X) = mix (W * V) X := by
26   funext k
27   simp only [mix, Finset.smul_sum, smul_smul, Matrix.mul_apply]
28   rw [Finset.sum_comm]
29   apply Finset.sum_congr rfl
30   intro l _
31   rw [Finset.sum_smul]
32 
33 theorem gram_inverse_reconstruction (D : ℕ) (X : Fin (D + 1) → E) :
34     mix (UniformIntervalMoments.gram D)⁻¹ (mix (UniformIntervalMoments.gram D) X) = X := by
35   rw [mix_mul, UniformIntervalMoments.gram_inv_mul, mix_one]
36 end Mixing
37 
38 section Population
39 variable {r D : ℕ}
40 abbrev State (r : ℕ) := Fin r → ℝ
41 abbrev Operator (r : ℕ) := Matrix (Fin r) (Fin r) ℝ
42 
43 /-- The declared degree bound is encoded by the supplied coefficient family. -/
44 def density (C : Fin (D + 1) → Operator r) (u : ℝ) : Operator r :=
45   ∑ k : Fin (D + 1), u ^ (k : ℕ) • C k
46 
47 /-- Actual population moment operators, evaluated entry by entry. -/
48 noncomputable def momentOperator (C : Fin (D + 1) → Operator r)
49     (j : Fin (D + 1)) : Operator r :=
50   fun i l => ∫ u, u ^ (j : ℕ) * density C u i l ∂UniformIntervalMoments.measure
51 
52 /-- Integrating the finite polynomial expansion derives the moment equation. -/
53 theorem momentOperator_eq (C : Fin (D + 1) → Operator r) :
54     momentOperator C = mix (UniformIntervalMoments.gram D) C := by
55   funext j
56   ext i l
57   simp only [momentOperator, density, Matrix.sum_apply, Matrix.smul_apply, smul_eq_mul,
58     Finset.mul_sum, mix, UniformIntervalMoments.gram]
59   rw [integral_finsetSum (f := fun (k : Fin (D + 1)) u => u ^ (j : ℕ) * (u ^ (k : ℕ) * C k i l)) _
60     (fun k _ => UniformIntervalMoments.continuous_integrable (by fun_prop))]
61   apply Finset.sum_congr rfl
62   intro k _
63   rw [show (fun u : ℝ => u ^ (j : ℕ) * (u ^ (k : ℕ) * C k i l)) =
64     (fun u => u ^ ((j : ℕ) + k) * C k i l) by funext u; rw [pow_add]; ring]
65   exact integral_mul_const _ _
66 
67 /-- Both factors are invertible because the population core is square in the
68 supplied true predictive dimension. -/
69 theorem core_factors_isUnit (T S : Operator r) (hH : IsUnit (T * S)) :
70     IsUnit T ∧ IsUnit S :=
71   ⟨isUnit_of_mul_isUnit_left hH, isUnit_of_mul_isUnit_right hH⟩
72 
73 noncomputable def populationCore (T S : Operator r) : Operator r := T * S
74 
75 noncomputable def populationTable (T S : Operator r)
76     (C : Fin (D + 1) → Operator r) (j : Fin (D + 1)) : Operator r :=
77   T * momentOperator C j * S
78 
79 /-- A population-table entry is the signed integral of the corresponding
80 history-to-future entry of the actual polynomial density. -/
81 theorem populationTable_entry_integral (T S : Operator r)
82     (C : Fin (D + 1) → Operator r) (j : Fin (D + 1)) (i l : Fin r) :
83     populationTable T S C j i l =
84       ∫ u, u ^ (j : ℕ) * (T * density C u * S) i l ∂UniformIntervalMoments.measure := by
85   simp only [populationTable, Matrix.mul_apply, momentOperator,
86     Finset.sum_mul, Finset.mul_sum]
87   rw [integral_finsetSum (f := fun b u =>
88     ∑ a, u ^ (j : ℕ) * (T i a * density C u a b * S b l)) _
89     (fun b _ => UniformIntervalMoments.continuous_integrable (by
90       unfold density
91       fun_prop))]
92   apply Finset.sum_congr rfl
93   intro b _
94   rw [integral_finsetSum (f := fun a u =>
95     u ^ (j : ℕ) * (T i a * density C u a b * S b l)) _
96     (fun a _ => UniformIntervalMoments.continuous_integrable (by
97       unfold density
98       fun_prop))]
99   apply Finset.sum_congr rfl
100   intro a _
101   rw [show (fun u : ℝ => u ^ (j : ℕ) * (T i a * density C u a b * S b l)) =
102       (fun u => T i a * (u ^ (j : ℕ) * density C u a b) * S b l) by
103     funext u
104     ring]
105   rw [integral_mul_const, integral_const_mul]
106 
107 noncomputable def reconstructedMoment (H : Operator r)
108     (G : Fin (D + 1) → Operator r) (j : Fin (D + 1)) : Operator r := G j * H⁻¹
109 
110 noncomputable def reconstructedCoefficient (H : Operator r)
111     (G : Fin (D + 1) → Operator r) (k : Fin (D + 1)) : Operator r :=
112   mix (UniformIntervalMoments.gram D)⁻¹ (reconstructedMoment H G) k
113 
114 /-- Exact population tables recover the moment operators in test coordinates. -/
115 theorem reconstructedMoment_eq (T S : Operator r) (hH : IsUnit (populationCore T S))
116     (C : Fin (D + 1) → Operator r) (j : Fin (D + 1)) :
117     reconstructedMoment (populationCore T S) (populationTable T S C) j =
118       T * momentOperator C j * T⁻¹ := by
119   have hS := (core_factors_isUnit T S hH).2
120   simp only [reconstructedMoment, populationTable, populationCore, Matrix.mul_inv_rev]
121   calc
122     T * momentOperator C j * S * (S⁻¹ * T⁻¹) =
123       T * momentOperator C j * (S * S⁻¹) * T⁻¹ := by simp only [Matrix.mul_assoc]
124     _ = _ := by rw [Matrix.mul_nonsing_inv S ((Matrix.isUnit_iff_isUnit_det S).mp hS)]; simp
125 
126 /-- Scalar inversion commutes with the single fixed change of predictive coordinates. -/
127 theorem mix_conjugate (W : Matrix (Fin (D + 1)) (Fin (D + 1)) ℝ)
128     (T : Operator r) (X : Fin (D + 1) → Operator r) (k : Fin (D + 1)) :
129     mix W (fun j => T * X j * T⁻¹) k = T * mix W X k * T⁻¹ := by
130   simp only [mix, Matrix.mul_sum, Matrix.sum_mul, Matrix.mul_smul, Matrix.smul_mul]
131 
132 /-- The inverse coefficient formula is derived from exact integration and the
133 proved Gram inverse; it is not an assumed moment relation. -/
134 theorem reconstructedCoefficient_eq (T S : Operator r)
135     (hH : IsUnit (populationCore T S)) (C : Fin (D + 1) → Operator r)
136     (k : Fin (D + 1)) :
137     reconstructedCoefficient (populationCore T S) (populationTable T S C) k =
138       T * C k * T⁻¹ := by
139   have hm : reconstructedMoment (populationCore T S) (populationTable T S C) =
140       (fun j => T * momentOperator C j * T⁻¹) := by
141     funext j
142     exact reconstructedMoment_eq T S hH C j
143   rw [reconstructedCoefficient, hm, mix_conjugate, momentOperator_eq]
144   rw [gram_inverse_reconstruction]
145 
146 /-- Reconstruction determines the entire declared polynomial density at every
147 observation value, rather than just its sampled moments. -/
148 theorem reconstructed_density_eq (T S : Operator r)
149     (hH : IsUnit (populationCore T S)) (C : Fin (D + 1) → Operator r) (u : ℝ) :
150     density (reconstructedCoefficient (populationCore T S) (populationTable T S C)) u =
151       T * density C u * T⁻¹ := by
152   simp only [density, reconstructedCoefficient_eq T S hH C,
153     Matrix.mul_sum, Matrix.sum_mul, Matrix.mul_smul, Matrix.smul_mul]
154 
155 /-- The supplied initial test vector is exactly the transformed initial state. -/
156 theorem initial_test_coordinates (T S : Operator r)
157     (hH : IsUnit (populationCore T S)) (z : State r) :
158     T⁻¹ *ᵥ (T *ᵥ z) = z := by
159   have hT := (core_factors_isUnit T S hH).1
160   rw [Matrix.mulVec_mulVec,
161     Matrix.nonsing_inv_mul T ((Matrix.isUnit_iff_isUnit_det T).mp hT), Matrix.one_mulVec]
162 
163 /-- Sequential polynomial operators on an unnormalized predictive state. -/
164 def wordState {A : Type*} (B : A → ℝ → Operator r) : List (A × ℝ) → State r → State r
165   | [], z => z
166   | (a, u) :: w, z => wordState B w (B a u *ᵥ z)
167 
168 /-- The recovered family and supplied initial test vector reproduce every
169 finite operator word in the same test coordinates. -/
170 theorem reconstructed_word_coordinates {A : Type*} (T S : Operator r)
171     (hH : IsUnit (populationCore T S)) (C : A → Fin (D + 1) → Operator r)
172     (w : List (A × ℝ)) (z : State r) :
173     wordState (fun a => density (reconstructedCoefficient (populationCore T S)
174       (populationTable T S (C a)))) w (T *ᵥ z) =
175       T *ᵥ wordState (fun a => density (C a)) w z := by
176   induction w generalizing z with
177   | nil => rfl
178   | cons s w ih =>
179     rcases s with ⟨a, u⟩
180     simp only [wordState, reconstructed_density_eq T S hH]
181     have hstep : (T * density (C a) u * T⁻¹) *ᵥ (T *ᵥ z) =
182         T *ᵥ (density (C a) u *ᵥ z) := by
183       rw [← Matrix.mulVec_mulVec, initial_test_coordinates T S hH,
184         ← Matrix.mulVec_mulVec]
185     rw [hstep, ih]
186 end Population
187 
188 section Error
189 variable {n : ℕ} {E : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
190 
191 /-- The inverse-row estimate holds for one fixed norm on the moment and
192 coefficient space; in particular it holds for any fixed induced operator norm. -/
193 theorem fixed_norm_reconstruction_error (W : Matrix (Fin n) (Fin n) ℝ)
194     (X Y : Fin n → E) (η : ℝ) (_hη : 0 ≤ η) (herror : ∀ j, ‖Y j - X j‖ ≤ η)
195     (k : Fin n) :
196     ‖mix W Y k - mix W X k‖ ≤ η * ∑ j, |W k j| := by
197   have heq : mix W Y k - mix W X k = ∑ j, W k j • (Y j - X j) := by
198     simp only [mix, smul_sub, Finset.sum_sub_distrib]
199   rw [heq]
200   calc
201     _ ≤ ∑ j, ‖W k j • (Y j - X j)‖ := norm_sum_le _ _
202     _ ≤ ∑ j, |W k j| * η := Finset.sum_le_sum (fun j _ => by
203       rw [norm_smul, Real.norm_eq_abs]
204       exact mul_le_mul_of_nonneg_left (herror j) (abs_nonneg _))
205     _ = _ := by rw [← Finset.sum_mul]; ring
206 
207 /-- P54's coefficient error bound for the same fixed induced norm throughout. -/
208 theorem induced_norm_coefficient_error {F : Type*} [NormedAddCommGroup F]
209     [NormedSpace ℝ F] (D : ℕ)
210     (X Y : Fin (D + 1) → (F →L[ℝ] F)) (η : ℝ) (hη : 0 ≤ η)
211     (herror : ∀ j, ‖Y j - X j‖ ≤ η) (k : Fin (D + 1)) :
212     ‖mix (UniformIntervalMoments.gram D)⁻¹ Y k - mix (UniformIntervalMoments.gram D)⁻¹ X k‖ ≤
213       η * ∑ j, |(UniformIntervalMoments.gram D)⁻¹ k j| :=
214   fixed_norm_reconstruction_error _ X Y η hη herror k
215 end Error
216 
217 section SignedReset
218 open ProbabilityTheory Real
219 variable {Ω : Type*} [MeasurableSpace Ω] {μ : Measure Ω} [IsProbabilityMeasure μ]
220 
221 /-- A reset event-word trial weighted by the designated observation power.
222 Odd powers are signed; no branch-probability normalization is applied. -/
223 noncomputable def signedResetValue (event : Set Ω) (observation : Ω → ℝ) (j : ℕ) : Ω → ℝ :=
224   event.indicator (fun ω => observation ω ^ j)
225 
226 omit [MeasurableSpace Ω] in
227 theorem signedResetValue_bounded (event : Set Ω) (observation : Ω → ℝ) (j : ℕ)
228     (hobs : ∀ ω, observation ω ∈ Set.Icc (-1 : ℝ) 1) (ω : Ω) :
229     signedResetValue event observation j ω ∈ Set.Icc (-1 : ℝ) 1 := by
230   classical
231   by_cases h : ω ∈ event
232   · rw [signedResetValue, Set.indicator_of_mem h]
233     apply abs_le.mp
234     rw [abs_pow]
235     exact pow_le_one₀ (abs_nonneg _) (abs_le.mpr (hobs ω))
236   · simp [signedResetValue, h]
237 
238 omit [IsProbabilityMeasure μ] in
239 /-- The population expectation is the event-restricted signed moment. -/
240 theorem signedResetValue_expectation (event : Set Ω) (hevent : MeasurableSet event)
241     (observation : Ω → ℝ) (j : ℕ) :
242     (∫ ω, signedResetValue event observation j ω ∂μ) =
243       ∫ ω in event, observation ω ^ j ∂μ := integral_indicator hevent
244 
245 theorem signed_centered_trial_subgaussian {X : Ω → ℝ} {p : ℝ}
246     (hm : AEMeasurable X μ) (hb : ∀ᵐ ω ∂μ, X ω ∈ Set.Icc (-1 : ℝ) 1)
247     (hp : ∫ ω, X ω ∂μ = p) :
248     HasSubgaussianMGF (fun ω => X ω - p) 1 μ := by
249   have h := hasSubgaussianMGF_of_mem_Icc hm hb
250   rw [hp] at h
251   convert h using 1
252   norm_num
253 
254 theorem signed_centered_sum_subgaussian {N : ℕ} (X : Fin N → Ω → ℝ) (p : ℝ)
255     (hind : iIndepFun X μ) (hm : ∀ i, AEMeasurable (X i) μ)
256     (hb : ∀ i, ∀ᵐ ω ∂μ, X i ω ∈ Set.Icc (-1 : ℝ) 1)
257     (hp : ∀ i, ∫ ω, X i ω ∂μ = p) :
258     HasSubgaussianMGF (fun ω => ∑ i, (X i ω - p)) (N : NNReal) μ := by
259   have hi : iIndepFun (fun i ω => X i ω - p) μ :=
260     hind.comp (fun _ x => x - p) (fun _ => by fun_prop)
261   simpa using HasSubgaussianMGF.sum_of_iIndepFun hi
262     (s := Finset.univ) (c := fun _ => (1 : NNReal))
263     (fun i _ => signed_centered_trial_subgaussian (hm i) (hb i) (hp i))
264 
265 /-- The two-sided signed-trial bound has range two and exponent `-N ε²/2`. -/
266 theorem signed_reset_error_bound {N : ℕ} (hN : 0 < N) (X : Fin N → Ω → ℝ) (p ε : ℝ)
267     (hind : iIndepFun X μ) (hm : ∀ i, AEMeasurable (X i) μ)
268     (hb : ∀ i, ∀ᵐ ω ∂μ, X i ω ∈ Set.Icc (-1 : ℝ) 1)
269     (hp : ∀ i, ∫ ω, X i ω ∂μ = p) (hε : 0 ≤ ε) :
270     μ.real {ω | ε ≤ |ResetCertificate.sampleMean X ω - p|} ≤ 2 * exp (-(N : ℝ) * ε^2 / 2) := by
271   have hn : (0 : ℝ) < N := by exact_mod_cast hN
272   have hs := signed_centered_sum_subgaussian X p hind hm hb hp
273   have hu := hs.measure_ge_le (mul_nonneg hn.le hε)
274   have hl := hs.neg.measure_ge_le (mul_nonneg hn.le hε)
275   have hexp : exp (-((N : ℝ) * ε)^2 / (2 * (N : NNReal))) =
276       exp (-(N : ℝ) * ε^2 / 2) := by
277     congr 1
278     push_cast
279     field_simp
280   rw [hexp] at hu hl
281   have hevent : {ω | ε ≤ |ResetCertificate.sampleMean X ω - p|} ⊆
282       {ω | (N : ℝ) * ε ≤ ∑ i, (X i ω - p)} ∪
283       {ω | (N : ℝ) * ε ≤ -(∑ i, (X i ω - p))} := by
284     intro ω h
285     change ε ≤ |ResetCertificate.sampleMean X ω - p| at h
286     rw [le_abs] at h
287     rcases h with h | h
288     · exact Or.inl (by
289         change (N : ℝ) * ε ≤ ∑ i, (X i ω - p)
290         rw [ResetCertificate.centered_sum_eq hN]
291         nlinarith)
292     · exact Or.inr (by
293         change (N : ℝ) * ε ≤ -(∑ i, (X i ω - p))
294         rw [ResetCertificate.centered_sum_eq hN]
295         nlinarith)
296   calc
297     _ ≤ μ.real ({ω | (N : ℝ) * ε ≤ ∑ i, (X i ω - p)} ∪
298         {ω | (N : ℝ) * ε ≤ -(∑ i, (X i ω - p))}) := measureReal_mono hevent
299     _ ≤ _ := measureReal_union_le _ _
300     _ ≤ _ := by simpa only [Pi.neg_apply, two_mul] using add_le_add hu hl
301 
302 /-- Entries are predeclared; independence is required within each entry's reset
303 trials, with no independence assumption between different entries. -/
304 theorem simultaneous_signed_error_bound {W : Type*} [Fintype W] {N : ℕ}
305     (hN : 0 < N) (X : W → Fin N → Ω → ℝ) (p : W → ℝ) (ε : ℝ)
306     (hind : ∀ w, iIndepFun (X w) μ) (hm : ∀ w i, AEMeasurable (X w i) μ)
307     (hb : ∀ w i, ∀ᵐ ω ∂μ, X w i ω ∈ Set.Icc (-1 : ℝ) 1)
308     (hp : ∀ w i, ∫ ω, X w i ω ∂μ = p w) (hε : 0 ≤ ε) :
309     μ.real {ω | ∃ w, ε < |ResetCertificate.sampleMean (X w) ω - p w|} ≤
310       2 * Fintype.card W * exp (-(N : ℝ) * ε^2 / 2) := by
311   classical
312   have hset : {ω | ∃ w, ε < |ResetCertificate.sampleMean (X w) ω - p w|} ⊆
313       ⋃ w, {ω | ε ≤ |ResetCertificate.sampleMean (X w) ω - p w|} := by
314     rintro ω ⟨w, hw⟩
315     exact Set.mem_iUnion.mpr ⟨w, hw.le⟩
316   calc
317     _ ≤ μ.real (⋃ w, {ω | ε ≤ |ResetCertificate.sampleMean (X w) ω - p w|}) :=
318       measureReal_mono hset
319     _ ≤ ∑ w, μ.real {ω | ε ≤ |ResetCertificate.sampleMean (X w) ω - p w|} :=
320       measureReal_iUnion_fintype_le _
321     _ ≤ ∑ _w : W, 2 * exp (-(N : ℝ) * ε^2 / 2) := Finset.sum_le_sum (fun w _ =>
322       signed_reset_error_bound hN (X w) (p w) ε (hind w) (hm w) (hb w) (hp w) hε)
323     _ = _ := by simp; ring
324 
325 noncomputable def signedRadius (M N : ℕ) (δ : ℝ) : ℝ :=
326   sqrt (2 * log (2 * M / δ) / N)
327 
328 omit [MeasurableSpace Ω] in
329 theorem signedRadius_budget {M N : ℕ} (hM : 0 < M) (hN : 0 < N) {δ : ℝ}
330     (hδ : 0 < δ) (hδ1 : δ < 1) :
331     2 * (M : ℝ) * exp (-(N : ℝ) * (signedRadius M N δ)^2 / 2) = δ := by
332   have hm : (1 : ℝ) ≤ M := by exact_mod_cast hM
333   have hn : (0 : ℝ) < N := by exact_mod_cast hN
334   have hratio : 0 < 2 * (M : ℝ) / δ := by positivity
335   have hrat1 : 1 ≤ 2 * (M : ℝ) / δ := (le_div_iff₀ hδ).mpr (by nlinarith)
336   have hlog : 0 ≤ log (2 * (M : ℝ) / δ) := log_nonneg hrat1
337   rw [signedRadius, sq_sqrt (div_nonneg (mul_nonneg (by norm_num) hlog) hn.le)]
338   have he : -(N : ℝ) * (2 * log (2 * M / δ) / N) / 2 = -log (2 * M / δ) := by
339     field_simp
340   rw [he, exp_neg, exp_log hratio]
341   field_simp
342 
343 /-- P54's simultaneous range-two confidence radius for actual bounded reset
344 random variables with the declared exact population expectations. -/
345 theorem signed_reset_probability_certificate {W : Type*} [Fintype W] [Nonempty W] {N : ℕ}
346     (hN : 0 < N) (X : W → Fin N → Ω → ℝ) (p : W → ℝ) {δ : ℝ}
347     (hδ : 0 < δ) (hδ1 : δ < 1)
348     (hind : ∀ w, iIndepFun (X w) μ) (hm : ∀ w i, Measurable (X w i))
349     (hb : ∀ w i, ∀ᵐ ω ∂μ, X w i ω ∈ Set.Icc (-1 : ℝ) 1)
350     (hp : ∀ w i, ∫ ω, X w i ω ∂μ = p w) :
351     1 - δ ≤ μ.real {ω | ∀ w, |ResetCertificate.sampleMean (X w) ω - p w| ≤
352       signedRadius (Fintype.card W) N δ} := by
353   have hbad := simultaneous_signed_error_bound hN X p
354     (signedRadius (Fintype.card W) N δ) hind
355     (fun w i => (hm w i).aemeasurable) hb hp (sqrt_nonneg _)
356   rw [signedRadius_budget Fintype.card_pos hN hδ hδ1] at hbad
357   have hmeas : MeasurableSet {ω | ∃ w, signedRadius (Fintype.card W) N δ <
358       |ResetCertificate.sampleMean (X w) ω - p w|} := by
359     simp only [Set.ofPred_exists]
360     apply MeasurableSet.iUnion
361     intro w
362     apply measurableSet_lt measurable_const
363     unfold ResetCertificate.sampleMean
364     fun_prop
365   have hcomp : {ω | ∀ w, |ResetCertificate.sampleMean (X w) ω - p w| ≤
366       signedRadius (Fintype.card W) N δ} =
367       {ω | ∃ w, signedRadius (Fintype.card W) N δ <
368         |ResetCertificate.sampleMean (X w) ω - p w|}ᶜ := by
369     ext ω
370     simp
371   rw [hcomp, measureReal_compl hmeas]
372   simpa using (by linarith : 1 - δ ≤ 1 - μ.real {ω | ∃ w,
373     signedRadius (Fintype.card W) N δ < |ResetCertificate.sampleMean (X w) ω - p w|})
374 end SignedReset
375 
376 end Komorebi.LSM.FiniteMomentIdentification

Back to source controls