# Efficient computation of Markov chain transition probability matrix

Consider a continuous Markov chain $X = (X_t)$ on a finite state space and let $Q$ be the (given) transition rate matrix. This matrix is very sparse, with non-zero values on 3 diagonals only (so from each state, there can be transitions to 2 other states only).

Let $P_t$ be the transition probability matrix, so $P_t(j,k) =$ Prob($X_t = k | X_0 = j$).

My question is: what is the best way to quickly compute the $j$th row of $P_1$?

Solving the Kolmogorov forward equations gives $P_t = e^{Qt}$, so one method is to perform this computation explicitly in matlab: expm(Q). But I'm thinking that there is perhaps a better way, particularly given the structure of $Q$ and since I'm only interested in one row of $P_1$. The actual instance of the problem I'm solving is small (120 states, say), but I would like the computation to be very fast.

Edit: following the suggestion by @meij in the comments below, I posted this question here and received some useful answers: http://scicomp.stackexchange.com/questions/6828/efficient-computation-of-markov-chain-transition-probability-matrix

-
If $Q$ is diagonalizable, say $Q = V\Lambda V^{-1}$, then using $e^{tQ} = V e^{t\Lambda} V^{-1}$ seems likely to be faster than expm. You can also approximate the exponential using Krylov subspace (Lanczos or Arnoldi) techniques. –  Steve Huntsman Apr 11 at 12:15
maybe a good question for scicomp stackexchange –  meij Apr 11 at 18:32
@meij, was not aware of its existence, thanks! –  Johannes Apr 12 at 7:18
add comment

## 2 Answers

What you need is called "computing the action of the matrix exponential" (that is, computing $\exp(A)b$ without forming $\exp(A)$ explicitly. There are techniques based on complex integrals and Krylov subspaces. See http://eprints.ma.man.ac.uk/1426/ and the references included there.

-
This was exactly what I was looking for, thanks. The link even contains the authors' matlab code, which works very nicely. –  Johannes Apr 16 at 13:40
add comment

Note sure if it helps but,

    h = 0.00001;
I = eye(size(Q));
A = (I + h*Q);
p = A^(t/h)


also solves the problem.

-
Hm this is of course one way of approximating $e^{Qt}$, but it seems unlikely that it is better than the algorithm implemented in expm. –  Johannes Apr 12 at 7:42
add comment