Hi, so today I found out that appearantly calling eigen on a hermitian 3x3 SMatrix{ComplexF64} allocates quite a bit of memory which I do not quite understand, since there are no allocations for a real matrix:

Yes, seems that you’re right: In the static arrays library it says:

# [...]
# TODO extend the method to complex hermitian
@inline function _eig(::Size{(3,3)}, A::LinearAlgebra.HermOrSym{T}, permute, scale) where {T <: Real}

Thats a little unfortunate, Ill try my luck with eigen!, although I believe it does not remove all allocations. Do you know a good way to get rid of allcations when repeatedly computing eigenvalues of many small matrices?

Some time ago I dealt with a similar problem, and found this paper (unpublished I think) describing a faster way of computing the eigendecomposition of 3x3 Hermitian matrices.

I wrote an implementation of this, feel free to use it:

Note that I only tested this with real Hermitian matrices, but it’s also supposed to work for complex ones. Maybe it’s worth adding this to StaticArrays, though I had small issues with numerical stability (I think) so I wasn’t sure it was a good idea. But that was probably because some of my matrices were badly conditioned…

thanks for that, I will check it out! I will probably still have have to use generic methods that allocate as little as possible since larger matrices could also occur

Nice move writing the eigenvalues in terms of (a-b)^2, etc. Seems like the authors of the paper should have done that if they were thinking about numerics at all. I don’t have a counter-example, but without an error analysis I wouldn’t trust the eigenvalue formulas for the 3\times 3 case. Even with your fix, it looks like there are a lot of places you could run into troubles.

It’s easy to demonstrate that the eigenvector formulas are unstable:

Both the measure or numerical orthogonality \|V^T V - I\| and the residual \|AV - VD\| should be on the order of the unit roundoff for a stable algorithm. Here, they are off that by 6 orders of magnitude. The eigenvectors are well conditioned in this case since the corresponding eigenvalues are well separated.