PositivePolynomialApproximation.lean

1 import KomorebiLSM.BernsteinKernel
2 import KomorebiLSM.BaseNormIntegration
3 import KomorebiLSM.ContinuousInstrument
4 import KomorebiLSM.PositiveCellApproximation
5 import KomorebiLSM.ContinuousTV
6 
7 /-! Positive polynomial approximation from actual integral coefficients. No
8 polynomial degree hypothesis is imposed on the original operator density. -/
9 set_option autoImplicit false
10 open MeasureTheory
11 open scoped BigOperators
12 
13 namespace Komorebi.LSM.PositivePolynomialApproximation
14 noncomputable section
15 
16 section Averaging
17 variable {F : Type*} [NormedAddCommGroup F] [NormedSpace ℝ F] [CompleteSpace F]
18 
19 /-- The exact positive integral coefficients, with the classical normalization. -/
20 def coefficient (n : ℕ) (B : ℝ → F) (k : Fin (n+1)) : F :=
21   ((n : ℝ)+1) • ∫ v, BernsteinKernel.basis n k v • B v ∂BernsteinKernel.measure
22 
23 /-- The approximating polynomial is given by its finite Bernstein expansion. -/
24 def approximation (n : ℕ) (B : ℝ → F) (x : ℝ) : F :=
25   ∑ k : Fin (n+1), BernsteinKernel.basis n k x • coefficient n B k
26 
27 omit [CompleteSpace F] in
28 theorem basis_weighted_integrable (n : ℕ) (B : ℝ → F) (hB : Integrable B BernsteinKernel.measure)
29     (k : Fin (n+1)) : Integrable (fun v => BernsteinKernel.basis n k v • B v) BernsteinKernel.measure := by
30   apply hB.bdd_smul 1 (BernsteinKernel.basis_continuous n k).aestronglyMeasurable
31   filter_upwards [BernsteinKernel.support_ae] with v hv
32   rw [Real.norm_eq_abs, abs_of_nonneg (BernsteinKernel.basis_nonneg n k hv)]
33   exact BernsteinKernel.basis_le_one n k hv
34 
35 omit [CompleteSpace F] in
36 theorem approximation_continuous (n : ℕ) (B : ℝ → F) : Continuous (approximation n B) := by
37   unfold approximation
38   exact continuous_finsetSum _ fun k _ => (BernsteinKernel.basis_continuous n k).smul continuous_const
39 
40 omit [CompleteSpace F] in
41 theorem approximation_integrable (n : ℕ) (B : ℝ → F) :
42     Integrable (approximation n B) BernsteinKernel.measure :=
43   BernsteinKernel.continuous_integrable (approximation_continuous n B)
44 
45 omit [CompleteSpace F] in
46 /-- Kernel multiplication preserves integrability by its finite basis expansion. -/
47 theorem kernel_weighted_integrable (n : ℕ) (B : ℝ → F) (hB : Integrable B BernsteinKernel.measure) (x : ℝ) :
48     Integrable (fun v => BernsteinKernel.kernel n x v • B v) BernsteinKernel.measure := by
49   have heq : (fun v => BernsteinKernel.kernel n x v • B v) =
50       (fun v => ((n : ℝ)+1) • ∑ k : Fin (n+1),
51         BernsteinKernel.basis n k x • (BernsteinKernel.basis n k v • B v)) := by
52     funext v
53     simp only [BernsteinKernel.kernel, mul_smul, Finset.sum_smul]
54   rw [heq]
55   exact Integrable.smul _ (integrable_finsetSum _ fun k _ =>
56     Integrable.smul _ (basis_weighted_integrable n B hB k))
57 
58 omit [CompleteSpace F] in
59 /-- Kernel averaging and the finite Bernstein coefficient formula coincide. -/
60 theorem approximation_eq_kernel (n : ℕ) (B : ℝ → F) (hB : Integrable B BernsteinKernel.measure) (x : ℝ) :
61     approximation n B x = ∫ v, BernsteinKernel.kernel n x v • B v ∂BernsteinKernel.measure := by
62   have heq : (fun v => BernsteinKernel.kernel n x v • B v) =
63       (fun v => ((n : ℝ)+1) • ∑ k : Fin (n+1),
64         BernsteinKernel.basis n k x • (BernsteinKernel.basis n k v • B v)) := by
65     funext v
66     simp only [BernsteinKernel.kernel, mul_smul, Finset.sum_smul]
67   calc
68     approximation n B x = ((n : ℝ)+1) • ∑ k : Fin (n+1),
69         BernsteinKernel.basis n k x • (∫ v, BernsteinKernel.basis n k v • B v ∂BernsteinKernel.measure) := by
70       simp only [approximation, coefficient, Finset.smul_sum, smul_smul]
71       apply Finset.sum_congr rfl
72       intro k _
73       rw [mul_comm]
74     _ = _ := by
75       rw [heq, integral_smul, integral_finsetSum
76         (f := fun (k : Fin (n+1)) v => BernsteinKernel.basis n k x • (BernsteinKernel.basis n k v • B v)) _
77         (fun k _ => Integrable.smul _ (basis_weighted_integrable n B hB k))]
78       congr 1
79       apply Finset.sum_congr rfl
80       intro k _
81       exact (integral_smul _ _).symm
82 
83 /-- Double stochasticity preserves the original integrated operator exactly. -/
84 theorem approximation_integral (n : ℕ) (B : ℝ → F) (hB : Integrable B BernsteinKernel.measure) :
85     (∫ x, approximation n B x ∂BernsteinKernel.measure) = ∫ v, B v ∂BernsteinKernel.measure := by
86   simp only [approximation]
87   rw [integral_finsetSum
88     (f := fun (k : Fin (n+1)) x => BernsteinKernel.basis n k x • coefficient n B k) _
89     (fun k _ => (BernsteinKernel.continuous_integrable (BernsteinKernel.basis_continuous n k)).smul_const _)]
90   have hi (k : Fin (n+1)) :
91       (∫ x, BernsteinKernel.basis n k x • coefficient n B k ∂BernsteinKernel.measure) =
92         ∫ v, BernsteinKernel.basis n k v • B v ∂BernsteinKernel.measure := by
93     rw [integral_smul_const, BernsteinKernel.basis_integral n k, coefficient, smul_smul]
94     have hm : (1 / ((n : ℝ)+1)) * ((n : ℝ)+1) = 1 := by field_simp
95     rw [hm, one_smul]
96   simp_rw [hi]
97   rw [← integral_finsetSum (f := fun (k : Fin (n+1)) v => BernsteinKernel.basis n k v • B v) _
98     (fun k _ => basis_weighted_integrable n B hB k)]
99   apply integral_congr_ae
100   filter_upwards [] with v
101   rw [← Finset.sum_smul, BernsteinKernel.sum_basis, one_smul]
102 /-- Constant operator densities retain their exact coefficient normalization. -/
103 theorem coefficient_constant (n : ℕ) (c : F) (k : Fin (n+1)) :
104     coefficient n (fun _ => c) k = c := by
105   rw [coefficient, integral_smul_const, BernsteinKernel.basis_integral n k, smul_smul]
106   have hn : ((n : ℝ)+1) * (1 / ((n : ℝ)+1)) = 1 := by field_simp
107   rw [hn, one_smul]
108 
109 theorem approximation_constant (n : ℕ) (c : F) (x : ℝ) :
110     approximation n (fun _ => c) x = c := by
111   simp only [approximation, coefficient_constant]
112   rw [← Finset.sum_smul, BernsteinKernel.sum_basis, one_smul]
113 end Averaging
114 
115 /-- Every Bernstein basis polynomial has the declared degree bound. -/
116 theorem basis_degree_le (n : ℕ) (k : Fin (n+1)) :
117     (bernsteinPolynomial ℝ n k).natDegree ≤ n := by
118   have hlin : (1 - Polynomial.X : Polynomial ℝ).natDegree ≤ 1 := by
119     simpa using Polynomial.natDegree_sub_le (1 : Polynomial ℝ) Polynomial.X
120   unfold bernsteinPolynomial
121   calc
122     _ ≤ (n.choose k : Polynomial ℝ).natDegree + (Polynomial.X^(k : ℕ) : Polynomial ℝ).natDegree +
123         ((1 - Polynomial.X)^(n-(k : ℕ)) : Polynomial ℝ).natDegree :=
124       Polynomial.natDegree_mul_le.trans (Nat.add_le_add_right Polynomial.natDegree_mul_le _)
125     _ ≤ 0 + (k : ℕ) + (n-(k : ℕ)) := by
126       apply Nat.add_le_add
127       · simp
128       · exact (Polynomial.natDegree_pow_le_of_le _ hlin).trans_eq (Nat.mul_one _)
129     _ = n := by omega
130 
131 theorem basis_power_expansion (n : ℕ) (k : Fin (n+1)) (x : ℝ) :
132     BernsteinKernel.basis n k x = ∑ m : Fin (n+1),
133       (bernsteinPolynomial ℝ n k).coeff m * x^(m : ℕ) := by
134   have h := congrArg (Polynomial.eval x)
135     (Polynomial.as_sum_range_C_mul_X_pow' (bernsteinPolynomial ℝ n k)
136       (Nat.lt_succ_of_le (basis_degree_le n k)))
137   simpa [BernsteinKernel.basis, Polynomial.eval_finsetSum, Finset.sum_range] using h
138 
139 section Degree
140 variable {F : Type*} [NormedAddCommGroup F] [NormedSpace ℝ F] [CompleteSpace F]
141 
142 def powerCoefficient (n : ℕ) (B : ℝ → F) (m : Fin (n+1)) : F :=
143   ∑ k : Fin (n+1), (bernsteinPolynomial ℝ n k).coeff m • coefficient n B k
144 
145 omit [CompleteSpace F] in
146 /-- An explicit power expansion certifies degree at most n in the vector space. -/
147 theorem approximation_power_expansion (n : ℕ) (B : ℝ → F) (x : ℝ) :
148     approximation n B x = ∑ m : Fin (n+1), x^(m : ℕ) • powerCoefficient n B m := by
149   simp only [approximation, powerCoefficient, basis_power_expansion, Finset.sum_smul,
150     Finset.smul_sum, smul_smul]
151   rw [Finset.sum_comm]
152   apply Finset.sum_congr rfl
153   intro m _
154   apply Finset.sum_congr rfl
155   intro k _
156   rw [mul_comm]
157 end Degree
158 
159 section Operators
160 variable {E A : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E] [CompleteSpace E]
161 variable [FiniteDimensional ℝ E] [MeasurableSpace E] [BorelSpace E]
162 variable (C : BaseNorm.ProbabilityCone E)
163 
164 omit [FiniteDimensional ℝ E] [MeasurableSpace E] [BorelSpace E] in
165 /-- Every integral coefficient is a positive operator on the original cone. -/
166 theorem coefficient_positive (n : ℕ) (B : ℝ → E →L[ℝ] E)
167     (hB : Integrable B BernsteinKernel.measure)
168     (hpos : ∀ᵐ v ∂BernsteinKernel.measure, ∀ z ∈ C.carrier, B v z ∈ C.carrier)
169     (k : Fin (n+1)) (z : E) (hz : z ∈ C.carrier) : coefficient n B k z ∈ C.carrier := by
170   have hi := basis_weighted_integrable n B hB k
171   simp only [coefficient, smul_apply, ContinuousLinearMap.integral_apply hi]
172   apply C.smul_mem (by positivity)
173   apply C.convex_carrier.integral_mem C.closed _ ((ContinuousLinearMap.apply ℝ E z).integrable_comp hi)
174   filter_upwards [BernsteinKernel.support_ae, hpos] with v hv hBv
175   exact C.smul_mem (BernsteinKernel.basis_nonneg n k hv) (hBv z hz)
176 
177 omit [FiniteDimensional ℝ E] [MeasurableSpace E] [BorelSpace E] in
178 /-- The entire approximating polynomial is cone positive at every unit input. -/
179 theorem approximation_positive (n : ℕ) (B : ℝ → E →L[ℝ] E)
180     (hB : Integrable B BernsteinKernel.measure)
181     (hpos : ∀ᵐ v ∂BernsteinKernel.measure, ∀ z ∈ C.carrier, B v z ∈ C.carrier)
182     {x : ℝ} (hx : x ∈ Set.Icc (0 : ℝ) 1) (z : E) (hz : z ∈ C.carrier) :
183     approximation n B x z ∈ C.carrier := by
184   simp only [approximation, sum_apply, smul_apply]
185   apply Finset.sum_induction
186   · exact fun _ _ => C.add_mem
187   · exact C.zero_mem
188   · intro k _
189     exact C.smul_mem (BernsteinKernel.basis_nonneg n k hx) (coefficient_positive C n B hB hpos k z hz)
190 
191 omit [FiniteDimensional ℝ E] [MeasurableSpace E] [BorelSpace E] in
192 /-- The error subtracts the original operator inside the normalized kernel. -/
193 theorem approximation_error_identity (n : ℕ) (B : ℝ → E →L[ℝ] E)
194     (hB : Integrable B BernsteinKernel.measure) (x : ℝ) (z : E) :
195     approximation n B x z - B x z =
196       ∫ v, BernsteinKernel.kernel n x v • (B v z - B x z) ∂BernsteinKernel.measure := by
197   have hi := kernel_weighted_integrable n B hB x
198   have ha : Integrable (fun v => BernsteinKernel.kernel n x v • B v z) BernsteinKernel.measure :=
199     (ContinuousLinearMap.apply ℝ E z).integrable_comp hi
200   have hc : Integrable (fun v => BernsteinKernel.kernel n x v • B x z) BernsteinKernel.measure :=
201     BernsteinKernel.continuous_integrable
202       (((BernsteinKernel.kernel_continuous n).comp (continuous_const.prodMk continuous_id)).smul continuous_const)
203   simp_rw [smul_sub]
204   rw [integral_sub ha hc, integral_smul_const, BernsteinKernel.kernel_integral_right, one_smul,
205     approximation_eq_kernel n B hB, ContinuousLinearMap.integral_apply hi z]
206   simp only [smul_apply]
207 
208 omit [MeasurableSpace E] [BorelSpace E] in
209 /-- A fixed normalized input obeys the kernel displacement bound pointwise. -/
210 theorem approximation_error_point (n : ℕ) (B : ℝ → E →L[ℝ] E)
211     (hB : Integrable B BernsteinKernel.measure) (L : ℝ) (_hL : 0 ≤ L)
212     {x : ℝ} (hx : x ∈ Set.Icc (0 : ℝ) 1) (z : E)
213     (hLip : ∀ v ∈ Set.Icc (0 : ℝ) 1, C.value (B v z-B x z) ≤ L*|v-x|) :
214     C.value (approximation n B x z-B x z) ≤
215       L * ∫ v, BernsteinKernel.kernel n x v * |v-x| ∂BernsteinKernel.measure := by
216   have hf : Integrable (fun v => B v z-B x z) BernsteinKernel.measure :=
217     ((ContinuousLinearMap.apply ℝ E z).integrable_comp hB).sub (integrable_const _)
218   have hi := kernel_weighted_integrable n (fun v => B v z-B x z) hf x
219   rw [approximation_error_identity n B hB x z]
220   calc
221     _ ≤ ∫ v, C.value (BernsteinKernel.kernel n x v • (B v z-B x z)) ∂BernsteinKernel.measure :=
222       C.value_integral_le_integral_value hi
223     _ ≤ ∫ v, L * (BernsteinKernel.kernel n x v * |v-x|) ∂BernsteinKernel.measure := by
224       apply integral_mono_ae (C.integrable_value hi)
225         (BernsteinKernel.continuous_integrable (by
226           have hc : Continuous (fun v : ℝ => BernsteinKernel.kernel n x v) :=
227             (BernsteinKernel.kernel_continuous n).comp (continuous_const.prodMk continuous_id)
228           exact continuous_const.mul (hc.mul (continuous_id.sub continuous_const).abs)))
229       filter_upwards [BernsteinKernel.support_ae] with v hv
230       rw [C.value_smul, abs_of_nonneg (BernsteinKernel.kernel_nonneg n hx hv)]
231       have h := mul_le_mul_of_nonneg_left (hLip v hv) (BernsteinKernel.kernel_nonneg n hx hv)
232       nlinarith [h]
233     _ = _ := integral_const_mul _ _
234 
235 omit [MeasurableSpace E] [BorelSpace E] in
236 /-- The integrated error has the universal Bernstein–Durrmeyer rate. -/
237 theorem approximation_integrated_error (n : ℕ) (B : ℝ → E →L[ℝ] E)
238     (hB : Integrable B BernsteinKernel.measure) (L : ℝ) (hL : 0 ≤ L) (z : E)
239     (hLip : ∀ x ∈ Set.Icc (0 : ℝ) 1, ∀ v ∈ Set.Icc (0 : ℝ) 1,
240       C.value (B v z-B x z) ≤ L*|v-x|) :
241     (∫ x, C.value (approximation n B x z-B x z) ∂BernsteinKernel.measure) ≤
242       L / Real.sqrt (3*((n : ℝ)+2)) := by
243   have hleft := C.integrable_value
244     (((ContinuousLinearMap.apply ℝ E z).integrable_comp (approximation_integrable n B)).sub
245       ((ContinuousLinearMap.apply ℝ E z).integrable_comp hB))
246   have hj : Integrable (fun p : ℝ × ℝ => BernsteinKernel.kernel n p.1 p.2 * |p.2-p.1|)
247       (BernsteinKernel.measure.prod BernsteinKernel.measure) :=
248     BernsteinKernel.continuous_joint_integrable
249       ((BernsteinKernel.kernel_continuous n).mul (continuous_snd.sub continuous_fst).abs)
250   have hright : Integrable (fun x => ∫ v, BernsteinKernel.kernel n x v * |v-x| ∂BernsteinKernel.measure)
251       BernsteinKernel.measure := hj.integral_prod_left
252   calc
253     _ ≤ ∫ x, L * (∫ v, BernsteinKernel.kernel n x v * |v-x| ∂BernsteinKernel.measure)
254         ∂BernsteinKernel.measure := by
255       apply integral_mono_ae hleft (hright.const_mul L)
256       filter_upwards [BernsteinKernel.support_ae] with x hx
257       exact approximation_error_point C n B hB L hL hx z (hLip x hx)
258     _ = L * (∫ x, ∫ v, BernsteinKernel.kernel n x v * |v-x|
259         ∂BernsteinKernel.measure ∂BernsteinKernel.measure) := integral_const_mul _ _
260     _ ≤ L * (1 / Real.sqrt (3*((n : ℝ)+2))) :=
261       mul_le_mul_of_nonneg_left (BernsteinKernel.kernel_distance_bound n) hL
262     _ = _ := by ring
263 
264 variable {C}
265 
266 /-- A genuine probability instrument, with exact mass rather than a correction. -/
267 def approximatingInstrument (n : ℕ)
268     (D : ContinuousInstrument.DensityInstrument (A := A) C BernsteinKernel.measure) :
269     ContinuousInstrument.DensityInstrument (A := A) C BernsteinKernel.measure where
270   op a := approximation n (D.op a)
271   measurable_op a := (approximation_continuous n (D.op a)).measurable
272   integrable_op a := approximation_integrable n (D.op a)
273   positive a := BernsteinKernel.support_ae.mono fun _ hx z hz =>
274     approximation_positive C n (D.op a) (D.integrable_op a) (D.positive a) hx z hz
275   mass_preserving a z := by
276     rw [approximation_integral n (D.op a) (D.integrable_op a)]
277     exact D.mass_preserving a z
278 
279 omit [MeasurableSpace E] [BorelSpace E] in
280 /-- The source's Lipschitz premise uses the actual induced base norm. -/
281 theorem induced_lipschitz_approximation_error (n : ℕ) (B : ℝ → E →L[ℝ] E)
282     (hB : Integrable B BernsteinKernel.measure) (L : ℝ) (hL : 0 ≤ L)
283     (hLip : ∀ x ∈ Set.Icc (0 : ℝ) 1, ∀ v ∈ Set.Icc (0 : ℝ) 1,
284       PositiveCellApproximation.inducedNorm (C := C) (B v-B x) ≤ L*|v-x|)
285     (z : E) (hz : z ∈ C.base) :
286     (∫ x, C.value (approximation n B x z-B x z) ∂BernsteinKernel.measure) ≤
287       L / Real.sqrt (3*((n : ℝ)+2)) := by
288   apply approximation_integrated_error C n B hB L hL z
289   intro x hx v hv
290   have h := PositiveCellApproximation.value_apply_le_inducedNorm (C := C) (B v-B x) z
291   change C.value (B v z-B x z) ≤ _ at h
292   rw [C.value_on_cone hz.1, hz.2, mul_one] at h
293   exact h.trans (hLip x hx v hv)
294 
295 omit [MeasurableSpace E] [BorelSpace E] in
296 /-- Uniform integrated discrepancy of the actual valid instruments, in the
297 orientation used by the common-policy law comparison. -/
298 theorem approximating_instrument_discrepancy (n : ℕ)
299     (D : ContinuousInstrument.DensityInstrument (A := A) C BernsteinKernel.measure)
300     (L : ℝ) (hL : 0 ≤ L)
301     (hLip : ∀ a x, x ∈ Set.Icc (0 : ℝ) 1 → ∀ v ∈ Set.Icc (0 : ℝ) 1,
302       PositiveCellApproximation.inducedNorm (C := C) (D.op a v-D.op a x) ≤ L*|v-x|)
303     (a : A) (z : E) (hz : z ∈ C.base) :
304     (∫ x, C.value (D.op a x z-(approximatingInstrument n D).op a x z) ∂BernsteinKernel.measure) ≤
305       L / Real.sqrt (3*((n : ℝ)+2)) := by
306   have heq : (fun x => C.value (D.op a x z-(approximatingInstrument n D).op a x z)) =
307       (fun x => C.value (approximation n (D.op a) x z-D.op a x z)) := by
308     funext x
309     change C.value (D.op a x z-approximation n (D.op a) x z) = _
310     rw [← C.value_neg (D.op a x z-approximation n (D.op a) x z), neg_sub]
311   rw [heq]
312   exact induced_lipschitz_approximation_error n (D.op a) (D.integrable_op a) L hL (hLip a) z hz
313 /-- P64's length-T bound compares the actual laws under the same reported-history
314 policy and identical initial states. It invokes the full P59 comparison theorem. -/
315 theorem common_policy_law_error [MeasurableSpace A] [Fintype A] [MeasurableSingletonClass A]
316     (n : ℕ) (D : ContinuousInstrument.DensityInstrument (A := A) C BernsteinKernel.measure)
317     (L : ℝ) (hL : 0 ≤ L)
318     (hLip : ∀ a x, x ∈ Set.Icc (0 : ℝ) 1 → ∀ v ∈ Set.Icc (0 : ℝ) 1,
319       PositiveCellApproximation.inducedNorm (C := C) (D.op a v-D.op a x) ≤ L*|v-x|)
320     (π : ContinuousRecords.Policy (A := A) (Y := ℝ))
321     (z : ContinuousInstrument.DensityInstrument.State C) (T : ℕ) :
322     SequentialLowerCertificate.lawTV
323       (TimeRecords.recordLaw (fun _ => D) π z.1 T)
324       (TimeRecords.recordLaw (fun _ => approximatingInstrument n D) π z.1 T) ≤
325         min 1 ((T : ℝ)*L/(2*Real.sqrt (3*((n : ℝ)+2)))) := by
326   have h := ContinuousTV.record_lawTV_le (fun _ => D) (fun _ => approximatingInstrument n D)
327     π z z (fun _ => L / Real.sqrt (3*((n : ℝ)+2)))
328     (fun _ => div_nonneg hL (Real.sqrt_nonneg _))
329     (fun _ a u hu => approximating_instrument_discrepancy n D L hL hLip a u hu)
330     0 (by simp) T
331   convert h using 1
332   congr 1
333   simp only [Finset.sum_const, Finset.card_range, nsmul_eq_mul, zero_add]
334   ring
335 end Operators
336 
337 end
338 end Komorebi.LSM.PositivePolynomialApproximation

Back to source controls