6
$\begingroup$

Given an anti-Hermitian and sparse matrix $X$, I am using Python (NumPy and SciPy) to compute the matrix exponential $f(t) := e^{tX}$ for many values of $t$. The method I am currently using is to essentially compute the exponent with the following function:

X = SomeSparseMatrix
def f(t):
    return scipy.sparse.linalg.expm(t * X)

I'm wondering if there is a way to make this faster. For example, if I compute this function once to obtain $e^X$, can I obtain $e^{tX}$ from $e^X$ without having to recompute the exponent?

$\endgroup$
8
  • 3
    $\begingroup$ Does your matrix have any useful properties? The link below provides some ideas for solving your problem. Even powers of antihermitian matrices are hermitian, and odd powers are antihermitian. Perhaps you can do something with that? sites.millersville.edu/bikenaga/linear-algebra/… $\endgroup$ – Charlie S Dec 8 '20 at 21:21
  • 2
    $\begingroup$ If $t$ is an integer, you should be able to compute $e^{tX}$ using repeated squaring. Since $e^X$ is unitary, the terms of the matrix won't blow up $\endgroup$ – ogogmad Dec 8 '20 at 21:23
  • 3
    $\begingroup$ Have you searched for the two papers about "dubious ways to compute the matrix exponential"? $\endgroup$ – Wolfgang Bangerth Dec 8 '20 at 21:24
  • 3
    $\begingroup$ This is the paper that @WolfgangBangerth is referring to and I also highly recommend it: doi.org/10.1137/S00361445024180 $\endgroup$ – Daniel Shapero Dec 8 '20 at 23:09
  • 3
    $\begingroup$ Another question -- do you actually need the matrix exponential itself, or do you only need its action $e^{tX}\cdot v$ for some vector $v$? In many applications it's the latter and there are better ways to do that without explicitly computing the exponential. In general the exponential of a sparse matrix is dense, just like the inverse. If you can avoid explicitly forming this matrix then you should do so. $\endgroup$ – Daniel Shapero Dec 8 '20 at 23:53
9
$\begingroup$

An anti-Hermitian matrix is diagonalizable, with orthogonal eigenvectors (ref). Hence you can write $X = PDP^{-1}$, where $D$ is a diagonal matrix. Therefore the exponential can be calculated as $e^X=Pe^DP^{-1}$, and $e^{tX} = Pe^{tD}P^{-1}$.

If $d_1$, $d_2$, etc.... are the diagonal elements of $D$, then for each value of $t_i \in t$, $e^{t_iD}$ is just a diagonal matrix with elements $e^{t_id_1}$, $e^{t_id_2}$, etc....Now your calculation of the exponential for each value $t_i$ is just reduced to $N$ scalar exponentiations for an $N\times N$ matrix $X$.

You still have the task of diagonalizing $X$,but you only have to do it once. If $t$ has $M$ elements, then you also have $2M$ matrix multiplications to do.

So you would be trading $M$ matrix exponentiations for one diagonalization of $X$, as well as $N\cdot M$ scalar exponentiations, and $2M$ matrix multiplications. I don't know what the complexity of matrix exponentiation is, but I think this is likely to be a good trade.

As Daniel Shapero points out in his comment above, if you only need the action of $e^{tX}v$ for some vector $v$ and don't really need the matrix, then you can often save a whole lot of time. That savings is easy to realize with this approach, since $e^{tX}v = Pe^{tD}P^{-1}v=Pe^{tD}(P^{-1}v)$. You only do the $P^{-1}v$ calculation once, and the $Pe^{tD}$ calculations can be done faster than an ordinary matrix multiplication since $e^{t_iD}$ is a diagonal matrix. In this case, the diagonalization is the only thing that scales as $N^3$ (or maybe a bit less, depending on the implementation of the diagonalization), and it is done just once, not once per each point in $t$.

$\endgroup$

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.