[Merged by Bors] - feat(LinearAlgebra/Matrix): Woodbury Identity

Mathlib/Data/Matrix/Invertible.lean

@@ -1,9 +1,10 @@
 /-
 Copyright (c) 2023 Eric Wieser. All rights reserved.
 Released under Apache 2.0 license as described in the file LICENSE.
-Authors: Eric Wieser
+Authors: Eric Wieser, Ahmad Alkhalawi
 -/
 import Mathlib.Data.Matrix.Basic
+import Mathlib.Tactic.Abel
 
 /-! # Extra lemmas about invertible matrices
 
@@ -47,10 +48,12 @@ protected theorem invOf_mul_cancel_right (A : Matrix m n α) (B : Matrix n n α)
 protected theorem mul_invOf_cancel_right (A : Matrix m n α) (B : Matrix n n α) [Invertible B] :
     A * B * ⅟ B = A := by rw [Matrix.mul_assoc, mul_invOf_self, Matrix.mul_one]
 
-@[deprecated (since := "2024-09-07")] alias invOf_mul_self_assoc := invOf_mul_cancel_left
-@[deprecated (since := "2024-09-07")] alias mul_invOf_self_assoc := mul_invOf_cancel_left
-@[deprecated (since := "2024-09-07")] alias mul_invOf_mul_self_cancel := invOf_mul_cancel_right
-@[deprecated (since := "2024-09-07")] alias mul_mul_invOf_self_cancel := mul_invOf_cancel_right
+@[deprecated (since := "2024-09-07")] alias invOf_mul_self_assoc := Matrix.invOf_mul_cancel_left
+@[deprecated (since := "2024-09-07")] alias mul_invOf_self_assoc := Matrix.mul_invOf_cancel_left
+@[deprecated (since := "2024-09-07")]
+alias mul_invOf_mul_self_cancel := Matrix.invOf_mul_cancel_right
+@[deprecated (since := "2024-09-07")]
+alias mul_mul_invOf_self_cancel := Matrix.mul_invOf_cancel_right
 
 section ConjTranspose
 variable [StarRing α] (A : Matrix n n α)
@@ -106,4 +109,68 @@ def transposeInvertibleEquivInvertible : Invertible Aᵀ ≃ Invertible A where
 
 end CommSemiring
 
+section Ring
+
+section Woodbury
+
+variable [Fintype m] [DecidableEq m] [Ring α]
+    (A : Matrix n n α) (U : Matrix n m α) (C : Matrix m m α) (V : Matrix m n α)
+    [Invertible A] [Invertible C] [Invertible (⅟C + V * ⅟A * U)]
+
+-- No spaces around multiplication signs for better clarity
+lemma add_mul_mul_invOf_mul_eq_one :
+    (A + U*C*V)*(⅟A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A) = 1 := by
+  calc
+    (A + U*C*V)*(⅟A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A)
+    _ = A*⅟A - A*⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A + U*C*V*⅟A - U*C*V*⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A := by
+      simp_rw [add_sub_assoc, add_mul, mul_sub, Matrix.mul_assoc]
+    _ = (1 + U*C*V*⅟A) - (U*⅟(⅟C + V*⅟A*U)*V*⅟A + U*C*V*⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A) := by
+      rw [mul_invOf_self, Matrix.one_mul]
+      abel
+    _ = 1 + U*C*V*⅟A - (U + U*C*V*⅟A*U)*⅟(⅟C + V*⅟A*U)*V*⅟A := by
+      rw [sub_right_inj, Matrix.add_mul, Matrix.add_mul, Matrix.add_mul]
+    _ = 1 + U*C*V*⅟A - U*C*(⅟C + V*⅟A*U)*⅟(⅟C + V*⅟A*U)*V*⅟A := by
+      congr
+      simp only [Matrix.mul_add, Matrix.mul_invOf_cancel_right, ← Matrix.mul_assoc]
+    _ = 1 := by
+      rw [Matrix.mul_invOf_cancel_right]
+      abel
+
+/-- Like `add_mul_mul_invOf_mul_eq_one`, but with multiplication reversed. -/
+lemma add_mul_mul_invOf_mul_eq_one' :
+    (⅟A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A)*(A + U*C*V) = 1 := by
+  calc
+    (⅟A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A)*(A + U*C*V)
+    _ = ⅟A*A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A*A + ⅟A*U*C*V - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A*U*C*V := by
+      simp_rw [add_sub_assoc, _root_.mul_add, _root_.sub_mul, Matrix.mul_assoc]
+    _ = (1 + ⅟A*U*C*V) - (⅟A*U*⅟(⅟C + V*⅟A*U)*V + ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A*U*C*V) := by
+      rw [invOf_mul_self, Matrix.invOf_mul_cancel_right]
+      abel
+    _ = 1 + ⅟A*U*C*V - ⅟A*U*⅟(⅟C + V*⅟A*U)*(V + V*⅟A*U*C*V) := by
+      rw [sub_right_inj, Matrix.mul_add]
+      simp_rw [Matrix.mul_assoc]
+    _ = 1 + ⅟A*U*C*V - ⅟A*U*⅟(⅟C + V*⅟A*U)*(⅟C + V*⅟A*U)*C*V := by
+      congr 1
+      simp only [Matrix.mul_add, Matrix.add_mul, ← Matrix.mul_assoc,
+        Matrix.invOf_mul_cancel_right]
+    _ = 1 := by
+      rw [Matrix.invOf_mul_cancel_right]
+      abel
+
+/-- If matrices `A`, `C`, and `C⁻¹ + V * A⁻¹ * U` are invertible, then so is `A + U * C * V`-/
+def invertibleAddMulMul : Invertible (A + U*C*V) where
+  invOf := ⅟A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A
+  invOf_mul_self := add_mul_mul_invOf_mul_eq_one' _ _ _ _
+  mul_invOf_self := add_mul_mul_invOf_mul_eq_one _ _ _ _
+
+/-- The **Woodbury Identity** (`⅟` version). -/
+theorem invOf_add_mul_mul [Invertible (A + U*C*V)] :
+    ⅟(A + U*C*V) = ⅟A - ⅟A*U*⅟(⅟C + V*⅟A*U)*V*⅟A := by
+  letI := invertibleAddMulMul A U C V
+  convert (rfl : ⅟(A + U*C*V) = _)
+
+end Woodbury
+
+end Ring
+
 end Matrix

Mathlib/LinearAlgebra/Matrix/NonsingularInverse.lean

@@ -607,6 +607,24 @@ theorem inv_diagonal (v : n → α) : (diagonal v)⁻¹ = diagonal (Ring.inverse
 
 end Diagonal
 
+section Woodbury
+
+variable [Fintype m] [DecidableEq m]
+variable (A : Matrix n n α) (U : Matrix n m α) (C : Matrix m m α) (V : Matrix m n α)
+
+/-- The **Woodbury Identity** (`⁻¹` version). -/
+theorem add_mul_mul_inv_eq_sub (hA : IsUnit A) (hC : IsUnit C) (hAC : IsUnit (C⁻¹ + V * A⁻¹ * U)) :
+    (A + U * C * V)⁻¹ = A⁻¹ - A⁻¹ * U * (C⁻¹ + V * A⁻¹ * U)⁻¹ * V * A⁻¹ := by
+  obtain ⟨_⟩ := hA.nonempty_invertible
+  obtain ⟨_⟩ := hC.nonempty_invertible
+  obtain ⟨iAC⟩ := hAC.nonempty_invertible
+  simp only [← invOf_eq_nonsing_inv] at iAC
+  letI := invertibleAddMulMul A U C V
+  simp only [← invOf_eq_nonsing_inv]
+  apply invOf_add_mul_mul
+
+end Woodbury
+
 @[simp]
 theorem inv_inv_inv (A : Matrix n n α) : A⁻¹⁻¹⁻¹ = A⁻¹ := by
   by_cases h : IsUnit A.det