CoordinateCovector.lean

1 import KomorebiSpace.CoordinateBasisChange
2 import KomorebiSpace.ProtectedRegularizedRepair
3 import Mathlib.Analysis.InnerProductSpace.PiL2
4 import Mathlib.LinearAlgebra.Matrix.PosDef
5 import Mathlib.LinearAlgebra.Matrix.ToLin
6 
7 set_option autoImplicit false
8 
9 /-!
10 # Coordinate covectors in an arbitrary full Hilbert basis
11 
12 The basis is a basis of the entire correction space.  Coefficients and
13 covectors are finite real coordinate functions; every norm in the certificate
14 is the original Hilbert norm.  No nonemptiness assumption is made on the basis
15 index, so the zero-dimensional correction space is included.
16 -/
17 
18 noncomputable section
19 universe u v w
20 namespace KomorebiSpace.CoordinateCovector
21 open Module
22 open scoped RealInnerProductSpace Matrix
23 
24 variable {ι : Type u} [Fintype ι]
25 variable {H : Type v} [NormedAddCommGroup H] [InnerProductSpace ℝ H]
26 
27 /-- Synthesis in a full, possibly nonorthonormal basis. -/
28 def synthesis (B : Basis ι ℝ H) (c : ι → ℝ) : H := B.equivFun.symm c
29 
30 /-- The coordinate covector of a vector consists of its pairings with the basis. -/
31 def pairings (B : Basis ι ℝ H) (x : H) : ι → ℝ := fun i => ⟪B i, x⟫
32 
33 /-- The Gram matrix uses the actual Hilbert inner product. -/
34 def gram (B : Basis ι ℝ H) : Matrix ι ι ℝ := fun i j => ⟪B i, B j⟫
35 
36 theorem synthesis_eq_sum (B : Basis ι ℝ H) (c : ι → ℝ) :
37     synthesis B c = ∑ i, c i • B i := by
38   exact B.equivFun_symm_apply c
39 
40 theorem synthesis_injective (B : Basis ι ℝ H) : Function.Injective (synthesis B) :=
41   B.equivFun.symm.injective
42 
43 theorem pairings_synthesis (B : Basis ι ℝ H) (c : ι → ℝ) :
44     pairings B (synthesis B c) = (gram B) *ᵥ c := by
45   ext i
46   simp [pairings, synthesis_eq_sum, gram, Matrix.mulVec, dotProduct, inner_sum,
47     real_inner_smul_right, mul_comm]
48 
49 theorem inner_synthesis (B : Basis ι ℝ H) (c : ι → ℝ) (x : H) :
50     ⟪synthesis B c, x⟫ = dotProduct c (pairings B x) := by
51   simp [synthesis_eq_sum, pairings, dotProduct, sum_inner, real_inner_smul_left]
52 
53 theorem gram_posDef (B : Basis ι ℝ H) : (gram B).PosDef := by
54   classical
55   refine Matrix.PosDef.of_dotProduct_mulVec_pos ?_ ?_
56   · ext i j
57     simp [Matrix.conjTranspose_apply, gram, real_inner_comm]
58   · intro c hc
59     have hz : synthesis B c ≠ 0 := by
60       intro h
61       apply hc
62       exact B.equivFun.symm.injective (h.trans (map_zero B.equivFun.symm).symm)
63     have hp := real_inner_self_pos.mpr hz
64     rw [inner_synthesis, pairings_synthesis] at hp
65     simpa using hp
66 
67 variable [DecidableEq ι]
68 
69 theorem gram_isUnit (B : Basis ι ℝ H) : IsUnit (gram B) := (gram_posDef B).isUnit
70 
71 theorem inverse_gram_pairings (B : Basis ι ℝ H) (x : H) :
72     (gram B)⁻¹ *ᵥ pairings B x = B.equivFun x := by
73   have h : pairings B x = gram B *ᵥ B.equivFun x := by
74     rw [← pairings_synthesis]
75     simp [synthesis]
76   rw [h, Matrix.mulVec_mulVec,
77     Matrix.nonsing_inv_mul _ ((gram B).isUnit_iff_isUnit_det.mp (gram_isUnit B)),
78     Matrix.one_mulVec]
79 
80 theorem reconstruct (B : Basis ι ℝ H) (x : H) :
81     x = synthesis B ((gram B)⁻¹ *ᵥ pairings B x) := by
82   rw [inverse_gram_pairings]
83   simp [synthesis]
84 
85 theorem norm_sq_pairings (B : Basis ι ℝ H) (x : H) :
86     ‖x‖ ^ 2 = dotProduct (pairings B x) ((gram B)⁻¹ *ᵥ pairings B x) := by
87   rw [← real_inner_self_eq_norm_sq]
88   conv_lhs => lhs; rw [reconstruct B x]
89   rw [inner_synthesis, dotProduct_comm]
90 
91 omit [DecidableEq ι] in
92 /-- The coordinate quadratic uses the Euclidean inner product, never the Pi norm. -/
93 theorem dotProduct_eq_euclidean_inner (c d : ι → ℝ) :
94     dotProduct c d = inner ℝ (WithLp.toLp 2 c : EuclideanSpace ℝ ι)
95       (WithLp.toLp 2 d) := by
96   simp [PiLp.inner_apply, dotProduct, mul_comm]
97 
98 /-- An empty full basis gives the zero scalar certificate in the original norm. -/
99 theorem norm_sq_eq_zero_of_empty_basis [IsEmpty ι] (B : Basis ι ℝ H) (x : H) :
100     ‖x‖ ^ 2 = 0 := by
101   rw [norm_sq_pairings B]
102   simp [dotProduct]
103 
104 /-- Pairings with the changed full basis transform covariantly. -/
105 theorem pairings_changedBasis (B : Basis ι ℝ H) (C : Matrix ι ι ℝ)
106     (hC : IsUnit C) (x : H) :
107     pairings (changedBasis B C hC) x = Cᵀ *ᵥ pairings B x := by
108   ext i
109   change ⟪changedBasis B C hC i, x⟫ = _
110   rw [changedBasis_apply]
111   change ⟪synthesis B (fun j => C j i), x⟫ = _
112   rw [inner_synthesis]
113   rfl
114 
115 /-- The Gram of the actual changed full basis is the matrix congruence CᵀGC. -/
116 theorem gram_changedBasis (B : Basis ι ℝ H) (C : Matrix ι ι ℝ)
117     (hC : IsUnit C) :
118     gram (changedBasis B C hC) = Cᵀ * gram B * C := by
119   ext i j
120   change ⟪changedBasis B C hC i, changedBasis B C hC j⟫ = _
121   rw [changedBasis_apply, changedBasis_apply]
122   change ⟪synthesis B (fun k => C k i), synthesis B (fun k => C k j)⟫ = _
123   rw [inner_synthesis, pairings_synthesis]
124   simp only [dotProduct, Matrix.mulVec, Matrix.mul_apply, Matrix.transpose_apply]
125   simp_rw [Finset.mul_sum, Finset.sum_mul]
126   rw [Finset.sum_comm]
127   apply Finset.sum_congr rfl
128   intro k _
129   apply Finset.sum_congr rfl
130   intro l _
131   ring
132 
133 section OrthonormalCoordinates
134 variable {κ : Type w} [Fintype κ]
135 
136 /-- Coordinates in an orthonormal basis, used only for the residual space. -/
137 def orthonormalCoordinates {Y : Type*} [NormedAddCommGroup Y]
138     [InnerProductSpace ℝ Y] (e : OrthonormalBasis κ ℝ Y) (y : Y) : κ → ℝ :=
139   fun j => ⟪e j, y⟫
140 
141 /-- Columns of an arbitrary basis in orthonormal Hilbert coordinates. -/
142 def basisMatrix (e : OrthonormalBasis κ ℝ H) (B : Basis ι ℝ H) : Matrix κ ι ℝ :=
143   fun j i => ⟪e j, B i⟫
144 
145 omit [Fintype ι] [DecidableEq ι] in
146 theorem gram_eq_transpose_mul (e : OrthonormalBasis κ ℝ H) (B : Basis ι ℝ H) :
147     gram B = (basisMatrix e B)ᵀ * basisMatrix e B := by
148   ext i j
149   change ⟪B i, B j⟫ = ∑ k, ⟪e k, B i⟫ * ⟪e k, B j⟫
150   rw [← e.sum_inner_mul_inner (B i) (B j)]
151   apply Finset.sum_congr rfl
152   intro k _
153   rw [real_inner_comm (B i) (e k)]
154 
155 variable {Y : Type*} [NormedAddCommGroup Y] [InnerProductSpace ℝ Y]
156 
157 /-- The literal observation matrix has columns equal to the observed basis vectors. -/
158 def operatorMatrix (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
159     (T : H →L[ℝ] Y) : Matrix κ ι ℝ := fun j i => ⟪e j, T (B i)⟫
160 
161 omit [DecidableEq ι] in
162 theorem operatorMatrix_synthesis (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
163     (T : H →L[ℝ] Y) (c : ι → ℝ) :
164     orthonormalCoordinates e (T (synthesis B c)) = operatorMatrix B e T *ᵥ c := by
165   ext j
166   simp [orthonormalCoordinates, synthesis_eq_sum, operatorMatrix, map_sum,
167     Matrix.mulVec, dotProduct, inner_sum, real_inner_smul_right, mul_comm]
168 
169 /-- The observation matrix of the changed full basis is SC. -/
170 theorem operatorMatrix_changedBasis (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
171     (T : H →L[ℝ] Y) (C : Matrix ι ι ℝ) (hC : IsUnit C) :
172     operatorMatrix (changedBasis B C hC) e T = operatorMatrix B e T * C := by
173   ext j i
174   change ⟪e j, T (changedBasis B C hC i)⟫ = _
175   rw [changedBasis_apply]
176   exact congrFun (operatorMatrix_synthesis B e T (fun k => C k i)) j
177 
178 /-- The computed residual in the changed basis is exactly Cᵀq. -/
179 theorem computed_normalResidual_changedBasis
180     (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y) (T : H →L[ℝ] Y)
181     (C : Matrix ι ι ℝ) (hC : IsUnit C) (α : ℝ) (c : ι → ℝ) (r : Y) :
182     normalResidual (operatorMatrix (changedBasis B C hC) e T)
183         (gram (changedBasis B C hC)) α (C⁻¹ *ᵥ c) (orthonormalCoordinates e r) =
184       Cᵀ *ᵥ normalResidual (operatorMatrix B e T) (gram B) α c
185         (orthonormalCoordinates e r) := by
186   rw [operatorMatrix_changedBasis, gram_changedBasis]
187   exact normalResidual_basis_change _ _ _ hC _ _ _
188 
189 /-- WR15 invariance for two actual full bases of the same Hilbert space. -/
190 theorem computed_certificate_changedBasis
191     (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y) (T : H →L[ℝ] Y)
192     (C : Matrix ι ι ℝ) (hC : IsUnit C) (α : ℝ) (c : ι → ℝ) (r : Y) :
193     let B' := changedBasis B C hC
194     let q' := normalResidual (operatorMatrix B' e T) (gram B') α (C⁻¹ *ᵥ c)
195       (orthonormalCoordinates e r)
196     let q := normalResidual (operatorMatrix B e T) (gram B) α c
197       (orthonormalCoordinates e r)
198     dotProduct q' ((gram B')⁻¹ *ᵥ q') = dotProduct q ((gram B)⁻¹ *ᵥ q) := by
199   dsimp only
200   rw [operatorMatrix_changedBasis, gram_changedBasis]
201   exact normalResidual_certificate_basis_change _ _ _ hC _ _ _
202 
203 variable [CompleteSpace H] [CompleteSpace Y]
204 
205 omit [Fintype ι] [DecidableEq ι] in
206 theorem pairings_adjoint (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
207     (T : H →L[ℝ] Y) (y : Y) :
208     pairings B (T.adjoint y) = (operatorMatrix B e T)ᵀ *ᵥ orthonormalCoordinates e y := by
209   ext i
210   change ⟪B i, T.adjoint y⟫ = ∑ j, ⟪e j, T (B i)⟫ * ⟪e j, y⟫
211   rw [ContinuousLinearMap.adjoint_inner_right, ← e.sum_inner_mul_inner (T (B i)) y]
212   apply Finset.sum_congr rfl
213   intro j _
214   rw [real_inner_comm (T (B i)) (e j)]
215 
216 omit [DecidableEq ι] in
217 /-- The displayed normal-equation residual is derived to be the covector of the
218 accepted Hilbert normal residual; this is not an additional assumption. -/
219 theorem normalResidual_eq_pairings (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
220     (T : H →L[ℝ] Y) (α : ℝ) (c : ι → ℝ) (r : Y) :
221     normalResidual (operatorMatrix B e T) (gram B) α c (orthonormalCoordinates e r) =
222       pairings B (RegularizedHilbertRepair.gradient T α r (synthesis B c)) := by
223   have hT := pairings_adjoint B e T (T (synthesis B c))
224   rw [operatorMatrix_synthesis] at hT
225   have hr := pairings_adjoint B e T r
226   have hG := pairings_synthesis B c
227   ext i
228   simp only [normalResidual, Matrix.add_mulVec, Matrix.smul_mulVec,
229     ← Matrix.mulVec_mulVec, Pi.sub_apply, Pi.add_apply,
230     Pi.smul_apply, smul_eq_mul]
231   change _ = ⟪B i, T.adjoint (T (synthesis B c)) + α • synthesis B c - T.adjoint r⟫
232   rw [inner_sub_right, inner_add_right, real_inner_smul_right]
233   rw [← congrFun hT i, ← congrFun hr i, ← congrFun hG i]
234   rfl
235 
236 /-- WR15: the actual normal residual is obtained by the inverse Gram metric. -/
237 theorem gradient_eq_synthesis_inverse_normalResidual
238     (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
239     (T : H →L[ℝ] Y) (α : ℝ) (c : ι → ℝ) (r : Y) :
240     RegularizedHilbertRepair.gradient T α r (synthesis B c) =
241       synthesis B ((gram B)⁻¹ *ᵥ
242         normalResidual (operatorMatrix B e T) (gram B) α c (orthonormalCoordinates e r)) := by
243   rw [normalResidual_eq_pairings]
244   exact reconstruct B _
245 
246 /-- WR15: the Gram-corrected scalar is the squared original Hilbert norm. -/
247 theorem gradient_norm_sq_eq_normalResidual
248     (B : Basis ι ℝ H) (e : OrthonormalBasis κ ℝ Y)
249     (T : H →L[ℝ] Y) (α : ℝ) (c : ι → ℝ) (r : Y) :
250     ‖RegularizedHilbertRepair.gradient T α r (synthesis B c)‖ ^ 2 =
251       dotProduct
252         (normalResidual (operatorMatrix B e T) (gram B) α c (orthonormalCoordinates e r))
253         ((gram B)⁻¹ *ᵥ
254           normalResidual (operatorMatrix B e T) (gram B) α c (orthonormalCoordinates e r)) := by
255   rw [normalResidual_eq_pairings]
256   exact norm_sq_pairings B _
257 
258 end OrthonormalCoordinates
259 
260 section ProtectedKernel
261 variable {E Y Z : Type*}
262 variable [NormedAddCommGroup E] [InnerProductSpace ℝ E] [CompleteSpace E]
263 variable [NormedAddCommGroup Y] [InnerProductSpace ℝ Y] [CompleteSpace Y]
264 variable [NormedAddCommGroup Z] [NormedSpace ℝ Z]
265 variable {κ η : Type*} [Fintype κ] [Fintype η]
266 variable (D : E →L[ℝ] Y) (p : E →L[ℝ] Z)
267 
268 omit [Fintype ι] [DecidableEq ι] in
269 /-- The accepted gradient is the projected ambient normal residual on the whole kernel. -/
270 theorem protected_gradient_eq_projected (α : ℝ) (r : Y) (a : p.ker) :
271     RegularizedHilbertRepair.gradient (ProtectedHilbertRepair.restricted D p) α r a =
272       p.ker.orthogonalProjectionOnto (D.adjoint (D a.val - r)) + α • a := by
273   rw [← ProtectedRegularizedRepair.restricted_adjoint, map_sub]
274   unfold RegularizedHilbertRepair.gradient
275   rw [RegularizedHilbertRepair.normal_apply]
276   change _ = ((ProtectedHilbertRepair.restricted D p).adjoint
277     ((ProtectedHilbertRepair.restricted D p) a) -
278     (ProtectedHilbertRepair.restricted D p).adjoint r) + α • a
279   abel
280 
281 omit [Fintype ι] [DecidableEq ι] [CompleteSpace E] in
282 /-- The arbitrary full kernel basis has the literal ambient Gram matrix BᵀB. -/
283 theorem protected_gram_eq_transpose_mul (B : Basis ι ℝ p.ker)
284     (eE : OrthonormalBasis η ℝ E) :
285     gram B = (operatorMatrix B eE p.ker.subtypeL)ᵀ *
286       operatorMatrix B eE p.ker.subtypeL := by
287   ext i j
288   change ⟪(B i).val, (B j).val⟫ =
289     ∑ k, ⟪eE k, (B i).val⟫ * ⟪eE k, (B j).val⟫
290   rw [← eE.sum_inner_mul_inner (B i).val (B j).val]
291   apply Finset.sum_congr rfl
292   intro k _
293   rw [real_inner_comm (B i).val (eE k)]
294 
295 omit [Fintype ι] [DecidableEq ι] [CompleteSpace E] [CompleteSpace Y] in
296 /-- In ambient and residual orthonormal coordinates, the observation matrix is dB. -/
297 theorem protected_operatorMatrix_eq_mul (B : Basis ι ℝ p.ker)
298     (eE : OrthonormalBasis η ℝ E) (eY : OrthonormalBasis κ ℝ Y) :
299     operatorMatrix B eY (ProtectedHilbertRepair.restricted D p) =
300       operatorMatrix eE.toBasis eY D * operatorMatrix B eE p.ker.subtypeL := by
301   ext j i
302   change ⟪eY j, D (B i).val⟫ =
303     ∑ k, ⟪eY j, D (eE k)⟫ * ⟪eE k, (B i).val⟫
304   conv_lhs => rw [← eE.sum_repr' (B i).val]
305   simp only [map_sum, map_smul, inner_sum, real_inner_smul_right]
306   apply Finset.sum_congr rfl
307   intro k _
308   exact mul_comm _ _
309 
310 /-- WR15 on the actual protected kernel: g=B G⁻¹q, with q computed from S=dB. -/
311 theorem protected_gradient_eq_synthesis_inverse_normalResidual
312     (B : Basis ι ℝ p.ker) (eY : OrthonormalBasis κ ℝ Y)
313     (α : ℝ) (c : ι → ℝ) (r : Y) :
314     p.ker.orthogonalProjectionOnto (D.adjoint (D (synthesis B c).val - r)) +
315         α • synthesis B c =
316       synthesis B ((gram B)⁻¹ *ᵥ
317         normalResidual (operatorMatrix B eY (ProtectedHilbertRepair.restricted D p))
318           (gram B) α c (orthonormalCoordinates eY r)) := by
319   rw [← protected_gradient_eq_projected]
320   exact gradient_eq_synthesis_inverse_normalResidual B eY _ α c r
321 
322 /-- WR15 retains the ambient Hilbert norm, including for the zero kernel. -/
323 theorem protected_gradient_norm_sq_eq_normalResidual
324     (B : Basis ι ℝ p.ker) (eY : OrthonormalBasis κ ℝ Y)
325     (α : ℝ) (c : ι → ℝ) (r : Y) :
326     ‖(p.ker.orthogonalProjectionOnto (D.adjoint (D (synthesis B c).val - r))).val +
327         α • (synthesis B c).val‖ ^ 2 =
328       dotProduct
329         (normalResidual (operatorMatrix B eY (ProtectedHilbertRepair.restricted D p))
330           (gram B) α c (orthonormalCoordinates eY r))
331         ((gram B)⁻¹ *ᵥ
332           normalResidual (operatorMatrix B eY (ProtectedHilbertRepair.restricted D p))
333             (gram B) α c (orthonormalCoordinates eY r)) := by
334   have h := gradient_norm_sq_eq_normalResidual B eY
335     (ProtectedHilbertRepair.restricted D p) α c r
336   rw [protected_gradient_eq_projected] at h
337   exact h
338 
339 /-- Complete literal ambient-coordinate certificate.  The matrices d and B use
340 orthonormal ambient/residual coordinates, while B is an arbitrary full kernel
341 basis.  The vector g is the already defined restricted Hilbert normal residual. -/
342 theorem protected_coordinate_certificate
343     (B : Basis ι ℝ p.ker) (eE : OrthonormalBasis η ℝ E)
344     (eY : OrthonormalBasis κ ℝ Y) (α : ℝ) (c : ι → ℝ) (r : Y) :
345     let Bm := operatorMatrix B eE p.ker.subtypeL
346     let dm := operatorMatrix eE.toBasis eY D
347     let S := dm * Bm
348     let G := Bmᵀ * Bm
349     let q := normalResidual S G α c (orthonormalCoordinates eY r)
350     let g := RegularizedHilbertRepair.gradient (ProtectedHilbertRepair.restricted D p)
351       α r (synthesis B c)
352     g = p.ker.orthogonalProjectionOnto (D.adjoint (D (synthesis B c).val - r)) +
353         α • synthesis B c ∧
354       IsUnit G ∧
355       orthonormalCoordinates eE g.val = Bm *ᵥ (G⁻¹ *ᵥ q) ∧
356       ‖g.val‖ ^ 2 = dotProduct q (G⁻¹ *ᵥ q) := by
357   dsimp only
358   rw [← protected_gram_eq_transpose_mul p B eE,
359     ← protected_operatorMatrix_eq_mul D p B eE eY]
360   refine ⟨protected_gradient_eq_projected D p α r _, gram_isUnit B, ?_, ?_⟩
361   · rw [gradient_eq_synthesis_inverse_normalResidual B eY]
362     exact operatorMatrix_synthesis B eE p.ker.subtypeL _
363   · exact gradient_norm_sq_eq_normalResidual B eY _ α c r
364 
365 end ProtectedKernel
366 
367 end KomorebiSpace.CoordinateCovector

Back to source controls