LowRankAlgebra.lean

1 module
2 public import PacketIR
3 import all PacketIR
4 import all QuantumPacketCore
5 public import Mathlib.LinearAlgebra.Matrix.NonsingularInverse
6 public import Mathlib.LinearAlgebra.Matrix.Rank
7 public import Mathlib.Data.Matrix.ColumnRowPartitioned
8 public import Mathlib.Analysis.CStarAlgebra.Matrix
9 public import Mathlib.Tactic.Abel
10 
11 /-!
12 Author: Komorebi P.
13 Paper VIII items 5.1 and 5.3: actual rectangular factors and actual matrix inverse.
14 Woodbury specializes Mathlib's classical identity; no novelty is claimed for it.
15 The inverse is not used as an invertibility assumption: update invertibility is proved.
16 -/
17 
18 public section
19 set_option autoImplicit false
20 set_option backward.isDefEq.respectTransparency false
21 set_option backward.isDefEq.respectTransparency.types false
22 namespace KomorebiPaperVIII
23 open scoped BigOperators
24 open KomorebiPaperVII
25 
26 section Diagonal
27 variable {R I : Type*} [Ring R] [Fintype I] [DecidableEq I]
28 
29 /-- Diagonality is an exact vanishing condition, not a numerical threshold. -/
30 def PacketDiagonal (P : PacketFamily R I) (D : R) : Prop :=
31   ∀ b a, b ≠ a → P.block D b a = 0
32 
33 theorem packet_diagonal_product (P : PacketFamily R I) (D E : R)
34     (hD : PacketDiagonal P D) (hE : PacketDiagonal P E) : PacketDiagonal P (D * E) := by
35   intro b a hba
36   rw [P.compose_blocks]
37   apply Finset.sum_eq_zero
38   intro c _
39   by_cases hb : b = c
40   · subst c
41     rw [hE b a hba, mul_zero]
42   · rw [hD b c hb, zero_mul]
43 end Diagonal
44 
45 section Product
46 variable {K X RA RB : Type*} [Field K] [StarRing K]
47 variable [Fintype X] [Fintype RA] [Fintype RB]
48 
49 /-- The first factor has r_B+r_A columns, in the manuscript's displayed order. -/
50 def bdlrProductLeft (DA : Matrix X X K) (UA VA : Matrix X RA K) (UB : Matrix X RB K) :
51     Matrix X (RB ⊕ RA) K := Matrix.fromCols (DA * UB + UA * (VA.conjTranspose * UB)) UA
52 
53 def bdlrProductRight (DB : Matrix X X K) (VA : Matrix X RA K) (VB : Matrix X RB K) :
54     Matrix X (RB ⊕ RA) K := Matrix.fromCols VB (DB.conjTranspose * VA)
55 
56 /-- Theorem 5.1: multiplication without materializing an unstructured correction. -/
57 theorem bdlr_product (DA DB : Matrix X X K) (UA VA : Matrix X RA K)
58     (UB VB : Matrix X RB K) :
59     (DA + UA * VA.conjTranspose) * (DB + UB * VB.conjTranspose) =
60       DA * DB + bdlrProductLeft DA UA VA UB * (bdlrProductRight DB VA VB).conjTranspose := by
61   simp only [bdlrProductLeft, bdlrProductRight,
62     Matrix.conjTranspose_fromCols_eq_fromRows_conjTranspose,
63     Matrix.fromCols_mul_fromRows, Matrix.conjTranspose_mul, Matrix.conjTranspose_conjTranspose,
64     Matrix.add_mul, Matrix.mul_add, Matrix.mul_assoc]
65   abel
66 
67 /-- The rank statement is actual matrix rank, not only a count of allocated columns. -/
68 theorem bdlr_product_rank_bound (DA DB : Matrix X X K) (UA VA : Matrix X RA K)
69     (UB VB : Matrix X RB K) :
70     ((DA + UA * VA.conjTranspose) * (DB + UB * VB.conjTranspose) - DA * DB).rank ≤
71       Fintype.card RA + Fintype.card RB := by
72   rw [bdlr_product, add_sub_cancel_left]
73   calc
74     _ ≤ (bdlrProductLeft DA UA VA UB).rank := Matrix.rank_mul_le_left _ _
75     _ ≤ Fintype.card (RB ⊕ RA) := Matrix.rank_le_card_width _
76     _ = _ := by simp [Nat.add_comm]
77 
78 /-- An arbitrary batch of right-hand sides uses the two thin factors directly. -/
79 theorem bdlr_apply {C : Type*} (D : Matrix X X K) (U V : Matrix X RA K) (Z : Matrix X C K) :
80     (D + U * V.conjTranspose) * Z = D * Z + U * (V.conjTranspose * Z) := by
81   simp only [Matrix.add_mul, Matrix.mul_assoc]
82 end Product
83 
84 section Inverse
85 variable {K X Rk : Type*} [Field K] [StarRing K]
86 variable [Fintype X] [DecidableEq X] [Fintype Rk] [DecidableEq Rk]
87 
88 /-- Theorem 5.3: this is the ordinary nonsingular matrix inverse, with both needed hypotheses. -/
89 theorem bdlr_woodbury (D : Matrix X X K) (U V : Matrix X Rk K)
90     (hD : IsUnit D) (hS : IsUnit (1 + V.conjTranspose * D⁻¹ * U)) :
91     (D + U * V.conjTranspose)⁻¹ =
92       D⁻¹ - D⁻¹ * U * (1 + V.conjTranspose * D⁻¹ * U)⁻¹ * V.conjTranspose * D⁻¹ := by
93   simpa using Matrix.add_mul_mul_inv_eq_sub D U (1 : Matrix Rk Rk K)
94     V.conjTranspose hD isUnit_one (by simpa using hS)
95 
96 /-- The update really is invertible. The conclusion is not smuggled in as an assumption. -/
97 theorem bdlr_update_isUnit (D : Matrix X X K) (U V : Matrix X Rk K)
98     (hD : IsUnit D) (hS : IsUnit (1 + V.conjTranspose * D⁻¹ * U)) :
99     IsUnit (D + U * V.conjTranspose) := by
100   classical
101   obtain ⟨_⟩ := hD.nonempty_invertible
102   let _ := invertibleOne (α := Matrix Rk Rk K)
103   have hs : IsUnit (⅟(1 : Matrix Rk Rk K) + V.conjTranspose * ⅟D * U) := by
104     simpa only [invOf_one, Matrix.invOf_eq_nonsing_inv] using hS
105   obtain ⟨_⟩ := hs.nonempty_invertible
106   let _ := Matrix.invertibleAddMulMul D U (1 : Matrix Rk Rk K) V.conjTranspose
107   simpa using (isUnit_of_invertible (D + U * (1 : Matrix Rk Rk K) * V.conjTranspose))
108 
109 /-- The proposed inverse solves the actual updated system, for every batch Z. -/
110 theorem bdlr_woodbury_solve {C : Type*} (D : Matrix X X K) (U V : Matrix X Rk K)
111     (hD : IsUnit D) (hS : IsUnit (1 + V.conjTranspose * D⁻¹ * U)) (Z : Matrix X C K) :
112     (D + U * V.conjTranspose) *
113       ((D⁻¹ - D⁻¹ * U * (1 + V.conjTranspose * D⁻¹ * U)⁻¹ * V.conjTranspose * D⁻¹) * Z) = Z := by
114   rw [← bdlr_woodbury D U V hD hS, ← Matrix.mul_assoc]
115   have hu := bdlr_update_isUnit D U V hD hS
116   rw [Matrix.mul_nonsing_inv _ ((Matrix.isUnit_iff_isUnit_det _).mp hu), Matrix.one_mul]
117 end Inverse
118 
119 
120 section InverseNorm
121 open scoped Matrix.Norms.L2Operator
122 variable {K X Rk : Type*} [RCLike K]
123 variable [Fintype X] [DecidableEq X] [Fintype Rk] [DecidableEq Rk]
124 
125 /-- Corollary 5.4 in the Hilbert-space operator norm, including rectangular thin factors. -/
126 theorem bdlr_inverse_norm_bound (D : Matrix X X K) (U V : Matrix X Rk K)
127     (hD : IsUnit D) (hS : IsUnit (1 + V.conjTranspose * D⁻¹ * U)) :
128     ‖(D + U * V.conjTranspose)⁻¹‖ ≤ ‖D⁻¹‖ + ‖D⁻¹ * U‖ *
129       ‖(1 + V.conjTranspose * D⁻¹ * U)⁻¹‖ * ‖V.conjTranspose * D⁻¹‖ := by
130   rw [bdlr_woodbury D U V hD hS]
131   have hfactor : D⁻¹ * U * (1 + V.conjTranspose * D⁻¹ * U)⁻¹ * V.conjTranspose * D⁻¹ =
132       (D⁻¹ * U) * (1 + V.conjTranspose * D⁻¹ * U)⁻¹ * (V.conjTranspose * D⁻¹) := by
133     simp only [Matrix.mul_assoc]
134   rw [hfactor]
135   apply (norm_sub_le _ _).trans
136   apply add_le_add le_rfl
137   exact (Matrix.l2_opNorm_mul _ _).trans
138     (mul_le_mul_of_nonneg_right (Matrix.l2_opNorm_mul _ _) (norm_nonneg _))
139 end InverseNorm
140 
141 end KomorebiPaperVIII
142 
143 #print axioms KomorebiPaperVIII.packet_diagonal_product
144 #print axioms KomorebiPaperVIII.bdlr_product
145 #print axioms KomorebiPaperVIII.bdlr_product_rank_bound
146 #print axioms KomorebiPaperVIII.bdlr_apply
147 #print axioms KomorebiPaperVIII.bdlr_woodbury
148 #print axioms KomorebiPaperVIII.bdlr_update_isUnit
149 #print axioms KomorebiPaperVIII.bdlr_woodbury_solve
150 
151 #print axioms KomorebiPaperVIII.bdlr_inverse_norm_bound

Back to source controls