RegularizedHilbertRepair.lean

1 import Mathlib.Analysis.InnerProductSpace.LaxMilgram
2 import Mathlib.Analysis.InnerProductSpace.Adjoint
3 import Mathlib.Tactic.Linarith
4 import Mathlib.Tactic.Ring
5 import Mathlib.Tactic.Abel
6 
7 set_option autoImplicit false
8 
9 noncomputable section
10 universe u v
11 namespace KomorebiSpace.RegularizedHilbertRepair
12 open scoped RealInnerProductSpace
13 variable {H : Type u} {Y : Type v}
14 variable [NormedAddCommGroup H] [InnerProductSpace ℝ H] [CompleteSpace H]
15 variable [NormedAddCommGroup Y] [InnerProductSpace ℝ Y] [CompleteSpace Y]
16 variable (T : H →L[ℝ] Y) (α : ℝ)
17 
18 def normal : H →L[ℝ] H := T.adjoint.comp T + α • ContinuousLinearMap.id ℝ H
19 
20 @[simp] theorem normal_apply (a : H) : normal T α a = T.adjoint (T a) + α • a := rfl
21 
22 @[simp] theorem normal_inner (a b : H) :
23     ⟪normal T α a, b⟫ = ⟪T a, T b⟫ + α * ⟪a, b⟫ := by
24   simp only [normal_apply, inner_add_left, real_inner_smul_left,
25     ContinuousLinearMap.adjoint_inner_left]
26 
27 def energyBilin : H →L[ℝ] H →L[ℝ] ℝ := (innerSL ℝ).comp (normal T α)
28 
29 @[simp] theorem energyBilin_apply (a b : H) :
30     energyBilin T α a b = ⟪T a, T b⟫ + α * ⟪a, b⟫ := by
31   change ⟪normal T α a, b⟫ = _
32   exact normal_inner T α a b
33 
34 theorem energyBilin_coercive (hα : 0 < α) : IsCoercive (energyBilin T α) := by
35   refine ⟨α, hα, ?_⟩
36   intro a
37   rw [energyBilin_apply, real_inner_self_eq_norm_sq, real_inner_self_eq_norm_sq]
38   nlinarith [sq_nonneg ‖T a‖]
39 
40 def normalEquiv (hα : 0 < α) : H ≃L[ℝ] H :=
41   (energyBilin_coercive T α hα).continuousLinearEquivOfBilin
42 
43 @[simp] theorem normalEquiv_apply (hα : 0 < α) (a : H) :
44     normalEquiv T α hα a = normal T α a := by
45   apply ext_inner_right ℝ
46   intro b
47   exact (IsCoercive.continuousLinearEquivOfBilin_apply (energyBilin_coercive T α hα) a b).trans
48     (show energyBilin T α a b = ⟪normal T α a,b⟫ from rfl)
49 
50 def correction (hα : 0 < α) (r : Y) : H :=
51   (normalEquiv T α hα).symm (T.adjoint r)
52 
53 @[simp] theorem normal_correction (hα : 0 < α) (r : Y) :
54     normal T α (correction T α hα r) = T.adjoint r := by
55   rw [← normalEquiv_apply T α hα]
56   exact (normalEquiv T α hα).apply_symm_apply _
57 
58 def objective (r : Y) (a : H) : ℝ := ‖T a - r‖^2 + α * ‖a‖^2
59 
60 def gradient (r : Y) (a : H) : H := normal T α a - T.adjoint r
61 
62 @[simp] theorem gradient_correction (hα : 0 < α) (r : Y) :
63     gradient T α r (correction T α hα r) = 0 := by
64   rw [gradient, normal_correction, sub_self]
65 
66 theorem normal_equation (hα : 0 < α) (r : Y) :
67     T.adjoint (r - T (correction T α hα r)) = α • correction T α hα r := by
68   have h := normal_correction T α hα r
69   rw [normal_apply] at h
70   rw [map_sub, ← h]
71   abel
72 
73 theorem objective_expansion (r : Y) (a h : H) :
74     objective T α r (a+h) = objective T α r a +
75       2 * ⟪gradient T α r a,h⟫ + ‖T h‖^2 + α * ‖h‖^2 := by
76   have hm : T (a+h)-r = (T a-r)+T h := by rw [map_add]; abel
77   simp only [objective, hm, norm_add_sq_real, gradient, inner_sub_left,
78     normal_inner, ContinuousLinearMap.adjoint_inner_left, inner_sub_left]
79   ring
80 
81 theorem energy_identity (hα : 0 < α) (r : Y) (a : H) :
82     objective T α r a - objective T α r (correction T α hα r) =
83       ‖T (a-correction T α hα r)‖^2 + α * ‖a-correction T α hα r‖^2 := by
84   have he := objective_expansion T α r (correction T α hα r) (a-correction T α hα r)
85   rw [add_sub_cancel, gradient_correction, inner_zero_left] at he
86   linarith
87 
88 /-- The constructed minimizer is global, not only a stationary point. -/
89 theorem correction_best (hα : 0 < α) (r : Y) (a : H) :
90     objective T α r (correction T α hα r) ≤ objective T α r a := by
91   have he := energy_identity T α hα r a
92   have hn := mul_nonneg hα.le (sq_nonneg ‖a-correction T α hα r‖)
93   nlinarith [sq_nonneg ‖T (a-correction T α hα r)‖]
94 
95 theorem correction_unique (hα : 0 < α) (r : Y) (a : H)
96     (ha : objective T α r a ≤ objective T α r (correction T α hα r)) :
97     a = correction T α hα r := by
98   have he := energy_identity T α hα r a
99   have hz : ‖a-correction T α hα r‖^2 = 0 := by
100     nlinarith [sq_nonneg ‖T (a-correction T α hα r)‖,
101       sq_nonneg ‖a-correction T α hα r‖]
102   exact sub_eq_zero.mp (norm_eq_zero.mp (sq_eq_zero_iff.mp hz))
103 
104 theorem exists_unique_minimizer (hα : 0 < α) (r : Y) :
105     ∃! a : H, ∀ b : H, objective T α r a ≤ objective T α r b := by
106   refine ⟨correction T α hα r, correction_best T α hα r, ?_⟩
107   intro a ha
108   exact correction_unique T α hα r a (ha _)
109 
110 theorem normal_unique (hα : 0 < α) (r : Y) (a : H)
111     (ha : normal T α a = T.adjoint r) : a = correction T α hα r := by
112   apply (normalEquiv T α hα).injective
113   rw [normalEquiv_apply, normalEquiv_apply, normal_correction]
114   exact ha
115 
116 theorem gradient_zero_iff (hα : 0 < α) (r : Y) (a : H) :
117     gradient T α r a = 0 ↔ a = correction T α hα r := by
118   constructor
119   · intro ha
120     exact normal_unique T α hα r a (sub_eq_zero.mp ha)
121   · rintro rfl
122     exact gradient_correction T α hα r
123 
124 @[simp] theorem correction_zero (hα : 0 < α) : correction T α hα 0 = 0 := by
125   symm
126   apply normal_unique
127   simp
128 
129 theorem gradient_error (hα : 0 < α) (r : Y) (a : H) :
130     gradient T α r a = normal T α (a-correction T α hα r) := by
131   rw [map_sub, normal_correction]
132   rfl
133 
134 theorem normal_bound (a : H) :
135     α * ‖a‖ ≤ ‖normal T α a‖ := by
136   by_cases hz : a = 0
137   · simp [hz]
138   · have hp : 0 < ‖a‖ := norm_pos_iff.mpr hz
139     have hi := real_inner_le_norm (normal T α a) a
140     rw [normal_inner, real_inner_self_eq_norm_sq, real_inner_self_eq_norm_sq] at hi
141     nlinarith [sq_nonneg ‖T a‖]
142 
143 /-- Actual normal-equation error certifies correction error at the same parameter. -/
144 theorem correction_error_bound (hα : 0 < α) (r : Y) (a : H) :
145     ‖a-correction T α hα r‖ ≤ ‖gradient T α r a‖ / α := by
146   apply (le_div_iff₀ hα).2
147   rw [gradient_error]
148   simpa [mul_comm] using normal_bound T α (a-correction T α hα r)
149 
150 theorem gap_as_inner (hα : 0 < α) (r : Y) (a : H) :
151     objective T α r a - objective T α r (correction T α hα r) =
152       ⟪gradient T α r a, a-correction T α hα r⟫ := by
153   rw [energy_identity, gradient_error, normal_inner,
154     real_inner_self_eq_norm_sq, real_inner_self_eq_norm_sq]
155 
156 theorem objective_gap_bound (hα : 0 < α) (r : Y) (a : H) :
157     objective T α r a - objective T α r (correction T α hα r) ≤
158       ‖gradient T α r a‖^2 / α := by
159   rw [gap_as_inner]
160   calc
161     ⟪gradient T α r a, a-correction T α hα r⟫ ≤
162         ‖gradient T α r a‖ * ‖a-correction T α hα r‖ := real_inner_le_norm _ _
163     _ ≤ ‖gradient T α r a‖ * (‖gradient T α r a‖/α) :=
164       mul_le_mul_of_nonneg_left (correction_error_bound T α hα r a) (norm_nonneg _)
165     _ = ‖gradient T α r a‖^2/α := by ring
166 
167 /-- A norm budget does not silently become a penalty on protected data. -/
168 theorem residual_best_at_budget (hα : 0 < α) (r : Y) (a : H)
169     (ha : ‖a‖ ≤ ‖correction T α hα r‖) :
170     ‖T (correction T α hα r)-r‖ ≤ ‖T a-r‖ := by
171   have hb := correction_best T α hα r a
172   unfold objective at hb
173   have hm : ‖a‖^2 ≤ ‖correction T α hα r‖^2 := by
174     nlinarith [norm_nonneg a, norm_nonneg (correction T α hα r)]
175   have hp := mul_le_mul_of_nonneg_left hm hα.le
176   nlinarith [norm_nonneg (T (correction T α hα r)-r), norm_nonneg (T a-r)]
177 
178 theorem parameter_norm_antitone (hα : 0 < α) (γ : ℝ) (hγ : 0 < γ) (h : α < γ)
179     (r : Y) : ‖correction T γ hγ r‖ ≤ ‖correction T α hα r‖ := by
180   have h1 := correction_best T α hα r (correction T γ hγ r)
181   have h2 := correction_best T γ hγ r (correction T α hα r)
182   unfold objective at h1 h2
183   have hs : ‖correction T γ hγ r‖^2 ≤ ‖correction T α hα r‖^2 := by nlinarith
184   nlinarith [norm_nonneg (correction T γ hγ r), norm_nonneg (correction T α hα r)]
185 
186 theorem parameter_residual_monotone (hα : 0 < α) (γ : ℝ) (hγ : 0 < γ) (h : α < γ)
187     (r : Y) : ‖T (correction T α hα r)-r‖ ≤ ‖T (correction T γ hγ r)-r‖ := by
188   exact residual_best_at_budget T α hα r _ (parameter_norm_antitone T α hα γ hγ h r)
189 
190 theorem correction_data_energy (hα : 0 < α) (r s : Y) :
191     ‖T (correction T α hα r-correction T α hα s)‖^2 +
192       α * ‖correction T α hα r-correction T α hα s‖^2 =
193       ⟪r-s,T (correction T α hα r-correction T α hα s)⟫ := by
194   have hn : normal T α (correction T α hα r-correction T α hα s) =
195       T.adjoint (r-s) := by rw [map_sub, normal_correction, normal_correction, map_sub]
196   have hi := congrArg (fun z => ⟪z,correction T α hα r-correction T α hα s⟫) hn
197   simpa only [normal_inner, real_inner_self_eq_norm_sq,
198     ContinuousLinearMap.adjoint_inner_left] using hi
199 
200 theorem data_noise_squared (hα : 0 < α) (r s : Y) :
201     4*α*‖correction T α hα r-correction T α hα s‖^2 ≤ ‖r-s‖^2 := by
202   have he := correction_data_energy T α hα r s
203   have hi := real_inner_le_norm (r-s) (T (correction T α hα r-correction T α hα s))
204   nlinarith [sq_nonneg (2*‖T (correction T α hα r-correction T α hα s)‖-‖r-s‖)]
205 
206 theorem data_noise_bound (hα : 0 < α) (r s : Y) :
207     ‖correction T α hα r-correction T α hα s‖ ≤ ‖r-s‖/(2*Real.sqrt α) := by
208   have hsq := data_noise_squared T α hα r s
209   have hs := Real.sq_sqrt hα.le
210   have hp := Real.sqrt_pos.2 hα
211   apply (le_div_iff₀ (by positivity : 0 < 2*Real.sqrt α)).2
212   have hn := norm_nonneg (correction T α hα r-correction T α hα s)
213   have hd := norm_nonneg (r-s)
214   nlinarith [sq_nonneg (‖r-s‖-2*Real.sqrt α*‖correction T α hα r-correction T α hα s‖)]
215 
216 theorem data_fitted_bound (hα : 0 < α) (r s : Y) :
217     ‖T (correction T α hα r-correction T α hα s)‖ ≤ ‖r-s‖ := by
218   have he := correction_data_energy T α hα r s
219   have hi := real_inner_le_norm (r-s) (T (correction T α hα r-correction T α hα s))
220   have hn := mul_nonneg hα.le (sq_nonneg ‖correction T α hα r-correction T α hα s‖)
221   nlinarith [norm_nonneg (T (correction T α hα r-correction T α hα s)), norm_nonneg (r-s)]
222 
223 /-- Solver error in the fitted data has the sharp universal square-root gain. -/
224 theorem fitted_error_squared (hα : 0 < α) (r : Y) (a : H) :
225     4*α*‖T (a-correction T α hα r)‖^2 ≤ ‖gradient T α r a‖^2 := by
226   have he := gap_as_inner T α hα r a
227   rw [energy_identity] at he
228   have hi := real_inner_le_norm (gradient T α r a) (a-correction T α hα r)
229   nlinarith [sq_nonneg (2*α*‖a-correction T α hα r‖-‖gradient T α r a‖)]
230 
231 theorem fitted_error_bound (hα : 0 < α) (r : Y) (a : H) :
232     ‖T (a-correction T α hα r)‖ ≤ ‖gradient T α r a‖/(2*Real.sqrt α) := by
233   have hsq := fitted_error_squared T α hα r a
234   have hs := Real.sq_sqrt hα.le
235   have hp := Real.sqrt_pos.2 hα
236   apply (le_div_iff₀ (by positivity : 0 < 2*Real.sqrt α)).2
237   have hn := norm_nonneg (T (a-correction T α hα r))
238   have hd := norm_nonneg (gradient T α r a)
239   nlinarith [sq_nonneg (‖gradient T α r a‖-2*Real.sqrt α*‖T (a-correction T α hα r)‖)]
240 
241 /-- Numerical solver error and measured-data error are separate terms. -/
242 theorem solver_and_data_bound (hα : 0 < α) (r s : Y) (a : H) :
243     ‖a-correction T α hα r‖ ≤
244       ‖gradient T α s a‖/α + ‖s-r‖/(2*Real.sqrt α) := by
245   have he : a-correction T α hα r =
246       (a-correction T α hα s)+(correction T α hα s-correction T α hα r) := by abel
247   calc
248     ‖a-correction T α hα r‖ =
249         ‖(a-correction T α hα s)+(correction T α hα s-correction T α hα r)‖ := by rw [he]
250     _ ≤ ‖a-correction T α hα s‖+‖correction T α hα s-correction T α hα r‖ := norm_add_le _ _
251     _ ≤ ‖gradient T α s a‖/α + ‖s-r‖/(2*Real.sqrt α) :=
252       add_le_add (correction_error_bound T α hα s a) (data_noise_bound T α hα s r)
253 
254 end KomorebiSpace.RegularizedHilbertRepair

Back to source controls