4

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!

mepstein1218
  • 331
  • 5
  • 14

2 Answers2

3

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).

Community
  • 1
  • 1
ely
  • 74,674
  • 34
  • 147
  • 228
  • I just checked, the matrix is not symmetric, which is odd because it should be. I will have to check this out. – mepstein1218 Feb 13 '15 at 02:43
  • Transitions don't *have* to be symmetric, but they often are. The one in the example I linked to, for instance, is not symmetric because the problem states are not symmetric. – ely Feb 13 '15 at 02:45
  • Oh OK, after looking over the math, my matrix should not be symmetric. So how does this cause a problem with numpy/scipy? – mepstein1218 Feb 13 '15 at 02:48
  • I tried finding the steady state vector for the transposed transition matrix, T.T, but even still this steady state vector multiplied by T.T does not equal itself. – mepstein1218 Feb 13 '15 at 02:51
  • `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, – ely Feb 13 '15 at 02:51
  • Again, 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). – ely Feb 13 '15 at 02:54
0

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.

doug
  • 69,080
  • 24
  • 165
  • 199