RegressionGeometry.lean

1 import Mathlib.Analysis.Matrix.PosDef
2 import Mathlib.Tactic
3 
4 /-! Deterministic geometry for actual ridge-regression confidence ellipsoids. -/
5 set_option autoImplicit false
6 namespace Komorebi.LSM.RegressionConfidence
7 open Matrix
8 open scoped BigOperators ComplexOrder
9 
10 variable {I : Type*} [Fintype I] [DecidableEq I]
11 
12 /-- The literal quadratic norm used by the ridge confidence ellipsoid. -/
13 noncomputable def weightedNorm (V : Matrix I I ℝ) (x : I → ℝ) : ℝ :=
14   Real.sqrt (x ⬝ᵥ (V *ᵥ x))
15 
16 omit [DecidableEq I] in
17 theorem weightedNorm_nonneg (V : Matrix I I ℝ) (x : I → ℝ) :
18     0 ≤ weightedNorm V x := Real.sqrt_nonneg _
19 
20 omit [DecidableEq I] in
21 theorem induced_norm_eq {V : Matrix I I ℝ} (hV : V.PosDef) (x : I → ℝ) :
22     @norm (I → ℝ) (V.toNormedAddCommGroup hV).toNorm x = weightedNorm V x := by
23   change Real.sqrt ((V *ᵥ x) ⬝ᵥ x) = _
24   rw [dotProduct_comm]
25   rfl
26 
27 omit [DecidableEq I] in
28 /-- Triangle inequality for the actual positive-definite quadratic norm. -/
29 theorem weightedNorm_add_le {V : Matrix I I ℝ} (hV : V.PosDef) (x y : I → ℝ) :
30     weightedNorm V (x+y) ≤ weightedNorm V x + weightedNorm V y := by
31   rw [← induced_norm_eq hV (x+y), ← induced_norm_eq hV x, ← induced_norm_eq hV y]
32   exact @norm_add_le _ (V.toSeminormedAddCommGroup hV.posSemidef).toSeminormedAddGroup x y
33 
34 omit [DecidableEq I] in
35 theorem weightedNorm_sub_le {V : Matrix I I ℝ} (hV : V.PosDef) (x y : I → ℝ) :
36     weightedNorm V (x-y) ≤ weightedNorm V x + weightedNorm V y := by
37   rw [← induced_norm_eq hV (x-y), ← induced_norm_eq hV x, ← induced_norm_eq hV y]
38   exact @norm_sub_le _ (V.toSeminormedAddCommGroup hV.posSemidef).toSeminormedAddGroup x y
39 
40 /-- Weighted Cauchy--Schwarz uses the genuine inverse matrix norm. -/
41 theorem weighted_cauchy_schwarz {V : Matrix I I ℝ} (hV : V.PosDef) (e x : I → ℝ) :
42     |e ⬝ᵥ x| ≤ weightedNorm V e * weightedNorm V⁻¹ x := by
43   have hinv : V *ᵥ (V⁻¹ *ᵥ x) = x := by
44     rw [mulVec_mulVec, Matrix.mul_nonsing_inv V (isUnit_iff_ne_zero.mpr (ne_of_gt hV.det_pos)),
45       Matrix.one_mulVec]
46   have hquad : weightedNorm V (V⁻¹ *ᵥ x) = weightedNorm V⁻¹ x := by
47     unfold weightedNorm
48     congr 1
49     rw [hinv, dotProduct_comm]
50   have h := @abs_real_inner_le_norm _ (V.toSeminormedAddCommGroup hV.posSemidef)
51     (V.toInnerProductSpace hV.posSemidef) e (V⁻¹ *ᵥ x)
52   change |(V *ᵥ (V⁻¹ *ᵥ x)) ⬝ᵥ e| ≤ _ at h
53   rw [hinv, dotProduct_comm] at h
54   simpa only [induced_norm_eq hV, hquad] using h
55 
56 omit [DecidableEq I] in
57 theorem weightedNorm_sq {V : Matrix I I ℝ} (hV : V.PosSemidef) (x : I → ℝ) :
58     (weightedNorm V x)^2 = x ⬝ᵥ (V *ᵥ x) := by
59   unfold weightedNorm
60   apply Real.sq_sqrt
61   simpa only [star_trivial] using hV.dotProduct_mulVec_nonneg x
62 
63 omit [DecidableEq I] in
64 /-- Euclidean Cauchy--Schwarz in the literal coordinate convention. -/
65 theorem dotProduct_sq_le (x y : I → ℝ) :
66     (x ⬝ᵥ y)^2 ≤ (x ⬝ᵥ x)*(y ⬝ᵥ y) := by
67   simpa [dotProduct,pow_two] using Finset.sum_mul_sq_le_sq_mul_sq Finset.univ x y
68 
69 /-- A genuine quadratic lower bound reverses to the inverse quadratic bound. -/
70 theorem inverse_quadratic_le {V : Matrix I I ℝ} (hV : V.PosDef)
71     (γ : ℝ) (hγ : 0 < γ)
72     (hlower : ∀ y,γ*(y ⬝ᵥ y) ≤ y ⬝ᵥ (V *ᵥ y)) (x : I → ℝ) :
73     γ*(x ⬝ᵥ (V⁻¹ *ᵥ x)) ≤ x ⬝ᵥ x := by
74   let y := V⁻¹ *ᵥ x
75   have hinv : V *ᵥ y = x := by
76     dsimp [y]
77     rw [mulVec_mulVec,Matrix.mul_nonsing_inv V (isUnit_iff_ne_zero.mpr (ne_of_gt hV.det_pos)),
78       Matrix.one_mulVec]
79   have hl := hlower y
80   rw [hinv,dotProduct_comm y x] at hl
81   have hcs := dotProduct_sq_le x y
82   have ha : 0 ≤ x ⬝ᵥ x := by
83     exact Finset.sum_nonneg (fun i _ => mul_self_nonneg (x i))
84   have hr : 0 ≤ x ⬝ᵥ y := by
85     simpa only [y,star_trivial] using hV.inv.posSemidef.dotProduct_mulVec_nonneg x
86   have h1 := mul_le_mul_of_nonneg_left hcs hγ.le
87   have h2 := mul_le_mul_of_nonneg_left hl ha
88   change γ*(x ⬝ᵥ y) ≤ x ⬝ᵥ x
89   by_cases hz : x ⬝ᵥ y = 0
90   · simp [hz,ha]
91   · have hp : 0 < x ⬝ᵥ y := lt_of_le_of_ne hr (Ne.symm hz)
92     nlinarith
93 
94 /-- Euclidean length in coordinates, independent of the function-space sup norm. -/
95 noncomputable def euclideanNorm (x : I → ℝ) : ℝ := weightedNorm (1 : Matrix I I ℝ) x
96 
97 @[simp] theorem euclideanNorm_eq (x : I → ℝ) :
98     euclideanNorm x = Real.sqrt (x ⬝ᵥ x) := by simp [euclideanNorm,weightedNorm]
99 
100 theorem euclideanNorm_sq (x : I → ℝ) :
101     (euclideanNorm x)^2 = x ⬝ᵥ x := by
102   rw [euclideanNorm_eq]
103   exact Real.sq_sqrt (Finset.sum_nonneg (fun i _ => mul_self_nonneg (x i)))
104 
105 omit [DecidableEq I] in
106 theorem weightedNorm_smul (V : Matrix I I ℝ) (c : ℝ) (x : I → ℝ) :
107     weightedNorm V (c • x) = |c| * weightedNorm V x := by
108   unfold weightedNorm
109   have hq : (c • x) ⬝ᵥ (V *ᵥ (c • x)) = c^2*(x ⬝ᵥ (V *ᵥ x)) := by
110     simp only [Matrix.mulVec_smul,dotProduct_smul,smul_dotProduct,smul_eq_mul]
111     ring
112   rw [hq,Real.sqrt_mul (sq_nonneg c),Real.sqrt_sq_eq_abs]
113 
114 theorem weightedNorm_inverse_mulVec {V : Matrix I I ℝ} (hV : V.PosDef) (x : I → ℝ) :
115     weightedNorm V (V⁻¹ *ᵥ x) = weightedNorm V⁻¹ x := by
116   unfold weightedNorm
117   congr 1
118   rw [mulVec_mulVec,Matrix.mul_nonsing_inv V (isUnit_iff_ne_zero.mpr (ne_of_gt hV.det_pos)),
119     Matrix.one_mulVec,dotProduct_comm]
120 
121 /-- The ridge lower bound controls the inverse-matrix norm. -/
122 theorem inverse_weightedNorm_le {V : Matrix I I ℝ} (hV : V.PosDef)
123     (γ : ℝ) (hγ : 0 < γ)
124     (hlower : ∀ y,γ*(y ⬝ᵥ y) ≤ y ⬝ᵥ (V *ᵥ y)) (x : I → ℝ) :
125     weightedNorm V⁻¹ x ≤ euclideanNorm x / Real.sqrt γ := by
126   apply (le_div_iff₀ (Real.sqrt_pos.2 hγ)).2
127   have hs : (weightedNorm V⁻¹ x * Real.sqrt γ)^2 ≤ (euclideanNorm x)^2 := by
128     rw [mul_pow,Real.sq_sqrt hγ.le,weightedNorm_sq hV.inv.posSemidef,euclideanNorm_sq]
129     simpa only [mul_comm] using inverse_quadratic_le hV γ hγ hlower x
130   have hn : 0 ≤ euclideanNorm x := weightedNorm_nonneg _ _
131   nlinarith [mul_nonneg (weightedNorm_nonneg V⁻¹ x) (Real.sqrt_nonneg γ)]
132 
133 /-- The exact regularization contribution has the source's sqrt(gamma) factor. -/
134 theorem ridge_bias_le {V : Matrix I I ℝ} (hV : V.PosDef)
135     (γ : ℝ) (hγ : 0 < γ)
136     (hlower : ∀ y,γ*(y ⬝ᵥ y) ≤ y ⬝ᵥ (V *ᵥ y)) (x : I → ℝ) :
137     weightedNorm V (γ • (V⁻¹ *ᵥ x)) ≤ Real.sqrt γ * euclideanNorm x := by
138   rw [weightedNorm_smul,abs_of_pos hγ,weightedNorm_inverse_mulVec hV]
139   have hq := inverse_quadratic_le hV γ hγ hlower x
140   have hq' := mul_le_mul_of_nonneg_left hq hγ.le
141   have hs : (γ * weightedNorm V⁻¹ x)^2 ≤ (Real.sqrt γ * euclideanNorm x)^2 := by
142     rw [mul_pow,mul_pow,Real.sq_sqrt hγ.le,weightedNorm_sq hV.inv.posSemidef,
143       euclideanNorm_sq]
144     nlinarith
145   have hn : 0 ≤ euclideanNorm x := weightedNorm_nonneg _ _
146   nlinarith [mul_nonneg hγ.le (weightedNorm_nonneg V⁻¹ x),
147     mul_nonneg (Real.sqrt_nonneg γ) hn]
148 
149 /-- The actual accumulated rank-one design matrices. -/
150 noncomputable def designGram (φ : ℕ → I → ℝ) (n : ℕ) : Matrix I I ℝ :=
151   ∑ i∈Finset.range n,Matrix.vecMulVec (φ i) (φ i)
152 
153 /-- The literal ridge matrix in Proposition81. -/
154 noncomputable def ridgeMatrix (γ : ℝ) (φ : ℕ → I → ℝ) (n : ℕ) : Matrix I I ℝ :=
155   γ • (1 : Matrix I I ℝ) + designGram φ n
156 
157 omit [DecidableEq I] in
158 theorem designGram_posSemidef (φ : ℕ → I → ℝ) (n : ℕ) :
159     (designGram φ n).PosSemidef := by
160   apply Matrix.posSemidef_sum
161   intro i _
162   simpa only [star_trivial] using Matrix.posSemidef_vecMulVec_self_star (φ i)
163 
164 theorem ridgeMatrix_posDef (γ : ℝ) (hγ : 0 < γ) (φ : ℕ → I → ℝ) (n : ℕ) :
165     (ridgeMatrix γ φ n).PosDef :=
166   (Matrix.PosDef.one.smul hγ).add_posSemidef (designGram_posSemidef φ n)
167 
168 theorem ridgeMatrix_lower (γ : ℝ) (φ : ℕ → I → ℝ) (n : ℕ) (x : I → ℝ) :
169     γ*(x ⬝ᵥ x) ≤ x ⬝ᵥ (ridgeMatrix γ φ n *ᵥ x) := by
170   have hn : 0 ≤ x ⬝ᵥ (designGram φ n *ᵥ x) := by
171     simpa only [star_trivial] using (designGram_posSemidef φ n).dotProduct_mulVec_nonneg x
172   simpa [ridgeMatrix,Matrix.add_mulVec,Matrix.smul_mulVec,dotProduct_add,dotProduct_smul] using
173     le_add_of_nonneg_right hn (a := γ*(x ⬝ᵥ x))
174 
175 omit [DecidableEq I] in
176 /-- The actual Gram quadratic is the sum of squared predictable projections. -/
177 theorem designGram_quadratic (φ : ℕ → I → ℝ) (n : ℕ) (x : I → ℝ) :
178     x ⬝ᵥ (designGram φ n *ᵥ x) = ∑ i∈Finset.range n,(x ⬝ᵥ φ i)^2 := by
179   unfold designGram
180   rw [Matrix.sum_mulVec,dotProduct_sum]
181   apply Finset.sum_congr rfl
182   intro i _
183   have hout : Matrix.vecMulVec (φ i) (φ i) *ᵥ x = (x ⬝ᵥ φ i) • φ i := by
184     ext j
185     simp [Matrix.mulVec,dotProduct,Matrix.vecMulVec_apply,Finset.mul_sum,mul_comm,mul_left_comm]
186   rw [hout,dotProduct_smul,smul_eq_mul,pow_two]
187 
188 /-- Literal design-times-response sum; it also defines the accumulated noise vector. -/
189 noncomputable def responseVector (φ : ℕ → I → ℝ) (y : ℕ → ℝ) (n : ℕ) : I → ℝ :=
190   ∑ i∈Finset.range n,y i • φ i
191 
192 /-- The actual ridge estimator obtained by matrix inversion. -/
193 noncomputable def ridgeEstimate (γ : ℝ) (φ : ℕ → I → ℝ) (y : ℕ → ℝ) (n : ℕ) : I → ℝ :=
194   (ridgeMatrix γ φ n)⁻¹ *ᵥ responseVector φ y n
195 
196 omit [DecidableEq I] in
197 theorem responseVector_regression (φ : ℕ → I → ℝ) (y ξ : ℕ → ℝ) (θ : I → ℝ)
198     (n : ℕ) (hy : ∀ i<n,y i=θ ⬝ᵥ φ i+ξ i) :
199     responseVector φ y n = designGram φ n *ᵥ θ+responseVector φ ξ n := by
200   have hout (i : ℕ) : Matrix.vecMulVec (φ i) (φ i) *ᵥ θ = (θ ⬝ᵥ φ i) • φ i := by
201     ext j
202     simp [Matrix.mulVec,dotProduct,Matrix.vecMulVec_apply,Finset.mul_sum,mul_comm,mul_left_comm]
203   unfold responseVector designGram
204   rw [Matrix.sum_mulVec,← Finset.sum_add_distrib]
205   apply Finset.sum_congr rfl
206   intro i hi
207   rw [hy i (Finset.mem_range.mp hi),add_smul,hout]
208 
209 /-- The regression equation gives the exact ridge error identity before any
210 probabilistic bound is applied. -/
211 theorem ridge_error_identity (γ : ℝ) (hγ : 0 < γ) (φ : ℕ → I → ℝ)
212     (y ξ : ℕ → ℝ) (θ : I → ℝ) (n : ℕ)
213     (hy : ∀ i<n,y i=θ ⬝ᵥ φ i+ξ i) :
214     ridgeEstimate γ φ y n-θ =
215       (ridgeMatrix γ φ n)⁻¹ *ᵥ (responseVector φ ξ n-γ • θ) := by
216   let V := ridgeMatrix γ φ n
217   have hV := ridgeMatrix_posDef γ hγ φ n
218   have hr : responseVector φ y n = V *ᵥ θ+(responseVector φ ξ n-γ • θ) := by
219     rw [responseVector_regression φ y ξ θ n hy]
220     simp only [V,ridgeMatrix,Matrix.add_mulVec,Matrix.smul_mulVec,Matrix.one_mulVec]
221     abel
222   have hinv : V⁻¹ *ᵥ (V *ᵥ θ)=θ := by
223     rw [mulVec_mulVec,Matrix.nonsing_inv_mul V (isUnit_iff_ne_zero.mpr (ne_of_gt hV.det_pos)),
224       Matrix.one_mulVec]
225   change V⁻¹ *ᵥ responseVector φ y n-θ=V⁻¹ *ᵥ (responseVector φ ξ n-γ • θ)
226   rw [hr,Matrix.mulVec_add,hinv]
227   abel
228 
229 /-- The deterministic noise/regularization decomposition yields the true
230 ellipsoid error bound, before the anytime noise event is supplied. -/
231 theorem ridge_error_le (γ : ℝ) (hγ : 0 < γ) (φ : ℕ → I → ℝ)
232     (y ξ : ℕ → ℝ) (θ : I → ℝ) (n : ℕ)
233     (hy : ∀ i<n,y i=θ ⬝ᵥ φ i+ξ i) :
234     weightedNorm (ridgeMatrix γ φ n) (ridgeEstimate γ φ y n-θ) ≤
235       weightedNorm (ridgeMatrix γ φ n)⁻¹ (responseVector φ ξ n)+
236         Real.sqrt γ * euclideanNorm θ := by
237   have hV := ridgeMatrix_posDef γ hγ φ n
238   rw [ridge_error_identity γ hγ φ y ξ θ n hy,Matrix.mulVec_sub,Matrix.mulVec_smul]
239   apply (weightedNorm_sub_le hV _ _).trans
240   rw [weightedNorm_inverse_mulVec hV]
241   exact add_le_add le_rfl (ridge_bias_le hV γ hγ (ridgeMatrix_lower γ φ n) θ)
242 
243 end Komorebi.LSM.RegressionConfidence

Back to source controls