Solving large Markov Chains

Date:2008-12-28 (last modified), 2008-12-28 (created)

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:

In [ ]:
#!/usr/bin/env python

labda, mu1, mu2 = 1., 1.01, 1.001
N1, N2 = 50, 50
size = N1*N2

from numpy import ones, zeros, empty
import scipy.sparse as sp
import  pysparse
from pylab import matshow, savefig
from scipy.linalg import norm
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.

In [ ]:
def state(i,j):
    return j*N1 + i

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

In [ ]:
def fillOffDiagonal(Q):
    # labda
    for i in range(0,N1-1):
        for j in range(0,N2):
            Q[(state(i,j),state(i+1,j))]= labda
    # mu2
    for i in range(0,N1):
        for j in range(1,N2):
            Q[(state(i,j),state(i,j-1))]= mu2
    # mu1
    for i in range(1,N1):
        for j in range(0,N2-1):
            Q[(state(i,j),state(i-1,j+1))]= mu1
    print "ready filling"

In this function we use scipy.sparse

In [ ]:
def computePiMethod1():
    e0 = time.time()
    Q = sp.dok_matrix((size,size))
    # Set the diagonal of Q such that the row sums are zero
    Q.setdiag( -Q*ones(size) )
    # Compute a suitable stochastic matrix by means of uniformization
    l = min(Q.values())*1.001  # avoid periodicity, see the book of Bolch et al.
    P = sp.speye(size, size) - Q/l
    # compute Pi
    P =  P.tocsr()
    pi = zeros(size);  pi1 = zeros(size)
    pi[0] = 1;
    n = norm(pi - pi1,1); i = 0;
    while n > 1e-3 and i < 1e5:
        pi1 = pi*P
        pi = pi1*P   # avoid copying pi1 to pi
        n = norm(pi - pi1,1); i += 1
    print "Method 1: ", time.time() - e0, i
    return pi

Now use Pysparse.

In [ ]:
def computePiMethod2():
    e0 = time.time()
    Q = pysparse.spmatrix.ll_mat(size,size)
    # fill diagonal
    x =  empty(size)
    # uniformize
    l = min(Q.values())*1.001
    P = pysparse.spmatrix.ll_mat(size,size)
    P.shift(-1./l, Q)
    # Compute pi
    P = P.to_csr()
    pi = zeros(size);  pi1 = zeros(size)
    pi[0] = 1;
    n = norm(pi - pi1,1); i = 0;
    while n > 1e-3 and i < 1e5:
        n = norm(pi - pi1,1); i += 1
    print "Method 2: ", time.time() - e0, i
    return pi

Output the results.

In [ ]:
def plotPi(pi):
    pi = pi.reshape(N2,N1)
if __name__ == "__main__":
    pi = computePiMethod1()
    pi = computePiMethod2()

Here is the result:

Improvements of this Tutorial

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

Section author: nokfi