Skip to content

Commit

Permalink
feat(LinearAlgebra/Matrix): Woodbury Identity (#16325)
Browse files Browse the repository at this point in the history
This adds the [Woodbury Identity](https://en.wikipedia.org/wiki/Woodbury_matrix_identity).

[Zulip discussion](https://leanprover.zulipchat.com/#narrow/stream/287929-mathlib4/topic/Woodbury.20identity/near/462245284)

Also corrects some bad deprecations introduced in #16590, which affected development of this PR.

Co-authored-by: Mohanad Ahmad



Co-authored-by: Eric Wieser <wieser.eric@gmail.com>
  • Loading branch information
2 people authored and fbarroero committed Oct 2, 2024
1 parent 70ed733 commit 1965992
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 5 deletions.
77 changes: 72 additions & 5 deletions Mathlib/Data/Matrix/Invertible.lean
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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 α)
Expand Down Expand Up @@ -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
18 changes: 18 additions & 0 deletions Mathlib/LinearAlgebra/Matrix/NonsingularInverse.lean
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 1965992

Please sign in to comment.