BlockJacobi.lean

1 module
2 public import Mathlib.Analysis.Matrix.Order
3 public import Mathlib.Analysis.Real.Sqrt
4 public import Mathlib.Tactic
5 
6 /-!
7 Author: Komorebi P. Paper VIII items 7.4 and 7.6.
8 The normalization uses the actual positive matrix square root and actual inverse.
9 Condition number is the product of Euclidean operator norms. The PCG result
10 compares the classical convergence factors; it does not reprove the CG algorithm.
11 -/
12 public section
13 set_option autoImplicit false
14 namespace KomorebiPaperVIII
15 open scoped BigOperators Matrix Matrix.Norms.L2Operator ComplexOrder MatrixOrder
16 
17 section Ring
18 variable {R : Type*} [NormedRing R] [NormOneClass R]
19 
20 theorem inverse_near_one_bound (E B : R) (hB : (1+E)*B=1)
21     (η : ℝ) (hE : ‖E‖ ≤ η) (hη : η < 1) : ‖B‖ ≤ 1/(1-η) := by
22   have heq : B = 1-E*B := by
23     apply eq_sub_iff_add_eq.mpr
24     simpa only [add_mul,one_mul] using hB
25   have hb : ‖B‖ ≤ 1+η*‖B‖ := by
26     calc
27       ‖B‖ = ‖1-E*B‖ := congrArg norm heq
28       _ ≤ ‖(1:R)‖+‖E*B‖ := norm_sub_le _ _
29       _ ≤ 1+‖E‖*‖B‖ := by
30         simpa only [norm_one] using
31           (add_le_add (le_refl ‖(1:R)‖) (norm_mul_le E B))
32       _ ≤ 1+η*‖B‖ := by
33         gcongr
34   apply (le_div_iff₀ (by linarith : 0 < 1-η)).mpr
35   nlinarith
36 
37 theorem near_one_condition_bound (E B : R) (hB : (1+E)*B=1)
38     (η : ℝ) (hE : ‖E‖ ≤ η) (hη : η < 1) :
39     ‖1+E‖*‖B‖ ≤ (1+η)/(1-η) := by
40   have hη0 : 0 ≤ η := (norm_nonneg E).trans hE
41   have hn : ‖1+E‖ ≤ 1+η := by
42     have hadd := norm_add_le (1:R) E
43     rw [norm_one] at hadd
44     linarith
45   have hb := inverse_near_one_bound E B hB η hE hη
46   calc
47     ‖1+E‖*‖B‖ ≤ (1+η)*(1/(1-η)) :=
48       mul_le_mul hn hb (norm_nonneg B) (by linarith)
49     _ = (1+η)/(1-η) := by ring
50 end Ring
51 
52 section Matrix
53 variable {X : Type*} [Fintype X] [DecidableEq X]
54 
55 noncomputable def invSqrt (D : Matrix X X ℂ) : Matrix X X ℂ := (CFC.sqrt D)⁻¹
56 
57 noncomputable def matrixCondition (A : Matrix X X ℂ) : ℝ := ‖A‖*‖A⁻¹‖
58 
59 theorem invSqrt_isUnit (D : Matrix X X ℂ) (hD : D.PosDef) : IsUnit (invSqrt D) := by
60   unfold invSqrt
61   apply Matrix.isUnit_nonsing_inv_iff.mpr
62   exact (CFC.isUnit_sqrt_iff D hD.posSemidef.nonneg).mpr hD.isUnit
63 
64 /-- The actual positive inverse square root whitens D. -/
65 theorem invSqrt_normalizes (D : Matrix X X ℂ) (hD : D.PosDef) :
66     invSqrt D*D*invSqrt D = 1 := by
67   have hroot : IsUnit (CFC.sqrt D) :=
68     (CFC.isUnit_sqrt_iff D hD.posSemidef.nonneg).mpr hD.isUnit
69   have hd := (Matrix.isUnit_iff_isUnit_det _).mp hroot
70   have hs := CFC.sqrt_mul_sqrt_self D hD.posSemidef.nonneg
71   unfold invSqrt
72   calc
73     (CFC.sqrt D)⁻¹*D*(CFC.sqrt D)⁻¹ =
74         (CFC.sqrt D)⁻¹*(CFC.sqrt D*CFC.sqrt D)*(CFC.sqrt D)⁻¹ :=
75       congrArg (fun M => (CFC.sqrt D)⁻¹*M*(CFC.sqrt D)⁻¹) hs.symm
76     _ = ((CFC.sqrt D)⁻¹*CFC.sqrt D)*(CFC.sqrt D*(CFC.sqrt D)⁻¹) := by
77       simp only [Matrix.mul_assoc]
78     _ = 1 := by rw [Matrix.nonsing_inv_mul _ hd,Matrix.mul_nonsing_inv _ hd,one_mul]
79 
80 theorem normalized_decomposition (D K : Matrix X X ℂ) (hD : D.PosDef) :
81     invSqrt D*(D+K)*invSqrt D = 1+invSqrt D*K*invSqrt D := by
82   rw [Matrix.mul_add,Matrix.add_mul,invSqrt_normalizes D hD]
83 
84 variable [Nonempty X]
85 
86 /-- Actual condition bound. Hermiticity of the perturbation is not needed by this stronger estimate. -/
87 theorem block_jacobi_condition (D K : Matrix X X ℂ)
88     (hD : D.PosDef) (hA : (D+K).PosDef) (η : ℝ)
89     (hE : ‖invSqrt D*K*invSqrt D‖ ≤ η) (hη : η < 1) :
90     matrixCondition (invSqrt D*(D+K)*invSqrt D) ≤ (1+η)/(1-η) := by
91   have hunit := ((invSqrt_isUnit D hD).mul hA.isUnit).mul (invSqrt_isUnit D hD)
92   have hdet := (Matrix.isUnit_iff_isUnit_det _).mp hunit
93   have hinv := Matrix.mul_nonsing_inv (invSqrt D*(D+K)*invSqrt D) hdet
94   rw [normalized_decomposition D K hD] at hinv ⊢
95   exact near_one_condition_bound _ _ hinv η hE hη
96 
97 theorem matrix_condition_ge_one (A : Matrix X X ℂ) (hA : IsUnit A) :
98     1 ≤ matrixCondition A := by
99   have hn := norm_mul_le A A⁻¹
100   rw [Matrix.mul_nonsing_inv _ ((Matrix.isUnit_iff_isUnit_det _).mp hA),norm_one] at hn
101   exact hn
102 end Matrix
103 
104 noncomputable def cgFactor (κ : ℝ) : ℝ := (Real.sqrt κ-1)/(Real.sqrt κ+1)
105 
106 theorem cgFactor_mono {κ κmax : ℝ} (h : κ ≤ κmax) : cgFactor κ ≤ cgFactor κmax := by
107   have hs := Real.sqrt_le_sqrt h
108   have hκ := Real.sqrt_nonneg κ
109   have hκmax := Real.sqrt_nonneg κmax
110   unfold cgFactor
111   apply (div_le_div_iff₀ (by linarith) (by linarith)).mpr
112   nlinarith
113 
114 theorem cgFactor_nonneg {κ : ℝ} (h : 1 ≤ κ) : 0 ≤ cgFactor κ := by
115   have hs := Real.one_le_sqrt.mpr h
116   unfold cgFactor
117   positivity
118 
119 theorem cgFactor_lt_one (κ : ℝ) : cgFactor κ < 1 := by
120   have hs := Real.sqrt_nonneg κ
121   unfold cgFactor
122   apply (div_lt_one (by linarith)).mpr
123   linarith
124 
125 /-- Item 7.6 compares the actual normalized matrix condition-number factor. -/
126 theorem certified_packet_factor {X : Type*} [Fintype X] [DecidableEq X] [Nonempty X]
127     (D K : Matrix X X ℂ) (hD : D.PosDef) (hA : (D+K).PosDef)
128     (η : ℝ) (hE : ‖invSqrt D*K*invSqrt D‖ ≤ η) (hη : η < 1) :
129     cgFactor (matrixCondition (invSqrt D*(D+K)*invSqrt D)) ≤
130       cgFactor ((1+η)/(1-η)) :=
131   cgFactor_mono (block_jacobi_condition D K hD hA η hE hη)
132 
133 /-- Explicit propagation of an already established classical energy-error estimate. -/
134 theorem certified_energy_error {κ κmax e0 em : ℝ} (m : ℕ)
135     (hκ : 1 ≤ κ) (hbound : κ ≤ κmax) (he0 : 0 ≤ e0)
136     (hCG : em ≤ 2*(cgFactor κ)^m*e0) :
137     em ≤ 2*(cgFactor κmax)^m*e0 := by
138   have hp : (cgFactor κ)^m ≤ (cgFactor κmax)^m :=
139     pow_le_pow_left₀ (cgFactor_nonneg hκ) (cgFactor_mono hbound) m
140   exact hCG.trans (mul_le_mul_of_nonneg_right
141     (mul_le_mul_of_nonneg_left hp (by norm_num)) he0)
142 
143 end KomorebiPaperVIII
144 
145 -- FIFTH BLOCK AXIOM AUDIT
146 #print axioms KomorebiPaperVIII.inverse_near_one_bound
147 #print axioms KomorebiPaperVIII.near_one_condition_bound
148 #print axioms KomorebiPaperVIII.invSqrt_isUnit
149 #print axioms KomorebiPaperVIII.invSqrt_normalizes
150 #print axioms KomorebiPaperVIII.normalized_decomposition
151 #print axioms KomorebiPaperVIII.block_jacobi_condition
152 #print axioms KomorebiPaperVIII.matrix_condition_ge_one
153 #print axioms KomorebiPaperVIII.cgFactor_mono
154 #print axioms KomorebiPaperVIII.cgFactor_nonneg
155 #print axioms KomorebiPaperVIII.cgFactor_lt_one
156 #print axioms KomorebiPaperVIII.certified_packet_factor
157 #print axioms KomorebiPaperVIII.certified_energy_error
158 #print axioms KomorebiPaperVIII.invSqrt
159 #print axioms KomorebiPaperVIII.matrixCondition
160 #print axioms KomorebiPaperVIII.cgFactor

Back to source controls