This page shows how to compute the stationary distribution pi of a large Markov chain. The example is a tandem of two M/M/1 queues. Generally the transition matrix P of the Markov chain is sparse, so that we can either use scipy.sparse or Pysparse. Here we demonstrate how to use both of these tools.

## Power Method

In this section we find pi by means of iterative methods called the Power method. More specifically, given a (stochastic) transition matrix P, and an initial vector pi_0, compute iteratively pi_n = pi_{n-1} P until the distance (in some norm) between pi_n and pi_{n-1} is small enough.

Fist we build the generator matrix Q for the related Markov chain. Then we turn Q into a transition matrix P by the method of uniformization, that is, we define P as I - Q/l, where I is the identity matrix (of the same size as Q) and l is the smallest element on the diagonal of Q. Once we have P, we approximate pi (the left eigenvector of P that satisfies pi = pi P) by the iterates pi_n = pi_0 P^n, for some initial vector pi_0.

More details of the above approach can be found in (more or less) any book on probability and Markov Chains. A fine example, with many nice examples and attention to the numerical solution of Markov chains, is `Queueing networks and Markov Chains' by G. Bolch et al., John Wiley, 2nd edition, 2006.

You can get the source code for this tutorial here: tandemqueue.py

```
1 #!/usr/bin/env python
2
3 labda, mu1, mu2 = 1., 1.01, 1.001
4 N1, N2 = 50, 50
5 size = N1*N2
6
7 from numpy import ones, zeros, empty
8 import scipy.sparse as sp
9 import pysparse
10 from pylab import matshow, savefig
11 from scipy.linalg import norm
12 import time
```

This simple function converts the state (i,j), which represents that the first queue contains i jobs and the second queue j jobs, to a more suitable form to define a transition matrix.

```
1 def state(i,j):
2 return j*N1 + i
```

Build the off-diagonal elements of the generator matrix Q.

```
1 def fillOffDiagonal(Q):
2 # labda
3 for i in range(0,N1-1):
4 for j in range(0,N2):
5 Q[(state(i,j),state(i+1,j))]= labda
6 # mu2
7 for i in range(0,N1):
8 for j in range(1,N2):
9 Q[(state(i,j),state(i,j-1))]= mu2
10 # mu1
11 for i in range(1,N1):
12 for j in range(0,N2-1):
13 Q[(state(i,j),state(i-1,j+1))]= mu1
14 print "ready filling"
```

In this function we use scipy.sparse

```
1 def computePiMethod1():
2 e0 = time.time()
3 Q = sp.dok_matrix((size,size))
4 fillOffDiagonal(Q)
5 # Set the diagonal of Q such that the row sums are zero
6 Q.setdiag( -Q*ones(size) )
7 # Compute a suitable stochastic matrix by means of uniformization
8 l = min(Q.values())*1.001 # avoid periodicity, see the book of Bolch et al.
9 P = sp.speye(size, size) - Q/l
10 # compute Pi
11 P = P.tocsr()
12 pi = zeros(size); pi1 = zeros(size)
13 pi[0] = 1;
14 n = norm(pi - pi1,1); i = 0;
15 while n > 1e-3 and i < 1e5:
16 pi1 = pi*P
17 pi = pi1*P # avoid copying pi1 to pi
18 n = norm(pi - pi1,1); i += 1
19 print "Method 1: ", time.time() - e0, i
20 return pi
```

Now use Pysparse.

```
1 def computePiMethod2():
2 e0 = time.time()
3 Q = pysparse.spmatrix.ll_mat(size,size)
4 fillOffDiagonal(Q)
5 # fill diagonal
6 x = empty(size)
7 Q.matvec(ones(size),x)
8 Q.put(-x)
9 # uniformize
10 l = min(Q.values())*1.001
11 P = pysparse.spmatrix.ll_mat(size,size)
12 P.put(ones(size))
13 P.shift(-1./l, Q)
14 # Compute pi
15 P = P.to_csr()
16 pi = zeros(size); pi1 = zeros(size)
17 pi[0] = 1;
18 n = norm(pi - pi1,1); i = 0;
19 while n > 1e-3 and i < 1e5:
20 P.matvec_transp(pi,pi1)
21 P.matvec_transp(pi1,pi)
22 n = norm(pi - pi1,1); i += 1
23 print "Method 2: ", time.time() - e0, i
24 return pi
```

Output the results.

```
1 def plotPi(pi):
2 pi = pi.reshape(N2,N1)
3 matshow(pi)
4 savefig("pi.eps")
5 if __name__ == "__main__":
6 pi = computePiMethod1()
7 pi = computePiMethod2()
8 plotPi(pi)
```

Here is the result:

## Improvements of this Tutorial

Include other methods such as Jacobi's method or Gauss Seidel.