Numeric stability of Base.exp of static matrix of Float32

General Usage
, , , ,
1

Guys,

I want to use Base.exp in my model which uses Float32 with StaticArray and StaticVector. However, I found there is a numerical instability of Base.exp with the 3*3 matrix…see the example above:

using StaticArrays

julia> deriv_mat = Float32[-50.75532 69.86392 0.0; 50.75532 -139.72784 54.483025; 0.0 69.86392 -108.96605]
3×3 Matrix{Float32}:
 -50.7553    69.8639     0.0
  50.7553  -139.728     54.483
   0.0       69.8639  -108.966

julia> exp(deriv_mat)
3×3 Matrix{Float32}:
 8.57472f-6  6.68494f-6  3.72485f-6
 4.85653f-6  3.78621f-6  2.10967f-6
 3.46999f-6  2.70524f-6  1.50736f-6

julia> deriv_mat_s = SMatrix{3, 3, Float32}(Float32[-50.75532 69.86392 0.0; 50.75532 -139.72784 54.483025; 0.0 69.86392 -108.96605])
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
 -50.7553    69.8639     0.0
  50.7553  -139.728     54.483
   0.0       69.8639  -108.966

julia> exp(deriv_mat_s)
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
 NaN  NaN  NaN
 NaN  NaN  NaN
 NaN  NaN  NaN

I need to use STRICTKY the Float32 with static matrix…so is there any way to overcome such numerical instability? if I convert it to Float64 then it creates extra allocations…

I am thinking to use something like this…@btime SMatrix{3, 3, Float32}(exp(SMatrix{3, 3, Float64}($deriv_mat_s_64))) but very ugly…

A3 = rand(3, 3)
As3 = SMatrix{3, 3}(A3)

julia> @btime exp($As3)
  93.400 ns (0 allocations: 0 bytes)
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.52765   0.53686  0.282278
 0.792425  1.74997  0.504976
 1.28887   1.30304  1.93264

julia> @btime SMatrix{3, 3, Float32}(exp(SMatrix{3, 3, Float64}($deriv_mat_s_64)))
  113.689 ns (0 allocations: 0 bytes)
3×3 SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
 8.5747f-6   6.68493f-6  3.72484f-6
 4.85652f-6  3.7862f-6   2.10967f-6
 3.46999f-6  2.70524f-6  1.50736f-6
2

I don’t know about the numerical instability, but this is less ugly:

Float32.(exp(Float64.(x)))
3

Possibly related to [BUG] exp of SMatrix fails numerically · Issue #1295 · JuliaArrays/StaticArrays.jl, though the 2x2 example there works now.

Just checking, do you intend exponentiation of each element or the generalization to square matrices?

4

not the element-wise exponential, but the matrix exponential.

5

See issues #785 (which first reports the issue) and #1285 (which correctly diagnoses the problem).

This is exacerbated by the large coefficients used in the Pade approximation of StaticArrays’ version of exp. One improvement would be to scale all the coefficients down so that the largest (rather than smallest) ones are on the order of 1 (in fact, it would probably resolve most issues with overflow in exp). It wouldn’t be a very difficult pull request, just find each set of coefficients and scale them by 2^{-n} for a useful choice of n.

2 Likes
6

Not strictly in Base as you know by now, and I filed an issue just now:

but I see it’s a known problem, and seemingly the third duplicate, so I closed mine. But at least see the links to the code from mine. I didn’t move them to other issues.