# Hyperbolic tangents (of matrix) numerically inaccurate?

Hello everyone,

I need to calculate the tanh of a complex square matrix. The function needs to work for different matrixes and the numbers can get relative big (well, values around 20 don’t seem big but apparently they are big enough to cause trouble). E.g.

``````C = ComplexF64[20.532918596935836 + 13.825923408965787im 18.707372202499215 - 13.100775315493328im; 21.669638013670312 - 15.16565944942833im 19.776899867610187 + 14.368130886743469im]
``````

The standard version gives me an error `ERROR: SingularException(2)`.

Looking at the https://github.com/JuliaLang/julia/blob/ed4c44fbafbbace523ebdc9440bf1f542157048a/stdlib/LinearAlgebra/src/dense.jl#L1057 I figured this is essentially

``````tanh2(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
``````

which results in the same error.
This lead me to two alternatives:

``````tanh3(x) = (exp(2x) - I) / (exp(2x) + I)
tanh4(x) = (I - exp(-2x)) / (I + exp(-2x))
``````

which work both but give (as expected) very different results for C.
I expected Version 4 to be right in this case as C has a quite big positive real part. For the verification, I used

``````tanhdf(x) = tanh(ComplexDF64.(x))
``````

with DoubleFloats (and GenericLinearAlgebra that seems to be a missing dependency in the latest version) for verification and got similar results.

My questions:

• Has someone already implemented a better version of tanh somewhere that has higher accuracy so that I can simply use this version?
• I’m not sure of when to choose which version (tanh2,3,4). I would have chosen based on the common factor of real(C) (In this case approx. 20) bigger some threshhold → 4, negative threshold → 3 original version otherwise.

Thanks

I don’t think this is an entirely solved problem in which you can expect to find existing code that is flawlessly reliable in all cases. With that said, using the naive exponential formulas directly to produce intermediate results with large norm does look like it’s asking for trouble. If I needed to compute this, I would try one of two things:

• You might consider a Schur-Parlett approach. There is a Julia implementation with blocking at Funm.jl.

• A quick web search turned up this paper. One of the considered methods is to scale the matrix first, compute the exponential for a scaled matrix with relatively small norm (which should be easy) and then apply some identities to recover the hyperbolic tangent of the original matrix. It is analogous to scaling and squaring for the matrix exponential, which is a standard approach, but it also leverages existing code for computing matrix exponentials. If you want to stick to using exponentials, it would probably improve things. (The authors do report better results with scaling and direct Taylor approximation of the hyperbolic tangent, however. But that looks like a little more work to try.)

1 Like