find Markov steady state with left eigenvalues (using numpy or scipy)

Question

I need to find the steady state of Markov models using the left eigenvectors of their transition matrices using some python code.

It has already been established in this question that scipy.linalg.eig fails to provide actual left eigenvectors as described, but a fix is demonstrated there. The official documentation is mostly useless and incomprehensible as usual.

A bigger problem than than the incorrect format is that the eigenvalues produced are not in any particular order (not sorted and different each time). So if you want to find the left eigenvectors that correspond to the 1 eigenvalues you have to hunt for them, and this poses it's own problem (see below). The math is clear, but how to get python to compute this and return the correct eigenvectors is not clear. Other answers to this question, like this one, don't seem to be using the left eigenvectors, so those can't be correct solutions.

This question provides a partial solution, but it doesn't account for the unordered eigenvalues of larger transition matrices. So, just using

leftEigenvector = scipy.linalg.eig(A,left=True,right=False)[1][:,0]
leftEigenvector = leftEigenvector / sum(leftEigenvector)

is close, but doesn't work in general because the entry in the [:,0] position may not be the eigenvector for the correct eigenvalue (and in my case it is usually not).

Okay, but the output of scipy.linalg.eig(A,left=True,right=False) is an array in which the [0] element is a list of each eigenvalue (not in any order) and that is followed in position [1] by an array of eigenvectors in a corresponding order to those eigenvalues.

I don't know a good way to sort or search that whole thing by the eigenvalues to pull out the correct eigenvectors (all eigenvectors with eigenvalue 1 normalized by the sum of the vector entries.) My thinking is to get the indices of the eigenvalues that equal 1, and then pull those columns from the array of eigenvectors. My version of this is slow and cumbersome. First I have a function (that doesn't quite work) to find positions in a last that matches a value:

# Find the positions of the element a in theList
def findPositions(theList, a):
  return [i for i, x in enumerate(theList) if x == a]

Then I use it like this to get the eigenvectors matching the eigenvalues = 1.

M = transitionMatrix( G )
leftEigenvectors = scipy.linalg.eig(M,left=True,right=False)
unitEigenvaluePositions = findPositions(leftEigenvectors[0], 1.000)
steadyStateVectors = []
for i in unitEigenvaluePositions:
    thisEigenvector = leftEigenvectors[1][:,i]
    thisEigenvector / sum(thisEigenvector)
    steadyStateVectors.append(thisEigenvector)
print steadyStateVectors

But actually this doesn't work. There is one eigenvalue = 1.00000000e+00 +0.00000000e+00j that is not found even though two others are.

My expectation is that I am not the first person to use python to find stationary distributions of Markov models. Somebody who is more proficient/experienced probably has a working general solution (whether using numpy or scipy or not). Considering how popular Markov models are I expected there to be a library just for them to perform this task, and maybe it does exist but I couldn't find one.

Answer 1

You linked to How do I find out eigenvectors corresponding to a particular eigenvalue of a matrix? and said it doesn't compute the left eigenvector, but you can fix that by working with the transpose.

For example,

In [901]: import numpy as np

In [902]: import scipy.sparse.linalg as sla

In [903]: M = np.array([[0.5, 0.25, 0.25, 0], [0, 0.1, 0.9, 0], [0.2, 0.7, 0, 0.1], [0.2, 0.3, 0, 0.5]])

In [904]: M
Out[904]: 
array([[ 0.5 ,  0.25,  0.25,  0.  ],
       [ 0.  ,  0.1 ,  0.9 ,  0.  ],
       [ 0.2 ,  0.7 ,  0.  ,  0.1 ],
       [ 0.2 ,  0.3 ,  0.  ,  0.5 ]])

In [905]: eval, evec = sla.eigs(M.T, k=1, which='LM')

In [906]: eval
Out[906]: array([ 1.+0.j])

In [907]: evec
Out[907]: 
array([[-0.32168797+0.j],
       [-0.65529032+0.j],
       [-0.67018328+0.j],
       [-0.13403666+0.j]])

In [908]: np.dot(evec.T, M).T
Out[908]: 
array([[-0.32168797+0.j],
       [-0.65529032+0.j],
       [-0.67018328+0.j],
       [-0.13403666+0.j]])

To normalize the eigenvector (which you know should be real):

In [913]: u = (evec/evec.sum()).real

In [914]: u
Out[914]: 
array([[ 0.18060201],
       [ 0.36789298],
       [ 0.37625418],
       [ 0.07525084]])

In [915]: np.dot(u.T, M).T
Out[915]: 
array([[ 0.18060201],
       [ 0.36789298],
       [ 0.37625418],
       [ 0.07525084]])

If you don't know the multiplicity of eigenvalue 1 in advance, see @pv.'s comment showing code using scipy.linalg.eig. Here's an example:

In [984]: M
Out[984]: 
array([[ 0.9 ,  0.1 ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.3 ,  0.7 ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.25,  0.75,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.5 ,  0.5 ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  1.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  1.  ,  0.  ]])

In [985]: import scipy.linalg as la

In [986]: evals, lvecs = la.eig(M, right=False, left=True)

In [987]: tol = 1e-15

In [988]: mask = abs(evals - 1) < tol

In [989]: evals = evals[mask]

In [990]: evals
Out[990]: array([ 1.+0.j,  1.+0.j,  1.+0.j])

In [991]: lvecs = lvecs[:, mask]

In [992]: lvecs
Out[992]: 
array([[ 0.9486833 ,  0.        ,  0.        ],
       [ 0.31622777,  0.        ,  0.        ],
       [ 0.        , -0.5547002 ,  0.        ],
       [ 0.        , -0.83205029,  0.        ],
       [ 0.        ,  0.        ,  0.70710678],
       [ 0.        ,  0.        ,  0.70710678]])

In [993]: u = lvecs/lvecs.sum(axis=0, keepdims=True)

In [994]: u
Out[994]: 
array([[ 0.75, -0.  ,  0.  ],
       [ 0.25, -0.  ,  0.  ],
       [ 0.  ,  0.4 ,  0.  ],
       [ 0.  ,  0.6 ,  0.  ],
       [ 0.  , -0.  ,  0.5 ],
       [ 0.  , -0.  ,  0.5 ]])

In [995]: np.dot(u.T, M).T
Out[995]: 
array([[ 0.75,  0.  ,  0.  ],
       [ 0.25,  0.  ,  0.  ],
       [ 0.  ,  0.4 ,  0.  ],
       [ 0.  ,  0.6 ,  0.  ],
       [ 0.  ,  0.  ,  0.5 ],
       [ 0.  ,  0.  ,  0.5 ]])

Answer 2

All right, I had to make some changes while implementing Warren's solution and I've included those below. It's basically the same so he gets all the credit, but the realities of numerical approximations with numpy and scipy required more massaging which I thought would be helpful to see for others trying to do this in the future. I also changed the variable names to be super noob-friendly.

Please let me know if I got anything wrong or there are further recommended improvements (e.g. for speed).

# in this case my Markov model is a weighted directed graph, so convert that nx.graph (G) into it's transition matrix
M = transitionMatrix( G )   

#create a list of the left eigenvalues and a separate array of the left eigenvectors
theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True)  

# for stationary distribution the eigenvalues and vectors are always real, and this speeds it up a bit
theEigenvalues = theEigenvalues.real                 
leftEigenvectors = leftEigenvectors.real

# set how close to zero is acceptable as being zero...1e-15 was too low to find one of the actual eigenvalues
tolerance = 1e-10 

# create a filter to collect the eigenvalues that are near enough to zero                               
mask = abs(theEigenvalues - 1) < tolerance           

# apply that filter
theEigenvalues = theEigenvalues[mask]                

# filter out the eigenvectors with non-zero eigenvalues
leftEigenvectors = leftEigenvectors[:, mask] 

# convert all the tiny and negative values to zero to isolate the actual stationary distributions    
leftEigenvectors[leftEigenvectors < tolerance] = 0  

# normalize each distribution by the sum of the eigenvector columns
attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True)   

# this checks that the vectors are actually the left eigenvectors, but I guess it's not needed to usage 
#attractorDistributions = np.dot(attractorDistributions.T, M).T 

# convert the column vectors into row vectors (lists) for each attractor (the standard output for this kind of analysis)
attractorDistributions = attractorDistributions.T

# a list of the states in any attractor with the approximate stationary distribution within THAT attractor (e.g. for graph coloring)         
theSteadyStates = np.sum(attractorDistributions, axis=1)  

Putting that all together in an easy copy-and-paste format:

M = transitionMatrix( G ) 
theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True)  
theEigenvalues = theEigenvalues.real                 
leftEigenvectors = leftEigenvectors.real
tolerance = 1e-10            
mask = abs(theEigenvalues - 1) < tolerance 
theEigenvalues = theEigenvalues[mask]    
leftEigenvectors = leftEigenvectors[:, mask] 
leftEigenvectors[leftEigenvectors < tolerance] = 0  
attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True)   
attractorDistributions = attractorDistributions.T
theSteadyStates = np.sum(attractorDistributions, axis=0)  

Using this analysis on a generated Markov model produced one attractor (of three) with a steady state distribution of 0.19835218 and 0.80164782 compared to the mathematically accurate values of 0.2 and 0.8. So that's more than 0.1% off, kind of a big error for science. That's not a REAL problem because if accuracy is important then, now that the individual attractors have been identified, a more accurate analyses of behavior within each attractor using a matrix subset can be done.