Woodbury identity

Stream: mathlib4

[Kim Morrison]

Ahmad Alkhalawi said:

Woodbury identity

Hopefully you are planning on proving it in an arbitrary ring? Hopefully the proof is quite short, using e.g. the "algebraic proof" on the wikipedia page and the ring tactic.

[Ahmad Alkhalawi]

Kim Morrison said:

Ahmad Alkhalawi said:

Woodbury identity

Hopefully you are planning on proving it in an arbitrary ring? Hopefully the proof is quite short, using e.g. the "algebraic proof" on the wikipedia page and the ring tactic.

I have a proof for it in matrices over a commutative ring. Unfortunately I don't think ring helps since matrices are noncommutative, and abel/noncomm_ring helped but not much. I followed the direct proof in Wikipedia and it was about 22 lines.

[Kim Morrison]

Ahmad Alkhalawi said:

I have a proof for it in matrices over a commutative ring. Unfortunately I don't think ring helps since matrices are noncommutative, and abel/noncomm_ring helped but not much. I followed the direct proof in Wikipedia and it was about 22 lines.

But why restrict the result to matrices? From what I understand it's just true in any ring. Indeed I am wrong however about the ring tactic being helpful.

[Ahmad Alkhalawi]

Kim Morrison said:

Ahmad Alkhalawi said:

I have a proof for it in matrices over a commutative ring. Unfortunately I don't think ring helps since matrices are noncommutative, and abel/noncomm_ring helped but not much. I followed the direct proof in Wikipedia and it was about 22 lines.

But why restrict the result to matrices? From what I understand it's just true in any ring. Indeed I am wrong however about the ring tactic being helpful.

I proved it over general rings as well, but it seems like that doesnt imply the matrix version since the matrices don't have to be the same size. Maybe I could submit both?

[MohanadAhmed]

Ahmad Alkhalawi said:

Kim Morrison said:

Ahmad Alkhalawi said:

I have a proof for it in matrices over a commutative ring. Unfortunately I don't think ring helps since matrices are noncommutative, and abel/noncomm_ring helped but not much. I followed the direct proof in Wikipedia and it was about 22 lines.

But why restrict the result to matrices? From what I understand it's just true in any ring. Indeed I am wrong however about the ring tactic being helpful.

I proved it over general rings as well, but it seems like that doesnt imply the matrix version since the matrices don't have to be the same size. Maybe I could submit both?

Hi Ahamd,
Can you share both versions for the woodbury identity? I am interested in seeing how different they are.

[Ahmad Alkhalawi]

MohanadAhmed said:

Ahmad Alkhalawi said:

Kim Morrison said:

Ahmad Alkhalawi said:

I have a proof for it in matrices over a commutative ring. Unfortunately I don't think ring helps since matrices are noncommutative, and abel/noncomm_ring helped but not much. I followed the direct proof in Wikipedia and it was about 22 lines.

But why restrict the result to matrices? From what I understand it's just true in any ring. Indeed I am wrong however about the ring tactic being helpful.

I proved it over general rings as well, but it seems like that doesnt imply the matrix version since the matrices don't have to be the same size. Maybe I could submit both?

Hi Ahamd,
Can you share both versions for the woodbury identity? I am interested in seeing how different they are.

Hello, this is the matrix version:

theorem add_mul_inv
  (A: Matrix n n R) (U : Matrix n m R) (C : Matrix m m R) (V : Matrix m n R)
  [Invertible A] [Invertible C][Invertible (C⁻¹ + V * A⁻¹* U)] :
  (A+U*C*V)⁻¹=A⁻¹-A⁻¹*U*(C⁻¹+V*A⁻¹*U)⁻¹*V*A⁻¹:= by
    refine Matrix.inv_eq_right_inv ?h
    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
        rw[Matrix.add_mul, Matrix.mul_sub,Matrix.mul_sub]
        repeat rw[←Matrix.mul_assoc]
        rw[←@add_sub_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
        simp
        abel
      _ = 1+U*C*V*A⁻¹-(U+U*C*V*A⁻¹*U)*(C⁻¹+V*A⁻¹*U)⁻¹*V*A⁻¹ := by
        simp
        rw [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
        rw [Matrix.mul_add]
        simp
        rw [←Matrix.mul_assoc,←Matrix.mul_assoc]
      _ = 1 := by simp

Note R is commutative
This is the ring version

theorem add_mul_inv
  (A: S) (U : S) (C : S) (V : S)
  [Invertible A] [Invertible C][Invertible (⅟C + V * ⅟A* U)]:
  (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
        rw[add_mul, mul_sub,mul_sub]
        repeat rw[←mul_assoc]
        rw[←@add_sub_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
        simp
        abel
      _ = 1+U*C*V*⅟A-(U+U*C*V*⅟A*U)*⅟(⅟C+V*⅟A*U)*V*⅟A := by
        simp
        rw [add_mul,add_mul,add_mul]
      _ = 1+U*C*V*⅟A-U*C*(⅟C+V*⅟A*U)*⅟(⅟C+V*⅟A*U)*V*⅟A := by
        congr
        rw [mul_add]
        rw [mul_assoc U C ⅟C]
        simp
        rw [←mul_assoc,←mul_assoc]
      _ = 1 := by simp

Here \frac1 instead of \inv since \inv isnt defined, and i only prove it is a left inverse, since im not sure its two-sided

[Notification Bot]

11 messages were moved here from #mathlib4 > github permission by Eric Wieser.

[MohanadAhmed]

Thanks @Ahmad Alkhalawi

I tried to golf it but I could not do much. Compared to my earlier proof, I think this is much neater.
@Eric Wieser Any golfs you can think of?

import Mathlib.LinearAlgebra.Matrix.NonsingularInverse

theorem add_mul_inv {n m R: Type*}[CommRing R]
  [Fintype n][Fintype m][DecidableEq m][DecidableEq n]
  (A: Matrix n n R) (U : Matrix n m R) (C : Matrix m m R) (V : Matrix m n R)
  [Invertible A] [Invertible C][Invertible (C⁻¹ + V * A⁻¹* U)] :
  (A+U*C*V)⁻¹=A⁻¹-A⁻¹*U*(C⁻¹+V*A⁻¹*U)⁻¹*V*A⁻¹:= by
    refine Matrix.inv_eq_right_inv ?_
    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, Matrix.add_mul, Matrix.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 [Matrix.mul_inv_of_invertible, 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_inv_cancel_right_of_invertible, ←Matrix.mul_assoc]
      _ = 1 := by simp

[Eric Wieser]

Not a golf, but it would be great to show the inverse for the other side too, and bundle the resulting two-sided inverse up as Invertible (A + U*C*V)

[Eric Wieser]

(I guess the right-inverse follows quickly from the left-inverse)

[MohanadAhmed]

Should that not follow from: docs#Matrix.mul_eq_one_comm

[MohanadAhmed]

Eric Wieser said:

Not a golf, but it would be great to show the inverse for the other side too, and bundle the resulting two-sided inverse up as Invertible (A + U*C*V)

Is this what you meant?

lemma add_mul_mul_inv_mul_eq_one {n m R: Type*}[CommRing R]
  [Fintype n][Fintype m][DecidableEq m][DecidableEq n]
  (A: Matrix n n R) (U : Matrix n m R) (C : Matrix m m R) (V : Matrix m n R)
  [Invertible A] [Invertible C][Invertible (C⁻¹ + V * A⁻¹* U)] :
  (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, Matrix.add_mul, Matrix.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 [Matrix.mul_inv_of_invertible, 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_inv_cancel_right_of_invertible, ←Matrix.mul_assoc]
      _ = 1 := by simp

noncomputable def invertibleWoodburyLHS {n m R: Type*}[CommRing R]
  [Fintype n][Fintype m][DecidableEq m][DecidableEq n]
  (A: Matrix n n R) (U : Matrix n m R) (C : Matrix m m R) (V : Matrix m n R)
  [Invertible A] [Invertible C] [Invertible (C⁻¹ + V * A⁻¹* U)]  : Invertible (A+U*C*V) := by
  apply Matrix.invertibleOfRightInverse _ _ (add_mul_mul_inv_mul_eq_one A U C V)

theorem woodbury_identity {n m R: Type*}[CommRing R]
  [Fintype n][Fintype m][DecidableEq m][DecidableEq n]
  (A: Matrix n n R) (U : Matrix n m R) (C : Matrix m m R) (V : Matrix m n R)
  [Invertible A] [Invertible C][Invertible (C⁻¹ + V * A⁻¹* U)] :
  (A+U*C*V)⁻¹=A⁻¹-A⁻¹*U*(C⁻¹+V*A⁻¹*U)⁻¹*V*A⁻¹:=
  Matrix.inv_eq_right_inv (add_mul_mul_inv_mul_eq_one _ _ _ _)

[MohanadAhmed]

@Ahmad Alkhalawi Did you get non master mathlib access yet?
I think we should pr this before it becomes mathlib outdated like my earlier code

[MohanadAhmed]

@Eric Wieser Do you have any more thoughts on this?

[Ahmad Alkhalawi]

MohanadAhmed said:

Ahmad Alkhalawi Did you get non master mathlib access yet?
I think we should pr this before it becomes mathlib outdated like my earlier code

I did not

[MohanadAhmed]

@maintainers Maintainers, can someone please @Ahmad Alkhalawi (github name 4hma4d access to PR his proof of Woodbury (aka Sherman Morrison) identity to mathlib?

[Matthew Ballard]

@Ahmad Alkhalawi invitation sent!

[Ahmad Alkhalawi]

Thanks!

[Eric Wieser]

Kim Morrison said:

But why restrict the result to matrices? From what I understand it's just true in any ring.

Ahmad Alkhalawi said:

I proved it over general rings as well, but it seems like that doesnt imply the matrix version since the matrices don't have to be the same size. Maybe I could submit both?

I think the correct generality here is over a ring equipped with additive maps, trading the key steps using

(U*C*V)*A*(U*B*V) = U*(C*(V*A*U)*B)*V

for

Fuv C*A*Fuv B = Fuv (C * Fvu A * B)

where on matrices all the * are still square multiplications. This can be instantiated separately for Fuv C = U*C*V for regular and matrix multiplication.

I think @Ahmad Alkhalawi's #16325 can get merged without that generalization though, assuming everything else looks fine.