Julia Community 🟣

Cover image for Exploring the Matrix Inversion Lemma in Julia
Vinod V
Vinod V

Posted on

Exploring the Matrix Inversion Lemma in Julia

Let   AA be a square matrix of size  nn  that is invertible, and let  uu  and vv   be non-zero column vectors of size nn  . The matrix  A+uvTA + u v^T  is invertible if and only if the equation  vTA1u1v^TA^{−1}u \neq −1 is true. In such a case, the inverse of   A+uvTA + u v^T  can be determined using the following method:

(A+uvT)1=A1A1uvTA11+vTA1u(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^T A^{-1}}{1+v^TA^{−1}u}

The above formula can be coded in Julia as

using LinearAlgebra

N = 4
A = rand(N,N)
Ainv = inv(A)

println("A = ")
display(Ainv)

u = rand(N,1)
v = rand(N,1)

B = A + u*v'

println("A + u*v' = ")
display(B)

Binv = inv(B)

if !isequal(-1,dot(v,Ainv*u))
    C = Ainv - Ainv*u*v'*Ainv/(1 + dot(v,Ainv*u))
end

println("Crude Inverse")
display(Binv)

println("Inverse calculated with Matrix Inversion Lemma")
display(C)

println("Are the matrices equal")
isapprox(Binv,C)
Enter fullscreen mode Exit fullscreen mode

Top comments (0)