Matrix exponentiation with scipy: expm, expm2 and expm3

Question

Matrix exponentiation can be performed in python using functions within the scipy.linalg library, namely expm, expm2, expm3. expm makes use of a Pade approximation; expm2 uses the eigenvalue decomposition approach and expm3 makes use of a Taylor series with a default number of terms of 20.

In SciPy 0.13.0 release notes it is stated that:

The matrix exponential functions scipy.linalg.expm2 and scipy.linalg.expm3 are deprecated. All users should use the numerically more robust scipy.linalg.expm function instead.

Although expm2 and expm3 are deprecated since release version SciPy 0.13.0, I have found that in many situations these implementations are faster than expm. From this, some questions arise:

In what situations could expm2 and expm3 result in numerical instabilities?

In what situations (e.g. sparse matrices, symmetric, ...) is each of the algorithms faster/more precise?

Answer 1

This will depend a lot on the detail of the implementation of these different ways of exponentiating the matrix.

In general terms, I would expect the eigen-decomposition (expm2) to be poorly suited to sparse matrices, because it is likely to remove the sparseness. It will also be more difficult to apply to non-symmetric matrices, because this will require the use of complex arithmetic and more expensive algorithms to compute the eigen-decomposition.

For the Taylor-series approach (expm3), this sounds risky if there are a fixed number of terms independent of the norm of the matrix. When computing e^x for a scalar x, the largest terms in the Taylor series are around that for which n is close to x.

However, the implementation details of these (deprecated) functions may use tricks like diagonally loading the matrix so as to improve the stability of these series expansion.