Numerical issues with exp(::SMatrix{3, 3, Float32})

[willtebbutt]

MWE:

julia> using StaticArrays

julia> F = Float32[0.0 1.0 0.0; 0.0 0.0 1.0; -11.18034 -15.0 -6.708204]

# Works fine with `Matrix{Float32}`
julia> exp(F)
3×3 Array{Float32,2}:
  0.613059   0.345864   0.053439
 -0.597466  -0.188525  -0.0126152
  0.141042  -0.408237  -0.1039

julia> F_static = SMatrix{3, 3, Float32}(F);

# Breaks for `SMatrix{3, 3, Float32}`
julia> exp(F_static)
3×3 SArray{Tuple{3,3},Float32,2,9} with indices SOneTo(3)×SOneTo(3):
 NaN  NaN  NaN
 NaN  NaN  NaN
 NaN  NaN  NaN

Note that this matrix is a little weird (might be defective, I'm not actually sure).

It still happens if you pick a matrix that is diagonalisable:

julia> A = SMatrix{3, 3, Float32}(randn(3, 3));

julia> exp(A)
3×3 SArray{Tuple{3,3},Float32,2,9} with indices SOneTo(3)×SOneTo(3):
 NaN  NaN  NaN
 NaN  NaN  NaN
 NaN  NaN  NaN

Same thing if the matrix in question is symmetric.

Tested on version 0.12.3

[c42f]

So it appears this is particularly problematic with Float32, and choosing matrices with randn, for some reason.

The NaN's arise from the matrix solve at

StaticArrays.jl/src/expm.jl

Line 112 in e0fe5da

expA = (V - U) \ (V + U)

The obvious difference between Base and the StaticArrays implementation is the lack of the matrix balancing (the gebal call). This might be a place to start looking closer.