I’m currently using sqrt() from LinearAlgebra.jl to compute the square-root of 3 x 3 matrices repeatedly. But this has created a performance bottleneck. So I was wondering if there an implementation of a matrix square root that specializes for 3 x 3 matrices like the ones in StaticArrays.jl(which are limited to 2 x 2 matrices). FYI I’m dealing with real, symmetric and positive-definite matrices.

Depending on exactly what you need, you might be able to use the Cholesky decomposition, which computes a lower triangular matrix L or upper triangular matrix U such that L*L' == U'*U == A. I believe StaticArrays has specialized implementations for cholesky(A::SMatrix). This is similar to a sqrt in some ways, although not in the sense that S = sqrt(A) such that S*S == A (no transpose needed here and S is PSD if A is).

Right, it would be good to know the context here. In many problems where you conceptually might want a matrix square root, there is a way to reformulate it to use Cholesky, which is way cheaper to compute.

The Symmetric LinearAlgebra sqrt is implemented with eigen, so with a Symmetric SMatrix you don’t have to build it yourself.

If speed is more critical than precision you can do considerably better than that, however. This function only works for sufficiently distinct eigenvalues, so needs to be complemented with some special cases for colliding eigenvalues.

I need to factorize a 3 x 3 matrix A such that B*B’ = A and then B is used elsewhere in the code. The problem with the Cholesky decomposition is that now I’m not sure that my A will always be positive definite. But it is true that Cholesky is much faster than computing the square-root even with StaticArrays and Symmetric.