I am implementing a phylogenetic library in Scala and stumbled over numerical problems with the GTR-model (Yang 1994). The transition probabilities after a given time t of a rate-matrix Q are given by the matrix exponential exp(Qt). According to the literature, this is often computed using the Eigen-Decomposition (V * D * V^-1) of Q. While this decomposition is numerical stable, the exponential (V * exp(Dt) * V^-1) is not. For increasing values of t (e.g. t > 10) the probabilities do not sum to 1.

This problem holds true for standard numerical platforms like octave and several java-linalg-frameworks. In mathematica however, the result is computed correctly. Does mathematica use and arbitrary precision arithmetic?

Does someone experience the same problems or has any idea how to handle it (e.g. rescaling of probabilities to 1)?