Julia v1.11.1 with StaticArrays v1.9.8.
Discovered while attempting
julia> using StaticArrays
julia> M = @SMatrix [0 0 1f1; 0 0 0; 0 0 0]
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
0.0 0.0 10.0
0.0 0.0 0.0
0.0 0.0 0.0
julia> exp(M)
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
NaN NaN NaN
NaN NaN NaN
NaN NaN NaNThis is a duplicate of #785 up to this point, but some investigation reveals the underlying issue, which I discuss here.
The problematic calculation is in _exp where
julia> U = @SMatrix Float32[0.0 0.0 1.6191188f17; 0.0 0.0 0.0; 0.0 0.0 0.0]
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
0.0 0.0 1.61912f17
0.0 0.0 0.0
0.0 0.0 0.0
julia> V = @SMatrix Float32[6.4764752f16 0.0 0.0; 0.0 6.4764752f16 0.0; 0.0 0.0 6.4764752f16]
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
6.47648f16 0.0 0.0
0.0 6.47648f16 0.0
0.0 0.0 6.47648f16
julia> VmU = V - U; VpU = V + U;
julia> VmU \ VpU
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
NaN 0.0 NaN
0.0 NaN 0.0
0.0 0.0 NaNThis can be resolved by
julia> lu(VmU) \ (VpU)
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
1.0 0.0 5.0
0.0 1.0 0.0
0.0 0.0 1.0The issue is with _solve(::Size{(3,3)}, ::Size{(3,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}. The dependence of _solve on det makes it vulnerable to under/overflow. We can replicate the issue in Float64
julia> @SMatrix([1e200 0 0; 0 1e200 0; 0 0 1e200]) \ @SVector([1e0,0,0])
3-element SVector{3, Float64} with indices SOneTo(3):
NaN
NaN
NaNThis issue is closely related to #959, but that issue has to do with numerical instability and not overflow. Although it is likely that finding a resolution to that would improve the situation here. This issue is most acute in the 3x3 case, but also affects the 2x2 specialization.
I see that _solve is much faster than lu
julia> using BenchmarkTools
julia> @btime \($VmU, $VpU);
5.466 ns (0 allocations: 0 bytes)
julia> @btime \(lu($VmU), $VpU);
48.222 ns (0 allocations: 0 bytes)
but to what extent is this worth it? Is there a way we can keep the performance of the closed-form inverse without overflow? For example, can we scale the values in an attempt to keep det from overflowing? Can lu be made faster? At the very least, it seems that using lu in _exp would close #785.
Julia v1.11.1 with StaticArrays v1.9.8.
Discovered while attempting
This is a duplicate of #785 up to this point, but some investigation reveals the underlying issue, which I discuss here.
The problematic calculation is in
_expwhereThis can be resolved by
The issue is with
_solve(::Size{(3,3)}, ::Size{(3,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}. The dependence of_solveondetmakes it vulnerable to under/overflow. We can replicate the issue inFloat64This issue is closely related to #959, but that issue has to do with numerical instability and not overflow. Although it is likely that finding a resolution to that would improve the situation here. This issue is most acute in the 3x3 case, but also affects the 2x2 specialization.
I see that
_solveis much faster thanlubut to what extent is this worth it? Is there a way we can keep the performance of the closed-form inverse without overflow? For example, can we scale the values in an attempt to keep
detfrom overflowing? Canlube made faster? At the very least, it seems that usingluin_expwould close #785.