# How can I efficiently approximate the stationary distribution of an infinite CTMC with a sparse rate matrix?

I am looking for methods to approximate the stationary distribution of an infinite CTMC with a sparse rate matrix. Each row and column of the rate matrix has a finite number of non-zero elements. There is this paper Bounding the equilibrium distribution of Markov population models by Dayar et al. that identifies a region of the state space where most of the probability of the stationary distribution is concentrated and then finds upper and lower bounds on the stationary probabilities for the states in that region.

I wondered whether there was something similar for other classes of CTMC or something a bit more general.

A lot of what I've found has been applicable to various classes of queue systems but my case doesn't fall within that.

• you mention that "my case doesn't fall within that", but you do not really tell us what your case is; I guess you'll get a better response if you can be more specific. – Carlo Beenakker May 31 '16 at 14:05

Since each row $r$ has a finite number of non-zero elements, say $N_r \subset \{1, 2, \dots\}$, you could use the Gillespie algorithm to simulate the continuous time Markov chain. A stationary distribution can then be approximated by the normalized occupancy times. If the matrix is irreducible, this will be the stationary distribution.

Here is an outline of how this would work:

Gillespie algorithm

1. Given current state $X_t = r$, simulate an exponential random variable of rate $\lambda_r = \sum_{c \in N_r, c\neq r} q_{r,c}$. This corresponds to the time $t' = t + \lambda_r$ where the current state jumps to a different one.
2. Pick a new distinct state $X'_{t'}$ from a categorical random variable with probability $P(X_{t'} = c|X_{t'-} = r) = q_{r,c}/\lambda_r$ for $c \neq r$.
3. Repeat until a trajectory of fixed length $T$ has been produced.

Occupancy time

The approximation of a stationary distribution is then given by

$$\pi(x) \approx \frac{\text{length}\{t \in [0, T) : X_t = x\}}{T}.$$

• This Monte-Carlo technique is also known as the Doob-Gillespie or Stochastic Simulation Algorithm. If the process $\{ X_t \}$ has a stationary distribution and satisfies a strong LLN then: $\pi(x) = \lim_{t \to \infty} \frac{1}{t} \int_0^t \mathbf{1}_{\{X_s = x\}} ds$. This quantity basically gives the fraction of time spent by the (infinitely long) trajectory $\{ X_t \}$ in state $x$. – Nawaf Bou-Rabee Oct 2 '16 at 13:51

(I assume that you mean "approximate" as in "compute an approximate solution with a numerical method", but I am not sure if I am correct)

The "standard approach" for someone working in numerical linear algebra would be this one: let's assume that the Markov chain is irreducible and that $Q$ is the rate matrix; let us choose our notation so that $\mathbf{1}^T Q = \mathbf{0}^T.$ then, we can partition the system $Qv=0$ as $$\begin{bmatrix} q_{11} & q_{12}\\ q_{21} & Q_{22} \end{bmatrix} \begin{bmatrix} v_1 \\ v_{2} \end{bmatrix} = \begin{bmatrix} 0\\0 \end{bmatrix},$$ where $q_{11}$ is scalar and $Q_{22}$ is $(n-1)\times(n-1)$. We can assume, up to changing normalization, that $v_1=1$. Then, $v_2$ satisfies $-Q_{22}v_2 = q_{21}$, which is a linear system with a nonsingular sparse M-matrix as its matrix. Solution methods for sparse linear systems are well studied; one of the most famous algorithms is GMRES.

• Please note that the rate matrix $Q$ is an infinite matrix, so standard numerical linear algebra tools simply won't work unless $Q$ is truncated. – Nawaf Bou-Rabee Oct 2 '16 at 15:23
• @NawafBou-Rabee Thanks, I had missed the "infinite" part. – Federico Poloni Oct 2 '16 at 15:33