StaticArrays.jl
StaticArrays.jl copied to clipboard
[Feature Request] Eigenvalues of small asymmetric matrices without allocation?
Recently I found myself trying to find the eigenvalues of 1 million 3x3 rotations, and I was very frustrated by how slow it was going, so I turned to StaticArrays. It looks like StaticArrays provides methods to speed this up for hermitian matrices, but not for asymmetric matrices like rotations, so I was out of luck. Would anyone else be interested in seeing this feature added to StaticArrays?
One obstacle that could be a problem is the return type. Eigenvalues of general real matrices can be either real or complex. The standard library handles this by having eigen
, eigvals
, and eigvecs
return Union types over complex and real, e.g. Union{Vector{Float64}, Vector{ComplexF64}}
for eigvals
. I think this is a problem in general because of runtime type reflection and unnecessary allocations due to the unknown size of the result. If I'm not mistaken this will especially be a problem if we want to add this feature to StaticArrays since we can't have a result of unknown size without allocating.
So there's a few things we can do:
- Have the eigen functions accept an argument for the scalar type in the return value, i.e.
eigvals(ComplexF64, A)
would returnVector{ComplexF64}
. We could make this argument mandatory, or we could have the type default to complex, or have it default to real and throw aDomainError
when the eigenvalues are complex. - We could implement a no-allocation eigen decomposition algorithm just for complex matrices, and continue to default to the standard library for general real matrices. This way someone can cast their matrix to a complex matrix when they don't expect their eigenvalues to be real and they don't want to pay for allocations and type reflection.
- Have separate methods based on the return type:
eigvalsreal
,eigvalscomplex
, etc.
I think the first option would make the most sense. What do people think? I'm also wondering if for static arrays larger than 3x3 there's a way we can call the LAPACK routines without allocating, which would speed things up for medium sized matrices as well. But anyways, this all sounds like a big project. Would people want these features?