# Efficient ways to numerically evaluate matrix exponentials

What are some computationally efficient ways to solve matrix exponentials, i.e. functions of the form : $$f(X)=e^{X}$$, where $$X$$ is a square matrix?

So far I have been able to diagonalise some matrices and find the exponent of individual diagonal elements, but not all matrices I'm dealing with will be diagonalisable.

I am using Python with SciPy/NumPy, so solutions that can be implemented here will be most useful. If not, general solutions/solutions from other platforms are welcome too.

Notes:

• I need the exponential itself, not a solution using it.
• The matrix $$X$$ is dense, typically small ($$3\times 3$$ or $$4\times 4$$), might not be symmetric or Hermitian.
• For an overview of numerical methods, check out the "19 dubious ways..." paper: cs.cornell.edu/cv/ResearchPDF/19ways+.pdf – carlosvalderrama Jan 17 '20 at 15:55
• Do I understand your question correctly by saying: You have matrix $X$ and now you want to compute $e^{X}$ right? Have you looked at scipy.linalg.expm? It should pretty efficient I believe. What is the size of $X$? Is it dense or sparse? Is it symmetric/Hermitian/etc.? – Alone Programmer Jan 17 '20 at 15:55
• Do you need to actually compute the matrix exponential or do you only need matrix-vector products with it? If you only need the latter, then there are some efficient Krylov method for doing this that can also exploit sparsity epubs.siam.org/doi/abs/10.1137/… – whpowell96 Jan 17 '20 at 16:06
• @VoB- The matrices here represent deformation gradients of a deforming crystal. The change in distance between sets of parallel planes due to deformation is given by the individual elements of the inverted matrix exponential. Further information is obtained from more element-wise operations. I think for accessing the elements as well as inverting the matrix, we need the entire matrix. – ScientificPythonNovice Jan 17 '20 at 20:14
• In that case, one relatively trivial way to improve part of your computation is to compute the exponential of $-A$ instead of computing the exponential then inverting. If you need the inverse to solve linear systems, then a Krylov-based method would probably still be more efficient depending on the problem structure since $(e^A)^{-1} = e^{(-A)}$ – whpowell96 Jan 17 '20 at 20:25

As I mentioned in my comment, due to that you are searching for a method based on Python and ideally available in NumPy or SciPy, my suggestion is to use scipy.linalg.expm. It's really easy to use basically you have the $$X$$ matrix and you just pass it to the scipy.linalg.expm class and it would give you the exponential of $$X$$ (i.e. $$e^{X}$$).
Due to that you did not specify any specific feature of $$X$$, such as its sparsity or its symmetry, etc., it's really difficult to say when scipy.linalg.expm will fail to return the correct exponential of $$X$$. But as far as I know the the time complexity of algorithm used in scipy.linalg.expm is $$O(n^{3})$$, so, I expect if you have a really big dense matrix, it would be impractical to calculate its exponential and as other people are trying to say it would be better idea if you really want it to multiply this exponential to a vector and use other methods, which might be more efficient.
• Commenting just for the record: the main factor to assess the accuracy of this Padé method is the norm of $X$. If you check in Section 10.3 of Higham's book Functions of matrices, you will see explicit bounds such as the following: if $\|X\| \leq 5.4$, then the (13,13) Padé approximant $r(X)$ is such that $\|r(X)-\exp(X)\| / \|\exp(X)\|$ is smaller than machine precision. – Federico Poloni Apr 11 '20 at 9:43