3
$\begingroup$

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.
Improve this question
$\endgroup$
9
  • 5
    $\begingroup$ For an overview of numerical methods, check out the "19 dubious ways..." paper: cs.cornell.edu/cv/ResearchPDF/19ways+.pdf $\endgroup$
    – carlosvalderrama
    Jan 17 '20 at 15:55
  • 3
    $\begingroup$ 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.? $\endgroup$
    – Mithridates the Great
    Jan 17 '20 at 15:55
  • 2
    $\begingroup$ 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/… $\endgroup$
    – whpowell96
    Jan 17 '20 at 16:06
  • 2
    $\begingroup$ @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. $\endgroup$
    – ScientificPythonNovice
    Jan 17 '20 at 20:14
  • 1
    $\begingroup$ 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)}$ $\endgroup$
    – whpowell96
    Jan 17 '20 at 20:25
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks! The matrix in question is not large (mostly 3x3), and is dense, so time-wise it shouldn't be a problem. My concern was whether the approximation used in this function would fail for certain (3x3) matrices, i.e. does it require them to be symmetric etc. The documentation doesn't indicate such a limitation, so I think it is best to try it out on the matrices I have and see if there are any exceptions. $\endgroup$
    – ScientificPythonNovice
    Jan 17 '20 at 20:19
  • 1
    $\begingroup$ 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. $\endgroup$
    – Federico Poloni
    Apr 11 '20 at 9:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.