ProposalError.lean

1 import KomorebiLSM.SingularPerturbation
2 import KomorebiLSM.EuclideanMatrixNorms
3 
4 set_option autoImplicit false
5 
6 namespace Komorebi.LSM.ProposalError
7 open scoped Matrix
8 
9 section Abstract
10 variable {E F : Type*} [NormedAddCommGroup E] [InnerProductSpace ℝ E]
11   [NormedAddCommGroup F] [InnerProductSpace ℝ F]
12   [FiniteDimensional ℝ E] [FiniteDimensional ℝ F]
13 
14 /-- The last domain-indexed singular value. In dimension zero this is zero. -/
15 noncomputable def leastSingular (L : E →L[ℝ] F) : ℝ :=
16   L.toLinearMap.singularValues (Module.finrank ℝ E - 1)
17 
18 theorem leastSingular_nonneg (L : E →L[ℝ] F) : 0 ≤ leastSingular L :=
19   L.toLinearMap.singularValues_nonneg _
20 
21 /-- The least singular value controls every domain vector, including a zero-dimensional domain. -/
22 theorem leastSingular_mul_norm_le (L : E →L[ℝ] F) (x : E) :
23     leastSingular L * ‖x‖ ≤ ‖L x‖ := by
24   have ht : Spectral.rightHead L.toLinearMap (Module.finrank ℝ E) le_rfl = ⊤ :=
25     Submodule.eq_top_of_finrank_eq (Spectral.rightHead_finrank _ _ _)
26   exact Spectral.lower_on_rightHead _ _ le_rfl x (by rw [ht]; trivial)
27 
28 /-- Actual singular values obey the operator perturbation bound. -/
29 theorem leastSingular_sub_le (L K : E →L[ℝ] F) (e : ℝ) (he : ‖L-K‖ ≤ e) :
30     leastSingular L - e ≤ leastSingular K := by
31   have h := Spectral.singularValue_le_add_norm L K (Module.finrank ℝ E - 1)
32   change leastSingular L ≤ leastSingular K + ‖L-K‖ at h
33   linarith
34 
35 /-- Residual and operator error give an a posteriori solution error. -/
36 theorem residual_error (L K : E →L[ℝ] F) (x x₀ : E) (b : F) (e : ℝ)
37     (hpert : ‖L-K‖ ≤ e) (hgap : e < leastSingular L) (hexact : K x₀ = b) :
38     ‖x-x₀‖ ≤ (‖L x-b‖ + e*‖x‖) / (leastSingular L-e) := by
39   have hlower := leastSingular_mul_norm_le K (x-x₀)
40   have hs := leastSingular_sub_le L K e hpert
41   have hres : K (x-x₀) = (L x-b) - (L-K) x := by
42     simp only [map_sub, sub_apply, hexact]
43     abel
44   have hu := norm_sub_le (L x-b) ((L-K) x)
45   rw [← hres] at hu
46   have hn : ‖(L-K) x‖ ≤ e*‖x‖ :=
47     (ContinuousLinearMap.le_opNorm _ _).trans
48       (mul_le_mul_of_nonneg_right hpert (norm_nonneg _))
49   apply (le_div_iff₀ (sub_pos.mpr hgap)).mpr
50   have hh := mul_le_mul_of_nonneg_right hs (norm_nonneg (x-x₀))
51   nlinarith
52 
53 theorem reference_injective (L K : E →L[ℝ] F) (e : ℝ)
54     (hpert : ‖L-K‖ ≤ e) (hgap : e < leastSingular L) : Function.Injective K := by
55   intro x y hxy
56   have hl := leastSingular_mul_norm_le K (x-y)
57   have hs := leastSingular_sub_le L K e hpert
58   have hp : 0 < leastSingular K := by linarith
59   have hz : K (x-y) = 0 := by simp [hxy]
60   rw [hz, norm_zero] at hl
61   have : ‖x-y‖ = 0 := by nlinarith [norm_nonneg (x-y)]
62   exact sub_eq_zero.mp (norm_eq_zero.mp this)
63 end Abstract
64 
65 section Matrices
66 variable {D : Type*} [Fintype D] [DecidableEq D]
67 
68 /-- Symmetric matrices with their entrywise Euclidean (Frobenius) inner product. -/
69 def symmetricSpace (D : Type*) [Fintype D] :
70     Submodule ℝ (EuclideanSpace ℝ (D × D)) where
71   carrier := {H | ∀ i j, H (i,j) = H (j,i)}
72   zero_mem' := by simp
73   add_mem' := by intro H K hH hK i j; simp [hH i j, hK i j]
74   smul_mem' := by intro a H hH i j; simp [hH i j]
75 
76 abbrev Symmetric (D : Type*) [Fintype D] := symmetricSpace D
77 
78 /-- The literal matrix represented by a symmetric Frobenius vector. -/
79 def matrix (H : Symmetric D) : Matrix D D ℝ := fun i j => H.val (i,j)
80 
81 omit [DecidableEq D] in
82 @[simp] theorem matrix_add (H K : Symmetric D) : matrix (H+K) = matrix H+matrix K := rfl
83 omit [DecidableEq D] in
84 @[simp] theorem matrix_sub (H K : Symmetric D) : matrix (H-K) = matrix H-matrix K := rfl
85 omit [DecidableEq D] in
86 @[simp] theorem matrix_smul (a : ℝ) (H : Symmetric D) : matrix (a•H) = a•matrix H := rfl
87 
88 omit [DecidableEq D] in
89 theorem matrix_symmetric (H : Symmetric D) : (matrix H)ᵀ = matrix H := by
90   ext i j
91   exact H.property j i
92 
93 /-- The actual invariant-form/trace linear map, with the Euclidean product norm. -/
94 noncomputable def constraintLinear (A : Matrix D D ℝ) :
95     Symmetric D →ₗ[ℝ] EuclideanSpace ℝ ((D × D) ⊕ Unit) where
96   toFun H := WithLp.toLp 2 (Sum.elim
97     (fun ij => (Aᵀ * matrix H * A - matrix H) ij.1 ij.2)
98     (fun _ => Matrix.trace (matrix H)))
99   map_add' H K := by
100     ext x
101     cases x with
102     | inl ij => simp [Matrix.mul_add, Matrix.add_mul]; ring
103     | inr u => simp [Matrix.trace_add]
104   map_smul' a H := by
105     ext x
106     cases x with
107     | inl ij => simp; ring
108     | inr u => simp [Matrix.trace_smul]
109 
110 noncomputable def constraint (A : Matrix D D ℝ) :
111     Symmetric D →L[ℝ] EuclideanSpace ℝ ((D × D) ⊕ Unit) :=
112   (constraintLinear A).toContinuousLinearMap
113 
114 @[simp] theorem constraint_matrix (A : Matrix D D ℝ) (H : Symmetric D) (i j : D) :
115     constraint A H (Sum.inl (i,j)) = (Aᵀ * matrix H * A - matrix H) i j := rfl
116 
117 @[simp] theorem constraint_trace (A : Matrix D D ℝ) (H : Symmetric D) :
118     constraint A H (Sum.inr ()) = Matrix.trace (matrix H) := rfl
119 
120 /-- The trace-normalized invariant-form target `(0,1)`. -/
121 noncomputable def target : EuclideanSpace ℝ ((D × D) ⊕ Unit) :=
122   WithLp.toLp 2 (Sum.elim (fun _ => 0) (fun _ => 1))
123 
124 noncomputable def residual (A : Matrix D D ℝ) (H : Symmetric D) :=
125   constraint A H - (target : EuclideanSpace ℝ ((D × D) ⊕ Unit))
126 
127 open EuclideanMatrixNorms
128 
129 omit [DecidableEq D] in
130 theorem frobenius_matrix (H : Symmetric D) : frobenius (matrix H) = ‖H‖ := by
131   rfl
132 
133 omit [DecidableEq D] in
134 /-- A zero trace coordinate does not alter the Frobenius norm. -/
135 theorem zero_trace_norm (B : Matrix D D ℝ) :
136     ‖(WithLp.toLp 2 (Sum.elim (fun ij : D × D => B ij.1 ij.2)
137       (fun _ : Unit => (0 : ℝ))) : EuclideanSpace ℝ ((D × D) ⊕ Unit))‖ = frobenius B := by
138   have hs : ‖(WithLp.toLp 2 (Sum.elim (fun ij : D × D => B ij.1 ij.2)
139       (fun _ : Unit => (0 : ℝ))) : EuclideanSpace ℝ ((D × D) ⊕ Unit))‖ ^ 2 =
140       frobenius B ^ 2 := by
141     rw [EuclideanSpace.real_norm_sq_eq, frobenius_sq]
142     simp [Fintype.sum_sum_type, Fintype.sum_prod_type]
143   nlinarith [norm_nonneg (WithLp.toLp 2 (Sum.elim (fun ij : D × D => B ij.1 ij.2)
144       (fun _ : Unit => (0 : ℝ))) : EuclideanSpace ℝ ((D × D) ⊕ Unit)), frobenius_nonneg B]
145 
146 theorem constraint_difference_norm (A B : Matrix D D ℝ) (H : Symmetric D) :
147     ‖(constraint A-constraint B) H‖ =
148       frobenius (Aᵀ * matrix H * A - Bᵀ * matrix H * B) := by
149   have he : (constraint A-constraint B) H =
150       WithLp.toLp 2 (Sum.elim
151         (fun ij : D × D => (Aᵀ * matrix H * A - Bᵀ * matrix H * B) ij.1 ij.2)
152         (fun _ : Unit => (0 : ℝ))) := by
153     ext x
154     cases x with
155     | inl ij =>
156       rcases ij with ⟨i,j⟩
157       simp
158     | inr u => cases u; simp
159   rw [he, zero_trace_norm]
160 
161 /-- The trace cancels, so the actual induced perturbation is bounded by Euclidean matrix errors. -/
162 theorem constraint_perturbation (A B : Matrix D D ℝ) :
163     ‖constraint A-constraint B‖ ≤
164       (spectralNorm A + spectralNorm B) * spectralNorm (A-B) := by
165   apply ContinuousLinearMap.opNorm_le_bound _
166     (mul_nonneg (add_nonneg (spectralNorm_nonneg A) (spectralNorm_nonneg B))
167       (spectralNorm_nonneg (A-B)))
168   intro H
169   rw [constraint_difference_norm]
170   simpa only [frobenius_matrix, mul_assoc] using frobenius_congruence_sub_le A B (matrix H)
171 
172 /-- The proposal error uses the literal invariant-form map, not an assumed residual identity. -/
173 theorem invariant_form_error (Ahat Astar : Matrix D D ℝ) (Hhat Hstar : Symmetric D)
174     (eM : ℝ) (hpert : ‖constraint Ahat - constraint Astar‖ ≤ eM)
175     (hgap : eM < leastSingular (constraint Ahat))
176     (hexact : constraint Astar Hstar = target) :
177     ‖Hhat-Hstar‖ ≤ (‖residual Ahat Hhat‖ + eM*‖Hhat‖) /
178       (leastSingular (constraint Ahat)-eM) :=
179   residual_error _ _ _ _ _ _ hpert hgap hexact
180 
181 /-- Literal Frobenius form of the a posteriori invariant-form error bound. -/
182 theorem invariant_form_frobenius_error (Ahat Astar : Matrix D D ℝ)
183     (Hhat Hstar : Symmetric D) (eM : ℝ)
184     (hpert : ‖constraint Ahat-constraint Astar‖ ≤ eM)
185     (hgap : eM < leastSingular (constraint Ahat))
186     (hexact : constraint Astar Hstar = target) :
187     frobenius (matrix Hhat-matrix Hstar) ≤
188       (‖residual Ahat Hhat‖+eM*frobenius (matrix Hhat)) /
189       (leastSingular (constraint Ahat)-eM) := by
190   simpa only [← matrix_sub, frobenius_matrix] using
191     invariant_form_error Ahat Astar Hhat Hstar eM hpert hgap hexact
192 
193 /-- The exact target equation is precisely invariance and trace one. -/
194 theorem constraint_eq_target_iff (A : Matrix D D ℝ) (H : Symmetric D) :
195     constraint A H = target ↔ Aᵀ*matrix H*A = matrix H ∧ Matrix.trace (matrix H) = 1 := by
196   constructor
197   · intro he
198     constructor
199     · have hz : Aᵀ*matrix H*A-matrix H = 0 := by
200         ext i j
201         have h := congrArg (fun x : EuclideanSpace ℝ ((D × D) ⊕ Unit) => x (Sum.inl (i,j))) he
202         simpa only [constraint_matrix, target, WithLp.toLp_ofLp, Sum.elim_inl,
203           Matrix.sub_apply, Matrix.zero_apply] using h
204       exact sub_eq_zero.mp hz
205     · have h := congrArg (fun x : EuclideanSpace ℝ ((D × D) ⊕ Unit) => x (Sum.inr ())) he
206       exact h
207   · rintro ⟨hinv,htr⟩
208     ext x
209     cases x with
210     | inl ij =>
211       rcases ij with ⟨i,j⟩
212       simp [hinv, target]
213     | inr u => cases u; simpa [target] using htr
214 
215 /-- A matrix-budget version derives the needed map budget from actual matrices. -/
216 theorem invariant_form_matrix_budget (Ahat Astar : Matrix D D ℝ)
217     (Hhat Hstar : Symmetric D) (eM : ℝ)
218     (hbudget : (spectralNorm Ahat+spectralNorm Astar)*spectralNorm (Ahat-Astar) ≤ eM)
219     (hgap : eM < leastSingular (constraint Ahat))
220     (hinv : Astarᵀ*matrix Hstar*Astar = matrix Hstar) (htrace : Matrix.trace (matrix Hstar)=1) :
221     frobenius (matrix Hhat-matrix Hstar) ≤
222       (‖residual Ahat Hhat‖+eM*frobenius (matrix Hhat)) /
223       (leastSingular (constraint Ahat)-eM) :=
224   invariant_form_frobenius_error _ _ _ _ _
225     ((constraint_perturbation _ _).trans hbudget) hgap
226     ((constraint_eq_target_iff _ _).mpr ⟨hinv,htrace⟩)
227 
228 /-- This is exactly the Euclidean product of Frobenius residual and trace. -/
229 theorem constraint_norm_sq (A : Matrix D D ℝ) (H : Symmetric D) :
230     ‖constraint A H‖ ^ 2 = frobenius (Aᵀ*matrix H*A-matrix H)^2 +
231       (Matrix.trace (matrix H))^2 := by
232   rw [EuclideanSpace.real_norm_sq_eq, frobenius_sq]
233   simp [Fintype.sum_sum_type, Fintype.sum_prod_type, constraint, constraintLinear]
234 
235 /-- The center is the actual nonsingular-inverse formula, in Euclidean coordinates. -/
236 noncomputable def fixedCenter (A : Matrix D D ℝ) (b : D → ℝ) : EuclideanSpace ℝ D :=
237   WithLp.toLp 2 ((1-A)⁻¹ *ᵥ b)
238 
239 theorem fixedCenter_equation (A : Matrix D D ℝ) (b : D → ℝ)
240     (hunit : IsUnit (1-A).det) :
241     Matrix.toEuclideanCLM (𝕜 := ℝ) (1-A) (fixedCenter A b) = WithLp.toLp 2 b := by
242   change WithLp.toLp 2 ((1-A) *ᵥ ((1-A)⁻¹ *ᵥ b)) = _
243   rw [Matrix.mulVec_mulVec, Matrix.mul_nonsing_inv _ hunit, Matrix.one_mulVec]
244 
245 theorem center_operator_perturbation (A B : Matrix D D ℝ) :
246     ‖Matrix.toEuclideanCLM (𝕜 := ℝ) (1-A) - Matrix.toEuclideanCLM (𝕜 := ℝ) (1-B)‖ =
247       spectralNorm (A-B) := by
248   rw [← map_sub]
249   have he : (1-A)-(1-B) = -(A-B) := by abel
250   rw [he, map_neg, norm_neg]
251   rfl
252 
253 /-- The fixed-center estimate uses actual inverse centers and actual spectral perturbation. -/
254 theorem fixed_center_error (Ahat Astar : Matrix D D ℝ) (bhat bstar : D → ℝ)
255     (hhat : IsUnit (1-Ahat).det) (hstar : IsUnit (1-Astar).det)
256     (eA eb : ℝ) (hA : spectralNorm (Ahat-Astar) ≤ eA)
257     (hb : ‖(WithLp.toLp 2 (bhat-bstar) : EuclideanSpace ℝ D)‖ ≤ eb)
258     (hgap : eA < leastSingular (Matrix.toEuclideanCLM (𝕜 := ℝ) (1-Ahat))) :
259     ‖fixedCenter Ahat bhat-fixedCenter Astar bstar‖ ≤
260       (eb+eA*‖fixedCenter Ahat bhat‖) /
261         (leastSingular (Matrix.toEuclideanCLM (𝕜 := ℝ) (1-Ahat))-eA) := by
262   have hpert : ‖Matrix.toEuclideanCLM (𝕜 := ℝ) (1-Ahat) -
263       Matrix.toEuclideanCLM (𝕜 := ℝ) (1-Astar)‖ ≤ eA := by
264     rw [center_operator_perturbation]
265     exact hA
266   have he := residual_error (E := EuclideanSpace ℝ D) (F := EuclideanSpace ℝ D)
267     (Matrix.toEuclideanCLM (𝕜 := ℝ) (1-Ahat))
268     (Matrix.toEuclideanCLM (𝕜 := ℝ) (1-Astar))
269     (fixedCenter Ahat bhat) (fixedCenter Astar bstar) (WithLp.toLp 2 bstar)
270     eA hpert hgap (fixedCenter_equation Astar bstar hstar)
271   rw [fixedCenter_equation Ahat bhat hhat] at he
272   have hbres : ‖(WithLp.toLp 2 bhat : EuclideanSpace ℝ D)-WithLp.toLp 2 bstar‖ ≤ eb := hb
273   exact he.trans ((div_le_div_iff_of_pos_right (sub_pos.mpr hgap)).mpr (by linarith))
274 
275 end Matrices
276 
277 end Komorebi.LSM.ProposalError

Back to source controls