title: Markov chain stationary distributions with scipy.sparse?
body: <p>I have a Markov chain given as a large sparse <code>scipy</code> matrix <code>A</code>. (I've constructed the matrix in <code>scipy.sparse.dok_matrix</code> format, but converting to other ones or constructing it as <code>csc_matrix</code> are fine.)</p>
  
  <p>I'd like to know any stationary distribution <code>p</code> of this matrix, which is an eigenvector to the eigenvalue <code>1</code>. All entries in this eigenvector should be positive and add up to 1, in order to represent a probability distribution.</p>
  
  <p>This means I want any solution for the system
  <code>(A-I) p = 0</code>, <code>p.sum()=1</code> (where <code>I=scipy.sparse.eye(*A.shape)</code> is the idententy matrix), but <code>(A-I)</code> will not be of full rank, and even the whole system may be under-determined. In addition, eigenvectors with negative entries can be generated, which cannot be normalized to be valid probability distributions. Preventing negative entries in <code>p</code> would be nice.</p>
  
  <ul>
  <li><p>Using <code>scipy.sparse.linalg.eigen.eigs</code> is not a solution:
  It does not permit specifying the additive constraint. (Normalizing does not help if the eigenvectors contain negative entries.) Furthermore, it deviates quite a bit from the true results, and sometimes has problems converging, behaving worse than <code>scipy.linalg.eig</code>. (Also, I used shift-invert mode, which improves finding the types of eigenvalues I want, but not their quality. If I don't use it, it's even more overkill, because I am only interested in one particular eigenvalue, <code>1</code>.)</p></li>
  <li><p>Converting to dense matrices and using <code>scipy.linalg.eig</code> is not a solution: In addition to the negative entry problem, the matrix is too large.</p></li>
  <li><p>Using <code>scipy.sparse.spsolve</code> is not an obvious solution:
  The matrix is either not square (when combining the additive constraint and the eigenvector condition) or not of full rank (when trying to specify them separately in some way), sometimes neither.</p></li>
  </ul>
  
  <p>Is there a good way to numerically get a stationary state of a Markov chain given as sparse matrix using python?
  If there is a way to get an exhaustive list (and possibly also nearly stationary states), that's appreciated, but not necessary. </p>