8
$\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
    Commented 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
    Commented 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$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .