Efficient computation of Markov chain transition probability matrix

Question

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.

Answer 1

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.

Answer 2

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.