Square root of the inverse of a matrix

Specific Domains Numerics
6

So to put this all together, here’s how I understand this should be approached.

You can use your favorite factorization as it suits you. I tested svd and eigen and found svd was a bit slower, but as Antoine said, it can be more stable:

using Random, RandomMatrices
function f1(A, Y, R)
    (op1=A\Y, 
     op2=real.(sqrt(inv(A)))*R)
end
function f2(A::Hermitian, Y, R)
    (op1 = A\Y, 
     op2 = sqrt(A)\R)
end
function f3(A::Hermitian, Y, R)
    λ, ϕ = eigen(A)
    
   (op1 = ϕ * Diagonal(1 ./ λ) * ϕ' * Y,
    op2 = ϕ * Diagonal(1 ./ sqrt.(λ)) * ϕ' * R)
end

let N = 100, Nx = 50
    U = randn(Nx, N) #switched Nx and N so the code runs
    A = I + U'*U
    
    Y = randn(N, Nx)  #switched Nx and N so the code runs

    # Generate Random Matrices
    d = Haar(1)
    R = rand(d, N)

    r1 = @btime f1($A, $Y, $R)
    r2 = @btime f2($(Hermitian(A)), $Y, $R)
    r3 = @btime f3($(Hermitian(A)), $Y, $R)
    Tuple(r1) .≈ Tuple(r2) .≈ Tuple(r3) #test that all the outputs are the same
end

#+RESULTS:
:   9.683 ms (31 allocations: 1.38 MiB)
:   1.781 ms (33 allocations: 705.03 KiB)
:   1.149 ms (31 allocations: 704.16 KiB)
: (true, true)
2 Likes