7
$\begingroup$

Is there a way to compute one matrix element of the exponential of a tridiagonal matrix without having to compute the rest of the elements?

Motivation: I'm trying to find the first passage time distribution from a master equation. I can impose an absorbing boundary at the threshold $n$, and the master equation with the new boundary condition is of the form $\frac{dp}{dt}=A p$ for a tridiagonal $n\times n$ matrix $A$. Then the first passage time distribution can be written as a particular matrix element of $\exp(A\,t)$. It takes forever for Mathematica to compute $\exp(A\,t)$ for large $n$, so I was wondering if there is a way to compute only the desired matrix element and not the whole matrix.

Improve this question
$\endgroup$
2
  • 4
    $\begingroup$ Do you wish to compute your matrix element for one time $t$ or for many times $t$? If the enormous computation time in your problem stems from a combination of $n$ being large and of evaluating $e^{tA}$ for many $t$'s, it could be very helpful to compute the diaganalisation of $A$ first (if it exists and is numerically stable - for this it would be helpful, of course, if $A$ was for instance symmetric). In this case, it would take some time to compute the diagonalisation at first, but afterwards you can compute your single matrix element of $e^{tA}$ with $\mathcal{O}(n)$ for each time $t$. $\endgroup$
    – Jochen Glueck
    Feb 17, 2018 at 8:58
  • $\begingroup$ In this context, I think a natural thing to do when $n$ is large is to use Monte Carlo instead of numerical linear algebra. Specifically compute the first passage times by simulating the underlying birth-death process associated to the tridiagonal $A$. Although one may need a large number of samples to resolve resulting the Monte Carlo error, typically this number is independent of $n$. $\endgroup$
    – Nawaf Bou-Rabee
    Feb 17, 2018 at 20:10

1 Answer 1

Reset to default
12
$\begingroup$

Yes! Most methods to compute exponentials of large sparse matrices are based on computing directly $\exp(A)b$ for a given vector $b$ rather than the full matrix $\exp(A)$. Just take $b$ as a vector of the canonical basis, and you're set.

The basic idea is that these algorithms are based on approximating the exponential with a rational function which is then expanded into a sum of partial fractions $$\exp(A)b \approx q(A)^{-1}p(A)b = \sum_{i=1}^k \omega_k (A-\tau_k I)^{-1} b.$$

Hence the problem is reduced to solving several shifted linear systems of the form $(A-\tau_k I)x_k=b$; this should work particularly well for a tridiagonal matrix, for which linear systems can be solved cheaply.

This technique is often combined with ideas from Krylov subspace methods, especially rational Krylov methods.

See for instance a chapter in the review https://doi.org/10.1017/S0962492910000036, https://doi.org/10.1137/100788860, or http://onlinelibrary.wiley.com/doi/10.1002/gamm.201310002/full, as well as the Matlab code in Stefan Güttel's webpage which you can use to compute it.

All these algorithms are approximate (but in the end, what is not approximate on a computer?), so you may encounter trouble if the element that you wish to compute is many orders of magnitudes smaller than $\|\exp(A)\|$. In that case, you may be interested in a bound on the decay rate of off-diagonal entries in https://arxiv.org/abs/1501.07376.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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