Question: 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!

Answer: 
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.](https://math.stackexchange.com/questions/121736/what-is-the-expected-number-of-dice-one-needs-to-roll-to-get-any-monotonically-i/121748#121748)

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