Computing Inverse of a stack of matrices

General Usage
,
2

If you have lots of tiny matrices, then StaticArrays is usually a good idea:

@btime compute_inv(A) setup=(A = randn(1000,4,3,3)); # 1.928 ms

function slice_inv(A) # expects matrix indices first
    B = similar(A)
    @inbounds for j in axes(A,4)
        for i in axes(A,3)
            B[:,:,i,j] .= inv(@view A[:,:,i,j]);
        end
    end
    B
end

@btime slice_inv(A) setup=(A = randn(3,3,4,1000)); # 1.260 ms

using StaticArrays

function slice_inv2(A::Array{T,4}) where {T}
    B = reinterpret(SArray{Tuple{3,3},T,2,9}, vec(A))
    C = map(inv, B)
    reshape(reinterpret(T, B), size(A))
end

@btime slice_inv2(A) setup=(A = randn(3,3,4,1000)); # 43.037 μs
3 Likes