Quickly and efficiently calculating an eigenvector for known eigenvalue

Question

Short version of my question:

What would be the optimal way of calculating an eigenvector for a matrix A, if we already know the eigenvalue belonging to the eigenvector?

Longer explanation:

I have a large stochastic matrix A which, because it is stochastic, has a non-negative left eigenvector x (such that A^Tx=x).

I'm looking for quick and efficient methods of numerically calculating this vector. (Preferrably in MATLAB or numpy/scipy - since both of these wrap around ARPACK/LAPACK, any one would be fine).

I know that 1 is the largest eigenvalue of A, so I know that calling something like this Python code:

from scipy.sparse.linalg import eigs
vals, vecs = eigs(A, k=1)

will result in vals = 1 and vecs equalling the vector I need.

However, the thing that bothers me here is that calculating eigenvalues is, in general, a more difficult operation than solving a linear system, and, in general, if a matrix M has eigenvalue l, then finding the appropriate eigenvector is a matter of solving the equation (M - 1 * I) * x = 0, which is, in theory at least, an operation that is simpler than calculating an eigenvalue, since we are only solving a linear system, more specifically, finding the nullspace of a matrix.

However, I find that all methods of nullspace calculation in MATLAB rely on svd calculation, a process I cannot afford to perform on a matrix of my size. I also cannot call solvers on the linear equation, because they all only find one solution, and that solution is 0 (which, yes, is a solution, but not the one I need).

Is there any way to avoid calls to eigs-like function to solve my problem more quickly than by calculating the largest eigenvalue and accompanying eigenvector?

Osmium

Are you concerned with only finding the largest eigenvector or eigenvalue? If that's the case, consider using the Power Method. It is an iterative algorithm that simultaneously determines the most prominent eigenvector and eigenvalue: math.ualberta.ca/~ewoolgar/labs/linalg/Lab15.pdf. It certainly is much more efficient than doing a full out SVD. – Trestle

@Trestle Actually, MATLAB uses ARPACK to find the largest eigenpair even more quickly than a regular Power Method. But the main focus of my question is this: Since I already know what the largest eigenvalue is, does there exist a method to find the corresponding eigenvector that is better than methods that simply find the largest eigenpair? – Osmium

That unfortunately I don't have the answer to. I'm also curious to know what the answer is! Hope someone is able to answer your question (Amro, chappjc, horchler). – Trestle

I've pinged Amro and chappjc to bring attention to your post. Hope they'll be able to solve your problem! – Trestle

To dig a little deeper: What is your matrix size? Does it have any special structure? (sparse, triangle, symmetric, positive definite,...) Do you want a single eigenvector corresponding to your eigenvalue or all of them (in case there are others)? Do you need to compute these for many matrices? – Dhow

@Dhow The matrix is, well, let's just say large. Anything from several 10k to one million rows (and columns, since it is square). It is a sparse stochastic matrix, not symmetric. I know the eigenspace for eigenvalue 1 has dimension 1. – Osmium

Answer 1

Here's one approach using Matlab:

1. Let x denote the (row) left^† eigenvector associated to eigenvalue 1. It satisfies the system of linear equations (or matrix equation) xA = x, or x(A−I)=0.
2. To avoid the all-zeros solution to that system of equations, remove the first equation and arbitrarily set the first entry of x to 1 in the remaining equations.
3. Solve those remaining equations (with x₁ = 1) to obtain the other entries of x.

Example using Matlab:

>> A = [.6 .1 .3
        .2 .7 .1
        .5 .1 .4]; %// example stochastic matrix
>> x = [1, -A(1, 2:end)/(A(2:end, 2:end)-eye(size(A,1)-1))]
x =
   1.000000000000000   0.529411764705882   0.588235294117647
>> x*A %// check
ans =
   1.000000000000000   0.529411764705882   0.588235294117647

Note that the code -A(1, 2:end)/(A(2:end, 2:end)-eye(size(A,1)-1)) is step 3.

In your formulation you define x to be a (column) right eigenvector of A^T (such that A^Tx = x). This is just x.' from the above code:

>> x = x.'
x =
   1.000000000000000
   0.529411764705882
   0.588235294117647
>> A.'*x %// check
ans =
   1.000000000000000
   0.529411764705882
   0.588235294117647

You can of course normalize the eigenvector to sum 1:

>> x = x/sum(x)
x =
   0.472222222222222
   0.250000000000000
   0.277777777777778
>> A.'*x %'// check
ans =
   0.472222222222222
   0.250000000000000
   0.277777777777778

---

^{† Following the usual convention. Equivalently, this corresponds to a right eigenvector of the transposed matrix.}

Racemose

Nice approach. I also thought about setting a random vector b, then solving the equation (A-I)(y) = (A-I)band then setting x=y-b, since (A-I)b = (A-I)(x+b) = (A-I)x + (A-I)b implies that (A-I)(x)=0. Your case is an example of this for a specific b. – Osmium