UniformIntervalMoments.lean

1 import KomorebiLSM.AveragedCongruence
2 import Mathlib.Algebra.Polynomial.OfFn
3 import Mathlib.Algebra.Polynomial.Roots
4 
5 /-! Exact power moments and positive polynomial Gram matrices for the actual
6 normalized Lebesgue measure on `[-1,1]`. Endpoints have zero Lebesgue mass. -/
7 set_option autoImplicit false
8 open MeasureTheory Matrix Polynomial
9 open scoped BigOperators ENNReal
10 
11 namespace Komorebi.LSM.UniformIntervalMoments
12 
13 noncomputable def measure : Measure ℝ := AveragedCongruence.uniformInterval (-1) 1
14 
15 instance measure_isProbability : IsProbabilityMeasure measure :=
16   AveragedCongruence.uniformInterval_isProbability (by norm_num)
17 
18 /-- Integration is against the specified probability measure, not a formal moment oracle. -/
19 theorem integral_eq (f : ℝ → ℝ) :
20     (∫ x, f x ∂measure) = (1 / 2 : ℝ) * ∫ x in (-1 : ℝ)..1, f x := by
21   rw [measure, AveragedCongruence.integral_uniformInterval (by norm_num)]
22   norm_num
23 
24 theorem continuous_integrable {f : ℝ → ℝ} (hf : Continuous f) : Integrable f measure := by
25   unfold measure AveragedCongruence.uniformInterval ProbabilityTheory.cond
26   exact (hf.intervalIntegrable (-1) 1).1.smul_measure (by simp)
27 
28 noncomputable def moment (k : ℕ) : ℝ := ∫ x, x ^ k ∂measure
29 
30 theorem moment_formula (k : ℕ) :
31     moment k = (1 - (-1 : ℝ) ^ (k + 1)) / (2 * (k + 1)) := by
32   rw [moment, integral_eq, integral_pow]
33   simp only [one_pow]
34   field_simp
35 
36 @[simp] theorem moment_even {k : ℕ} (hk : Even k) : moment k = 1 / (k + 1) := by
37   rw [moment_formula, pow_succ, hk.neg_pow]
38   simp only [one_pow, one_mul]
39   field_simp
40   norm_num
41 
42 @[simp] theorem moment_odd {k : ℕ} (hk : Odd k) : moment k = 0 := by
43   rw [moment_formula, pow_succ, hk.neg_pow]
44   simp
45 
46 @[simp] theorem moment_zero : moment 0 = 1 := by norm_num [moment_formula]
47 
48 @[simp] theorem moment_one : moment 1 = 0 := by norm_num [moment_formula]
49 
50 @[simp] theorem moment_two : moment 2 = 1 / 3 := by norm_num [moment_formula]
51 
52 /-- A finite power expansion has the expected exact integral. -/
53 theorem integral_power_sum {n : ℕ} (c : Fin n → ℝ) :
54     (∫ x, ∑ k : Fin n, c k * x ^ (k : ℕ) ∂measure) =
55       ∑ k : Fin n, c k * moment k := by
56   rw [integral_finsetSum (f := fun k x => c k * x ^ (k : ℕ)) _ (fun k _ =>
57     continuous_integrable (continuous_const.mul (continuous_id.pow k))) ]
58   simp only [integral_const_mul, moment]
59 
60 /-- Finite coefficient vectors are represented by genuine polynomials. -/
61 theorem eval_ofFn {n : ℕ} (c : Fin n → ℝ) (x : ℝ) :
62     (Polynomial.ofFn n c).eval x = ∑ k : Fin n, c k * x ^ (k : ℕ) := by
63   rw [Polynomial.ofFn_eq_sum_monomial]
64   simp only [Polynomial.eval_finsetSum, Polynomial.eval_monomial]
65 
66 /-- Every nonzero polynomial has a positive squared integral for this measure. -/
67 theorem polynomial_square_integral_pos (p : ℝ[X]) (hp : p ≠ 0) :
68     0 < ∫ x, (p.eval x) ^ 2 ∂measure := by
69   obtain ⟨x, hx, hroot⟩ := ((Set.Ioo_infinite (by norm_num : (-1 : ℝ) < 1)).sdiff
70     (Polynomial.finite_setOfPred_isRoot hp)).nonempty
71   rw [integral_eq]
72   apply mul_pos (by norm_num)
73   apply intervalIntegral.integral_pos (by norm_num)
74   · exact (p.continuous.pow 2).continuousOn
75   · intro y _
76     positivity
77   · exact ⟨x, ⟨hx.1.le, hx.2.le⟩, sq_pos_of_ne_zero hroot⟩
78 
79 /-- The `(D+1)` by `(D+1)` moment Gram matrix. -/
80 noncomputable def gram (D : ℕ) : Matrix (Fin (D + 1)) (Fin (D + 1)) ℝ :=
81   fun j k => moment ((j : ℕ) + k)
82 
83 /-- Its quadratic form is exactly the integral of the squared coefficient polynomial. -/
84 theorem gram_quadratic (D : ℕ) (c : Fin (D + 1) → ℝ) :
85     star c ⬝ᵥ (gram D *ᵥ c) =
86       ∫ x, ((Polynomial.ofFn (D + 1) c).eval x) ^ 2 ∂measure := by
87   simp only [dotProduct, Matrix.mulVec, gram, Pi.star_apply, star_trivial]
88   simp_rw [Finset.mul_sum]
89   rw [show (fun x => ((Polynomial.ofFn (D + 1) c).eval x) ^ 2) =
90       (fun x => ∑ j : Fin (D + 1), ∑ k : Fin (D + 1),
91         (c j * c k) * x ^ ((j : ℕ) + k)) by
92     funext x
93     rw [eval_ofFn, pow_two, Finset.sum_mul]
94     apply Finset.sum_congr rfl
95     intro j _
96     rw [Finset.mul_sum]
97     apply Finset.sum_congr rfl
98     intro k _
99     rw [pow_add]
100     ring]
101   rw [integral_finsetSum _ (fun j _ => continuous_integrable (by fun_prop))]
102   apply Finset.sum_congr rfl
103   intro j _
104   rw [integral_finsetSum (f := fun k x => (c j * c k) * x ^ ((j : ℕ) + k)) _
105     (fun k _ => continuous_integrable (by fun_prop))]
106   apply Finset.sum_congr rfl
107   intro k _
108   rw [integral_const_mul]
109   change c j * (moment ((j : ℕ) + k) * c k) = (c j * c k) * moment ((j : ℕ) + k)
110   ring
111 
112 /-- Positive definiteness is proved from the nonzero-polynomial square integral. -/
113 theorem gram_posDef (D : ℕ) : (gram D).PosDef := by
114   apply Matrix.PosDef.of_dotProduct_mulVec_pos
115   · ext j k
116     simp [gram, Nat.add_comm]
117   · intro c hc
118     rw [gram_quadratic]
119     apply polynomial_square_integral_pos
120     intro hz
121     apply hc
122     apply (Polynomial.injective_ofFn (D + 1))
123     simpa using hz
124 
125 theorem gram_isUnit (D : ℕ) : IsUnit (gram D) := (gram_posDef D).isUnit
126 
127 theorem gram_mul_inv (D : ℕ) : gram D * (gram D)⁻¹ = 1 :=
128   Matrix.mul_nonsing_inv _ ((Matrix.isUnit_iff_isUnit_det _).mp (gram_isUnit D))
129 
130 theorem gram_inv_mul (D : ℕ) : (gram D)⁻¹ * gram D = 1 :=
131   Matrix.nonsing_inv_mul _ ((Matrix.isUnit_iff_isUnit_det _).mp (gram_isUnit D))
132 
133 end Komorebi.LSM.UniformIntervalMoments

Back to source controls