# 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.

• Is $Q$ tridiagonal? That is, are the "3 diagonals only" the main diagonal and the ones immediately above and below that? – hardmath Apr 13 '13 at 8:19
• No, it's the main diagonal and two other diagonals (not the ones immediately above and below the main). – svangen Apr 15 '13 at 10:34

Computing an approximation of $P_1=e^Q$ is typically done by some kind of polynomial expansion of the matrix exponential. One doesn't take the Taylor expansion (it converges too slowly) but for the sake of the argument, let me assume that you do and that you approximate $$P_1 = e^Q \approx \sum_{k=0}^N \frac1{k!} Q^k = I+Q+\frac 1{2!} Q^2 + \frac 1{3!}Q^3 + \cdots = I+Q(I+\frac 12 Q(I+\frac 13 Q(\cdots)))$$ Then remember that you get the $j$th row of $P_1$ by forming the product $e_j^T P_1$ where $e_j$ is the $j$th unit vector $(0, \ldots 0, 1, 0, \ldots, 0)^T$ with the one in $j$th place of the vector. You can use this in the formula above to compute $$e_j^T P_1 \approx = e_j^T+(e_j^T Q)(I+\frac 12 Q(I+\frac 13 Q(\cdots)))$$ If you multiply this product out from left to right, you will find that all you have to ever do is form products of a vector (transposed) with a matrix. In other words, you never actually need any matrix-matrix products which would be expensive in your case because your matrix $Q$ is sparse but the powers of $Q$ will become less and less sparse. In the end, using this approach, all you need is $O(N)$ vector-matrix products to take the expansion to $N$th order.

As I said, in practice, one would not take a Taylor expansion but other polynomial expansions of the matrix exponential. However, the general approach will remain the same: reduce it all to a sequence of vector-matrix products.

• I second the suggestion not to compute the matrix exponential explicitly, but I cannot resist the opportunity to point out Moler and van Loan's Nineteen dubious ways to compute the matrix exponential here. Method 2, 3 and 4 are based on rational approximations (which are better suited for exponentials than polynomials) and could be used for this sort of approach as well. – Christian Clason Apr 12 '13 at 19:27
• Yes, this is clearly the paper to look at for this. I meant to include rational approximations -- let's just say they are "generalized polynomial expansions" ;-) Of course, in the current context, the matrix has only 3 entries per row which makes matrix-vector products particularly cheap, whereas rational approximations will suffer from having to do linear solves which are much more expensive. – Wolfgang Bangerth Apr 12 '13 at 20:08
• True, that's why I hesitated to bring it up (but it's just too nice a paper ;)). One would have to test whether the improved approximation properties outweigh the more expensive operations involved. It very much depends on the structure of the transition matrix. (From the question, it seems to be banded?) – Christian Clason Apr 12 '13 at 20:32

A software package mentioned in Moler and van Loan's "twenty-five years after" portion of their paper (Sec. 13, see link in Christian Clason's comment) is Roger Sidje's Expokit, with an accompanying paper that explains the approach.

For large sparse matrices a Taylor-like polynomial is applied to Krylov subspaces, obtaining the advantages of matrix-vector operations with sharper error bounds than a pure Taylor series approximation.