11

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?

eigenchris
  • 5,791
  • 2
  • 21
  • 30
5xum
  • 5,250
  • 8
  • 36
  • 56
  • 2
    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: http://www.math.ualberta.ca/~ewoolgar/labs/linalg/Lab15.pdf. It certainly is much more efficient than doing a full out SVD. – rayryeng Mar 26 '15 at 23:06
  • 1
    @rayryeng 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? – 5xum Mar 26 '15 at 23:19
  • 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). – rayryeng Mar 26 '15 at 23:27
  • I've pinged Amro and chappjc to bring attention to your post. Hope they'll be able to solve your problem! – rayryeng Mar 26 '15 at 23:31
  • 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? – knedlsepp Mar 27 '15 at 10:08
  • @knedlsepp 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. – 5xum Mar 27 '15 at 12:15

1 Answers1

8

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(AI)=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 x1 = 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 AT (such that ATx = 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.

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • 2
    Nice approach. I also thought about setting a random vector `b`, then solving the equation `(A-I)(y) = (A-I)b`and 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`. – 5xum Mar 27 '15 at 06:35