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