# Computing the largest eigenvalue of a very large sparse matrix

I am trying to compute the asymptotic growth-rate in a specific combinatorial problem depending on a parameter $w$, using the Transfer-Matrix method. This amounts to computing the largest eigenvalue of the corresponding matrix.

For small values of $w$, the corresponding matrix is small and I can use the so-called power method — start with some vector, and multiply it by the matrix over and over, and under certain conditions you'll get the eigenvector corresponding to the largest eigenvalue. However, for the values of $w$ I am interested in, the matrix becomes too large, and thus the vector becomes too large, say, $n > 10^7$ entries or so, and it can't be contained in the computer's memory anymore and I need extra programming tricks or a very powerful computer.

As for the matrix itself, I don't need to store it in memory — I can access it as a black box, i.e., given $i, j$ I can return $A_{ij}$ via a simple computation. Also, the matrix has only $0$ and $1$ entries, and I believe it to be sparse (i.e., only around $\log n$ of the entries are $1$'s, $n$ being the number of rows/columns). However, the matrix is not symmetric.

Is there some method more space-effective for computation of eigenvalues for a case like this?

• A duplicate of math.stackexchange.com/questions/4368/… where it has already had a reply. Sep 10, 2010 at 9:20
• This type of black box, where each entry of a huge matrix is determined by some computation, unfortunately seems pretty much useless for sparse matrix computations. More useful would be a list of the nonzero entries, or at least a nice small subset of the entries containing the nonzero ones as a subset. Sep 10, 2010 at 9:44
• Is it worth copy and pasting the answer in here or closing this question? Sep 10, 2010 at 9:52

You could use the Arnoldi Iteration algorithm. This algorithm only requires the matrix A for matrix-vector multiplication. I'm expecting that you will be able to black-box the function v→Av. What you generate is an upper Hessenberg matrix H whose eigenvalues whose can be computed cheaply (by a direct method or Rayleigh quotient iteration) and which approximate the eigenvalues of A. Arnoldi Iteration will give the best approximation to the dominant eigenvalue so I suspect you won't have to do many iterations before you have a good estimate.

An excellent introduction to this is: "Numerical Linear Algebra" by Trefethen and Bau. (p250)

The basic algorithm can be found here: http://en.wikipedia.org/wiki/Arnoldi_iteration

Now the only thing that is required to make this a fully functional algorithm is a termination condition. Since you don't seem to need the dominant eigenvalue to a high degree of accuracy I would not worry and just stop when the dominant eigenvalue estimate doesn't change too much.

If you have Matlab you can always use the built in function eigs(Afun,n,...) where Afun is the black-box function handle that computes Av.

• This is a copy of my answer from math.stackexchange.com/questions/4368/… Sep 10, 2010 at 10:00
• I think this combined with Ricky Demer's answer is a good answer to this question: first reduce the dimension of the space as much as possible (to about log(n)), then use an iterative method like Arnoldi. Sep 11, 2010 at 2:35
• I think you should put in a significant effort to use established packages -- in this case, probably ARPACK (caam.rice.edu/software/ARPACK). It sometimes feels like computer drudgery rather than research adventure, but it's likely to be a much quicker route to a reliable answer. I see in the Applications section of that website that someone got eigenvalues from a 2million-order matrix in 600 CPU hours. If nothing else, this should motivate you to use an efficient algorithm (and 600 CPU hours actually sounds really good for a problem this size). Sep 15, 2010 at 8:17

. nonzeroindices = set([])
. for i in range(n):
.. for j in range(n):
... if A(i,j) == 1:
.... nonzeroindices.union(set([i]))
. nonzeroindex = sorted(nonzeroindices)
. newmatrix = []
. for i in range(n):
.. newmatrix.append([])
.. for j in range(n):
... newmatrix[i].append(A(nonzeroindex[i],nonzeroindex[j]))

Then compute the eigenvalues of the around 2*log(n)-by-2*log(n) newmatrix.

• You don't even have to care about whether a column contains nonzero entries: if the row doesn't, then no eigenvector (except for the eigenvalue 0) can have a nonzero entry in that position. So you can just make a list of all the rows with nonzero entries and delete all the other ones (and the corresponding columns). Sep 10, 2010 at 9:55
• Assuming you have the transform operator as well, the sets of nonzeros can be computed in linear time (say by multiplying by random vectors, or doing the same thing symbolically). Jun 12, 2018 at 22:00