[BUG] `exp` of `SMatrix` fails numerically

[juliohm]

MWE:

julia> using StaticArrays

julia> M = 800.0 * [-1 1; 1 -1]
2×2 Matrix{Float64}:
 -800.0   800.0
  800.0  -800.0

julia> exp(M)
2×2 Matrix{Float64}:
 0.5  0.5
 0.5  0.5

julia> S = SMatrix{2,2}(M)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 -800.0   800.0
  800.0  -800.0

julia> exp(S)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 NaN  NaN
 NaN  NaN

[mikmoore]

It looks like the method is described in Corollary 2.4 of this paper.

It seems a fix will require re-doing the algebra. Check my math, but every term appears to be of the form

$$\exp(x) \left(a \cosh(y) + b \sinh(y)\right) = \frac{a + b}{2} \exp(x + y) + \frac{a - b}{2} \exp(x - y)$$

Using the right-hand version should avoid the range issues of computing the exp and hyperbolics separately. It should be a pretty simple PR.

[juliohm]

That is precisely the method, thank you for sharing @mikmoore.

I am trying to work my way to reach the result you shared, but no success so far. Could you perhaps share the steps so that we could both be sure of the derivation before a PR is submitted?

Which identities are used in your derivation?

[juliohm]

From the linked paper, the terms have a $\delta$ in the denominator though, so can't be expressed as $a\cosh(y) + b\sinh(y)$, see

That would be $a\cosh(y) + b\sinh(y) / y$

[mikmoore]

I just used the definitions $\cosh(x) = (\exp(x) + \exp(-x))/2$ and $\sinh(x) = (\exp(x) - \exp(-x))/2$. The equation above simply regroups the $\exp$'s from the hyperbolic functions with the leading $\exp$. The goal was to permit them to cancel to avoid the overflow issue you were seeing.

Sorry for recycling $a$ and $b$ in my statement when they were already used for matrix entries. Maybe that made it confusing. The (1,1) term uses $a = 1$, $b = \frac{a-d}{2\delta}$, $x = \frac{a+d}{2}$, $y=\delta$.

[juliohm]

You used the identity $\cosh(x) = (\exp(x) + \exp(-x)) / 2$. But isn't that missing the imaginary term?

The identities I know are displayed below:

Reference: https://hepweb.ucsd.edu/ph110b/110b_notes/node49.html

[mikmoore]

If $\cosh(ix) = \frac{\exp(ix) + \exp(-ix)}{2}$ then $\cosh(x) = \frac{\exp(x) + \exp(-x)}{2}$. The definitions you showed are intended to show the relationship between the trigonometric and hyperbolic functions, which is why they use imaginary arguments to the hyperbolics. For simpler definitions, see https://en.wikipedia.org/wiki/Hyperbolic_functions#Exponential_definitions.

[juliohm]

Thank you. I thought that the identity only worked with the imaginary part. I will try to check the expression you shared above and submit a PR.