BernsteinKernel.lean

1 import Mathlib.Analysis.SpecialFunctions.Bernstein
2 import Mathlib.Analysis.Calculus.Deriv.Polynomial
3 import Mathlib.Analysis.SpecialFunctions.Integrals.Basic
4 import Mathlib.MeasureTheory.Integral.Prod
5 import Mathlib.MeasureTheory.Measure.ProbabilityMeasure
6 import Mathlib.MeasureTheory.Function.LocallyIntegrable
7 
8 /-! The actual positive Bernstein–Durrmeyer kernel on the unit interval. -/
9 set_option autoImplicit false
10 open MeasureTheory Polynomial
11 open scoped BigOperators ENNReal
12 
13 namespace Komorebi.LSM.BernsteinKernel
14 noncomputable section
15 
16 /-- Lebesgue measure on [0,1], already normalized to total mass one. -/
17 def measure : Measure ℝ := volume.restrict (Set.Icc 0 1)
18 
19 instance measure_probability : IsProbabilityMeasure measure := by
20   constructor
21   norm_num [measure]
22 
23 theorem support_ae : ∀ᵐ x ∂measure, x ∈ Set.Icc (0 : ℝ) 1 :=
24   ae_restrict_mem measurableSet_Icc
25 
26 theorem integral_eq_interval (f : ℝ → ℝ) :
27     (∫ x, f x ∂measure) = ∫ x in (0 : ℝ)..1, f x := by
28   rw [intervalIntegral.integral_of_le (by norm_num)]
29   exact integral_Icc_eq_integral_Ioc
30 
31 theorem continuous_integrable {F : Type*} [NormedAddCommGroup F] [NormedSpace ℝ F]
32     {f : ℝ → F} (hf : Continuous f) : Integrable f measure := hf.integrableOn_Icc
33 
34 /-- The basis is defined by the literal polynomial, on all real inputs. -/
35 def basis (n k : ℕ) (x : ℝ) : ℝ := (bernsteinPolynomial ℝ n k).eval x
36 
37 theorem basis_formula (n k : ℕ) (x : ℝ) :
38     basis n k x = (n.choose k : ℝ) * x^k * (1-x)^(n-k) := by
39   simp [basis, bernsteinPolynomial]
40 
41 theorem basis_continuous (n k : ℕ) : Continuous (basis n k) := (bernsteinPolynomial ℝ n k).continuous
42 
43 theorem basis_nonneg (n k : ℕ) {x : ℝ} (hx : x ∈ Set.Icc (0 : ℝ) 1) : 0 ≤ basis n k x := by
44   rw [basis_formula]
45   exact mul_nonneg (mul_nonneg (Nat.cast_nonneg _) (pow_nonneg hx.1 _))
46     (pow_nonneg (sub_nonneg.mpr hx.2) _)
47 
48 theorem sum_basis (n : ℕ) (x : ℝ) : (∑ k : Fin (n+1), basis n k x) = 1 := by
49   have h := congrArg (Polynomial.eval x) (bernsteinPolynomial.sum ℝ n)
50   simpa [basis, Polynomial.eval_finsetSum, Finset.sum_range] using h
51 
52 theorem basis_le_one (n : ℕ) (k : Fin (n+1)) {x : ℝ} (hx : x ∈ Set.Icc (0 : ℝ) 1) :
53     basis n k x ≤ 1 := by
54   rw [← sum_basis n x]
55   exact Finset.single_le_sum (s := Finset.univ) (f := fun j : Fin (n+1) => basis n j x)
56     (fun j _ => basis_nonneg n j hx) (Finset.mem_univ k)
57 
58 /-- The derivative ladder equates the integrals of adjacent basis elements. -/
59 theorem basis_integral_adjacent (n k : ℕ) (hk : k < n) :
60     (∫ x, basis n k x ∂measure) = ∫ x, basis n (k+1) x ∂measure := by
61   let p := bernsteinPolynomial ℝ (n+1) (k+1)
62   have h := intervalIntegral.integral_eq_sub_of_hasDerivAt
63     (fun x _ => p.hasDerivAt x) (p.derivative.continuous.intervalIntegrable 0 1)
64   have he : p.eval 1 - p.eval 0 = 0 := by
65     simp [p, bernsteinPolynomial.eval_at_0, bernsteinPolynomial.eval_at_1,
66       Nat.ne_of_lt hk]
67   rw [he] at h
68   have hd : (fun x => p.derivative.eval x) =
69       (fun x => ((n : ℝ)+1) * (basis n k x - basis n (k+1) x)) := by
70     funext x
71     simp [p, bernsteinPolynomial.derivative_succ_aux, basis]
72   rw [hd, intervalIntegral.integral_const_mul,
73     intervalIntegral.integral_sub ((basis_continuous n k).intervalIntegrable 0 1)
74       ((basis_continuous n (k+1)).intervalIntegrable 0 1)] at h
75   rw [integral_eq_interval, integral_eq_interval]
76   have hn : (n : ℝ)+1 ≠ 0 := by positivity
77   exact sub_eq_zero.mp ((mul_eq_zero.mp h).resolve_left hn)
78 
79 theorem basis_integral_eq_zero_index (n k : ℕ) (hk : k ≤ n) :
80     (∫ x, basis n k x ∂measure) = ∫ x, basis n 0 x ∂measure := by
81   induction k with
82   | zero => rfl
83   | succ k ih =>
84     exact (basis_integral_adjacent n k (by omega)).symm.trans (ih (by omega))
85 
86 /-- Every basis element has exact integral 1/(n+1), proved by the ladder and sum identity. -/
87 theorem basis_integral (n : ℕ) (k : Fin (n+1)) :
88     (∫ x, basis n k x ∂measure) = 1 / ((n : ℝ)+1) := by
89   have hs : (∑ j : Fin (n+1), ∫ x, basis n j x ∂measure) = 1 := by
90     rw [← integral_finsetSum (f := fun (j : Fin (n+1)) x => basis n j x) _
91       (fun j _ => continuous_integrable (basis_continuous n j))]
92     simp_rw [sum_basis]
93     simp
94   have hz : ((n : ℝ)+1) * (∫ x, basis n 0 x ∂measure) = 1 := by
95     simpa only [basis_integral_eq_zero_index n _ (Fin.is_le _), Finset.sum_const,
96       Finset.card_univ, Fintype.card_fin, nsmul_eq_mul, Nat.cast_add, Nat.cast_one] using hs
97   rw [basis_integral_eq_zero_index n k k.is_le]
98   exact (eq_div_iff (by positivity : (n : ℝ)+1 ≠ 0)).mpr (by nlinarith [hz])
99 
100 /-- Degree elevation supplies the exact first weighted beta moment. -/
101 theorem basis_mul_x (n : ℕ) (k : Fin (n+1)) (x : ℝ) :
102     basis n k x * x = ((k : ℝ)+1) / ((n : ℝ)+1) * basis (n+1) ((k : ℕ)+1) x := by
103   have hc : ((n : ℝ)+1) * (n.choose k : ℝ) =
104       ((n+1).choose ((k : ℕ)+1) : ℝ) * ((k : ℝ)+1) := by
105     exact_mod_cast Nat.add_one_mul_choose_eq n (k : ℕ)
106   rw [basis_formula, basis_formula, Nat.add_sub_add_right]
107   field_simp
108   rw [pow_succ]
109   linear_combination (x^(k : ℕ) * x * (1-x)^(n-(k : ℕ))) * hc
110 
111 theorem basis_first_moment (n : ℕ) (k : Fin (n+1)) :
112     (∫ x, basis n k x * x ∂measure) =
113       ((k : ℝ)+1) / (((n : ℝ)+1) * ((n : ℝ)+2)) := by
114   simp_rw [basis_mul_x]
115   rw [integral_const_mul, basis_integral (n+1) ⟨(k : ℕ)+1, by omega⟩]
116   push_cast
117   field_simp
118   ring
119 
120 /-- The classical Bernstein–Durrmeyer kernel. -/
121 def kernel (n : ℕ) (x v : ℝ) : ℝ :=
122   ((n : ℝ)+1) * ∑ k : Fin (n+1), basis n k x * basis n k v
123 
124 theorem kernel_nonneg (n : ℕ) {x v : ℝ}
125     (hx : x ∈ Set.Icc (0 : ℝ) 1) (hv : v ∈ Set.Icc (0 : ℝ) 1) : 0 ≤ kernel n x v := by
126   exact mul_nonneg (by positivity) (Finset.sum_nonneg fun k _ =>
127     mul_nonneg (basis_nonneg n k hx) (basis_nonneg n k hv))
128 
129 theorem kernel_symmetric (n : ℕ) (x v : ℝ) : kernel n x v = kernel n v x := by
130   simp only [kernel, mul_comm]
131 
132 theorem kernel_continuous (n : ℕ) : Continuous (Function.uncurry (kernel n)) := by
133   unfold kernel basis bernsteinPolynomial
134   fun_prop
135 
136 /-- Both marginal kernel integrals equal one exactly. -/
137 theorem kernel_integral_right (n : ℕ) (x : ℝ) : (∫ v, kernel n x v ∂measure) = 1 := by
138   simp only [kernel]
139   rw [integral_const_mul, integral_finsetSum (f := fun (k : Fin (n+1)) v => basis n k x * basis n k v) _ (fun k _ =>
140     continuous_integrable (continuous_const.mul (basis_continuous n k)))]
141   simp_rw [integral_const_mul, basis_integral]
142   rw [← Finset.sum_mul, sum_basis]
143   field_simp
144 
145 theorem kernel_integral_left (n : ℕ) (v : ℝ) : (∫ x, kernel n x v ∂measure) = 1 := by
146   simp_rw [kernel_symmetric n _ v]
147   exact kernel_integral_right n v
148 
149 /-- Continuous functions have actual integrals on the unit square. -/
150 theorem continuous_joint_integrable {F : Type*} [NormedAddCommGroup F] [NormedSpace ℝ F]
151     {f : ℝ × ℝ → F} (hf : Continuous f) : Integrable f (measure.prod measure) := by
152   unfold measure
153   rw [Measure.prod_restrict]
154   exact hf.continuousOn.integrableOn_compact (isCompact_Icc.prod isCompact_Icc)
155 
156 /-- Joint averaging separates a product test into the exact basis moments. -/
157 theorem kernel_product_integral (n : ℕ) (f g : ℝ → ℝ) (hf : Continuous f) (hg : Continuous g) :
158     (∫ x, ∫ v, kernel n x v * (f x * g v) ∂measure ∂measure) =
159       ((n : ℝ)+1) * ∑ k : Fin (n+1),
160         (∫ x, basis n k x * f x ∂measure) * (∫ v, basis n k v * g v ∂measure) := by
161   have hi (x : ℝ) : (∫ v, kernel n x v * (f x * g v) ∂measure) =
162       ((n : ℝ)+1) * ∑ k : Fin (n+1), (basis n k x * f x) *
163         (∫ v, basis n k v * g v ∂measure) := by
164     have heq : (fun v => kernel n x v * (f x * g v)) =
165         (fun v => ((n : ℝ)+1) * ∑ k : Fin (n+1), (basis n k x * f x) * (basis n k v * g v)) := by
166       funext v
167       simp only [kernel, Finset.sum_mul, Finset.mul_sum]
168       congr 1
169       funext k
170       ring
171     rw [heq, integral_const_mul, integral_finsetSum (f := fun (k : Fin (n+1)) v =>
172       (basis n k x * f x) * (basis n k v * g v)) _ (fun k _ =>
173         continuous_integrable (continuous_const.mul ((basis_continuous n k).mul hg)))]
174     simp only [integral_const_mul]
175   simp_rw [hi]
176   rw [integral_const_mul, integral_finsetSum (f := fun (k : Fin (n+1)) x =>
177     (basis n k x * f x) * (∫ v, basis n k v * g v ∂measure)) _
178     (fun k _ => continuous_integrable (((basis_continuous n k).mul hf).mul continuous_const))]
179   simp only [integral_mul_const]
180 
181 /-- The sum of the first squares, stated directly over real casts. -/
182 theorem sum_succ_squares (m : ℕ) :
183     (∑ k ∈ Finset.range m, ((k : ℝ)+1)^2) = (m : ℝ) * (m+1) * (2*m+1) / 6 := by
184   induction m with
185   | zero => simp
186   | succ m ih =>
187     rw [Finset.sum_range_succ, ih]
188     push_cast
189     ring
190 
191 /-- The exact mixed moment under the actual joint kernel density. -/
192 theorem kernel_mixed_moment (n : ℕ) :
193     (∫ x, ∫ v, kernel n x v * (x*v) ∂measure ∂measure) =
194       (2*(n : ℝ)+3) / (6*((n : ℝ)+2)) := by
195   rw [kernel_product_integral n (fun x : ℝ => x) (fun v : ℝ => v) continuous_id continuous_id]
196   simp_rw [basis_first_moment, ← pow_two, div_pow]
197   rw [← Finset.sum_div, ← Finset.sum_range (fun k : ℕ => ((k : ℝ)+1)^2), sum_succ_squares]
198   push_cast
199   field_simp
200   ring
201 
202 /-- Each marginal is uniform, so its second moment is exactly one third. -/
203 theorem kernel_first_square (n : ℕ) :
204     (∫ x, ∫ v, kernel n x v * x^2 ∂measure ∂measure) = 1/3 := by
205   simp_rw [integral_mul_const, kernel_integral_right, one_mul]
206   rw [integral_eq_interval, integral_pow]
207   norm_num
208 
209 theorem kernel_second_square (n : ℕ) :
210     (∫ x, ∫ v, kernel n x v * v^2 ∂measure ∂measure) = 1/3 := by
211   have hj : Integrable (fun p : ℝ × ℝ => kernel n p.1 p.2 * p.2^2) (measure.prod measure) :=
212     continuous_joint_integrable ((kernel_continuous n).mul (continuous_snd.pow 2))
213   rw [integral_integral_swap (f := fun x v => kernel n x v * v^2) hj]
214   simp_rw [integral_mul_const, kernel_integral_left, one_mul]
215   rw [integral_eq_interval, integral_pow]
216   norm_num
217 
218 /-- The genuine joint squared-distance integral, including the exact constant. -/
219 theorem kernel_squared_distance (n : ℕ) :
220     (∫ x, ∫ v, kernel n x v * (v-x)^2 ∂measure ∂measure) = 1 / (3*((n : ℝ)+2)) := by
221   have hi (x : ℝ) : (∫ v, kernel n x v * (v-x)^2 ∂measure) =
222       (∫ v, kernel n x v * v^2 ∂measure) +
223         (∫ v, kernel n x v * x^2 ∂measure) -
224           2 * (∫ v, kernel n x v * (x*v) ∂measure) := by
225     have heq : (fun v => kernel n x v * (v-x)^2) =
226         (fun v => kernel n x v * v^2 + kernel n x v * x^2 - 2 * (kernel n x v * (x*v))) := by
227       funext v
228       ring
229     have hc : Continuous (fun v => kernel n x v) :=
230       (kernel_continuous n).comp (continuous_const.prodMk continuous_id)
231     rw [heq, integral_sub
232       (f := fun v => kernel n x v * v^2 + kernel n x v * x^2)
233       (g := fun v => 2 * (kernel n x v * (x*v)))
234       ((continuous_integrable (hc.mul (continuous_id.pow 2))).add
235         (continuous_integrable (hc.mul continuous_const)))
236       (continuous_integrable (continuous_const.mul (hc.mul (continuous_const.mul continuous_id)))),
237       integral_add (f := fun v => kernel n x v * v^2) (g := fun v => kernel n x v * x^2)
238         (continuous_integrable (hc.mul (continuous_id.pow 2)))
239         (continuous_integrable (hc.mul continuous_const)), integral_const_mul]
240   simp_rw [hi]
241   have h1 : Integrable (fun x => ∫ v, kernel n x v * v^2 ∂measure) measure :=
242     (continuous_joint_integrable ((kernel_continuous n).mul (continuous_snd.pow 2))).integral_prod_left
243   have h2 : Integrable (fun x => ∫ v, kernel n x v * x^2 ∂measure) measure :=
244     (continuous_joint_integrable ((kernel_continuous n).mul (continuous_fst.pow 2))).integral_prod_left
245   have h3 : Integrable (fun x => ∫ v, kernel n x v * (x*v) ∂measure) measure :=
246     (continuous_joint_integrable ((kernel_continuous n).mul (continuous_fst.mul continuous_snd))).integral_prod_left
247   rw [integral_sub
248     (f := fun x => (∫ v, kernel n x v * v^2 ∂measure) + (∫ v, kernel n x v * x^2 ∂measure))
249     (g := fun x => 2 * (∫ v, kernel n x v * (x*v) ∂measure)) (h1.add h2) (h3.const_mul 2),
250     integral_add (f := fun x => ∫ v, kernel n x v * v^2 ∂measure)
251       (g := fun x => ∫ v, kernel n x v * x^2 ∂measure) h1 h2, integral_const_mul,
252     kernel_first_square, kernel_second_square, kernel_mixed_moment]
253   field_simp
254   ring
255 
256 /-- The genuine joint density is supported on the unit square. -/
257 theorem joint_support_ae : ∀ᵐ p ∂measure.prod measure,
258     p ∈ Set.Icc (0 : ℝ) 1 ×ˢ Set.Icc (0 : ℝ) 1 := by
259   apply (Measure.ae_prod_iff_ae_ae (measurableSet_Icc.prod measurableSet_Icc)).mpr
260   filter_upwards [support_ae] with x hx
261   filter_upwards [support_ae] with v hv
262   exact ⟨hx,hv⟩
263 
264 theorem kernel_nonneg_ae (n : ℕ) :
265     ∀ᵐ p ∂measure.prod measure, 0 ≤ kernel n p.1 p.2 :=
266   joint_support_ae.mono fun _ hp => kernel_nonneg n hp.1 hp.2
267 
268 theorem kernel_joint_integrable (n : ℕ) :
269     Integrable (fun p : ℝ × ℝ => kernel n p.1 p.2) (measure.prod measure) :=
270   continuous_joint_integrable (kernel_continuous n)
271 
272 theorem joint_kernel_integral (n : ℕ) :
273     (∫ p : ℝ × ℝ, kernel n p.1 p.2 ∂measure.prod measure) = 1 := by
274   rw [integral_prod (fun p : ℝ × ℝ => kernel n p.1 p.2) (kernel_joint_integrable n)]
275   simp_rw [kernel_integral_right]
276   simp
277 
278 /-- An actual probability measure, not an assumed coupling. -/
279 def jointLaw (n : ℕ) : Measure (ℝ × ℝ) :=
280   (measure.prod measure).withDensity (fun p => ENNReal.ofReal (kernel n p.1 p.2))
281 
282 instance jointLaw_probability (n : ℕ) : IsProbabilityMeasure (jointLaw n) := by
283   constructor
284   rw [jointLaw, withDensity_apply _ MeasurableSet.univ, Measure.restrict_univ,
285     ← ofReal_integral_eq_lintegral_ofReal (kernel_joint_integrable n)
286       (kernel_nonneg_ae n), joint_kernel_integral, ENNReal.ofReal_one]
287 
288 theorem jointLaw_integral (n : ℕ) (f : ℝ × ℝ → ℝ) :
289     (∫ p, f p ∂jointLaw n) = ∫ p, kernel n p.1 p.2 * f p ∂measure.prod measure := by
290   have hm : Measurable (fun p : ℝ × ℝ => ENNReal.ofReal (kernel n p.1 p.2)) :=
291     (kernel_continuous n).measurable.ennreal_ofReal
292   rw [jointLaw, integral_withDensity_eq_integral_toReal_smul hm
293     (ae_of_all _ fun _ => ENNReal.ofReal_lt_top)]
294   apply integral_congr_ae
295   filter_upwards [kernel_nonneg_ae n] with p hp
296   rw [ENNReal.toReal_ofReal hp, smul_eq_mul]
297 
298 theorem jointLaw_support (n : ℕ) : ∀ᵐ p ∂jointLaw n,
299     p ∈ Set.Icc (0 : ℝ) 1 ×ˢ Set.Icc (0 : ℝ) 1 :=
300   (withDensity_absolutelyContinuous _ _).ae_le joint_support_ae
301 
302 def distance (p : ℝ × ℝ) : ℝ := |p.2-p.1|
303 
304 theorem distance_continuous : Continuous distance := (continuous_snd.sub continuous_fst).abs
305 
306 theorem distance_le_one {p : ℝ × ℝ} (hp : p ∈ Set.Icc (0 : ℝ) 1 ×ˢ Set.Icc (0 : ℝ) 1) :
307     distance p ≤ 1 := by
308   apply abs_le.mpr
309   constructor <;> linarith [hp.1.1,hp.1.2,hp.2.1,hp.2.2]
310 
311 theorem distance_integrable (n : ℕ) : Integrable distance (jointLaw n) := by
312   apply Integrable.of_bound distance_continuous.aestronglyMeasurable 1
313   filter_upwards [jointLaw_support n] with p hp
314   change ‖|p.2-p.1|‖ ≤ 1
315   rw [Real.norm_eq_abs, abs_abs]
316   exact distance_le_one hp
317 
318 theorem distance_square_integrable (n : ℕ) : Integrable (fun p => distance p ^ 2) (jointLaw n) := by
319   apply Integrable.of_bound (distance_continuous.pow 2).aestronglyMeasurable 1
320   filter_upwards [jointLaw_support n] with p hp
321   rw [Real.norm_eq_abs, abs_of_nonneg (sq_nonneg _)]
322   exact pow_le_one₀ (abs_nonneg _) (distance_le_one hp)
323 
324 theorem jointLaw_squared_distance (n : ℕ) :
325     (∫ p, distance p ^ 2 ∂jointLaw n) = 1 / (3*((n : ℝ)+2)) := by
326   have hi : Integrable (fun p : ℝ × ℝ => kernel n p.1 p.2 * distance p ^ 2) (measure.prod measure) :=
327     continuous_joint_integrable ((kernel_continuous n).mul (distance_continuous.pow 2))
328   rw [jointLaw_integral, integral_prod (fun p : ℝ × ℝ => kernel n p.1 p.2 * distance p ^ 2) hi]
329   simp only [distance, sq_abs]
330   exact kernel_squared_distance n
331 
332 /-- Cauchy–Schwarz for the actual joint law, proved by its nonnegative variance. -/
333 theorem jointLaw_distance_bound (n : ℕ) :
334     (∫ p, distance p ∂jointLaw n) ≤ 1 / Real.sqrt (3*((n : ℝ)+2)) := by
335   let m := ∫ p, distance p ∂jointLaw n
336   have hvar : 0 ≤ ∫ p, distance p ^ 2 - 2*m*distance p + m^2 ∂jointLaw n := by
337     apply integral_nonneg
338     intro p
339     change (0 : ℝ) ≤ distance p ^ 2 - 2*m*distance p + m^2
340     nlinarith [sq_nonneg (distance p-m)]
341   rw [integral_add (f := fun p => distance p^2 - 2*m*distance p) (g := fun _ => m^2)
342       ((distance_square_integrable n).sub ((distance_integrable n).const_mul (2*m)))
343       (integrable_const _), integral_sub (f := fun p => distance p^2) (g := fun p => 2*m*distance p)
344       (distance_square_integrable n)
345       ((distance_integrable n).const_mul (2*m)), integral_const_mul, integral_const] at hvar
346   simp only [probReal_univ, one_smul] at hvar
347   rw [jointLaw_squared_distance] at hvar
348   have hm : m^2 ≤ 1 / (3*((n : ℝ)+2)) := by nlinarith [hvar]
349   have hpos : 0 ≤ 1 / (3*((n : ℝ)+2)) := by positivity
350   have hbound : m ≤ Real.sqrt (1 / (3*((n : ℝ)+2))) := by
351     nlinarith [Real.sq_sqrt hpos, Real.sqrt_nonneg (1 / (3*((n : ℝ)+2)))]
352   simpa only [one_div, Real.sqrt_inv] using hbound
353 
354 /-- The integrated absolute kernel displacement has the stated exact-rate bound. -/
355 theorem kernel_distance_bound (n : ℕ) :
356     (∫ x, ∫ v, kernel n x v * |v-x| ∂measure ∂measure) ≤
357       1 / Real.sqrt (3*((n : ℝ)+2)) := by
358   have h := jointLaw_distance_bound n
359   have hi : Integrable (fun p : ℝ × ℝ => kernel n p.1 p.2 * distance p) (measure.prod measure) :=
360     continuous_joint_integrable ((kernel_continuous n).mul distance_continuous)
361   rw [jointLaw_integral, integral_prod (fun p : ℝ × ℝ => kernel n p.1 p.2 * distance p) hi] at h
362   exact h
363 
364 end
365 end Komorebi.LSM.BernsteinKernel

Back to source controls