CODEWITHSUNDEEP
pythonnumpylinear-algebraeigenvectormarkov-chains
mepstein1218
mepstein1218

Reputation: 331

Python numpy/scipy eigenvectors seemingly not correct for markov chain model

I have a large (351,351) numpy transition matrix. I would like to find the state state vector for this using numpy (I also tried scipy which has the same exact function).

sstate = np.linalg.eig(T)[1][:,0]

So this I believe should give me the eigenvector for the dominant left eigenvalue. The dominant left eigenvalue is 1+0j. This is somewhat correct, the dominant left eigenvalue should be 1, I'm a little new to this so I am not sure what to do with the imaginary numbers. Also, the sstate vector contains all complex numbers. Now, trying to check if this is correct, I do the following matrix multiplication:

np.dot(sstate,T)

If done correctly, this should return the same vector as 'sstate.' I'm not sure why this is not working. Could the imaginary numbers be the problem? Also, is it possible this transition matrix does not contain a steady state vector. Each row and column in my transition state matrix is supposed to sum to 1, however, I've found that rounding errors cause each row and column sum to only be approximately 1.

Any and all help is appreciated!

Upvotes: 4

Views: 2533

Answers (2)

ely
ely

Reputation: 77404

Is the transition matrix symmetric? If not, consider checking for T.T (the transpose), because you need to make sure you're looking at the right state transitions: you need the left eigenvector of your stochastic matrix, but almost all out-of-the-box scientific packages (numpy included) default to computing the right eigenvectors (this the same reason in textbooks and stuff you have to premultiply by row vectors instead of usual matrix-column multiplication when dealing with this stuff).

Maybe also sstate = sstate/sstate.sum() to make sure the probabilities sum to 1 despite roundoff.

Here's an example with numpy.

Added detail on right vs. left eigenvectors from comments:

eig and things like it will compute the right eigenvectors, as in vectors v such that Av = (lambda)v for scalar lambda. What you need though is the left eigenvector of A, so something that satisfies v.T*A = (lambda)v.T and this won't just be a transpose or conjugate of the right eigenvector.

So you will want to compute the eigenvector based on A.T but you will not want to compute with A.T later when checking whether the state vector really is stationary. You'll want to look at np.dot(sstate, T) (verify that sstate is a row vector, not column), and evaluate that (and possibly also the other point about renormalizing to help round off).

Upvotes: 3

doug
doug

Reputation: 70028

so look at the 1D array (the first element in the 2-tuple) returned from the call to eig; this is the eigenvalue array, and as you can see it is not in descending order, you need to sort it manually and then apply the same ordering to the eigenvector array. You might have already done that but it's not in your code snippet nor mentioned in the OP.

Once you do that, then you can select the first eigenvector:

>>> import numpy as NP
>>> from scipy import linalg as LA
>>> a = NP.random.rand(16).reshape(4, 4)

>>> E = LA.eig(a, left=True)

>>> evals, evecs = E

sort the eigenvalues in decreasing order

>>> idx = NP.argsort(evals)[::-1]

>>> eva = eva[idx]

apply the sort index to the eigenvector matrix

>>> eva[idx,]

select the first one

>>> eva[0].real

also, use the linalg from scipy; the NumPy installer includes a generic BLAS to build against if it's not on the machine

in addition, if the 2D array you are passing to eig is sparse, then use eig from scipy.sparse; it's much faster, particularly so in your case because you only want one eigenvector.

Upvotes: 0

Related Questions