Hyperbolic tangents (of matrix) numerically inaccurate?

General Usage
1

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

2

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
3

Thanks for your input.
The scaling and squaring looks interesting, I will have a look into this.

4

That sounds like the way to go. I didn’t look carefully enough at the Schur-Parlett code, which looks like a stalled work in progress implementation of the algorithm from this paper. It did look fairly complete, but the author didn’t release it and I didn’t see a license or tests. That does seem to be a gap in Julia support for numerical linear algebra. That algorithm should probably be the default choice for matrix functions of non-normal matrices for functions for which there are not well established alternative algorithms.