33
$\begingroup$

EDIT: I am testing if any eigenvalues have a magnitude of one or greater.

I need to find the largest absolute eigenvalue of a large sparse, non-symmetric matrix.

I have been using R's eigen() function, which uses the QR algo from either EISPACK or LAPACK to find all eigenvalues and then I use abs() to get the absolute values. However, I need to do it faster.

I have also tried using the ARPACK interface in igraph R package. However, it gave an error for one of my matrices.

The final implementation must be accessible from R.

There will probably be multiple eigenvalues of the same magnitude.

Do you have any suggestions?

EDIT: Accuracy only needs to be to 1e-11. A "typical" matrix has so far been $386\times 386$. I have been able to do a QR factorisation on it. However, it is also possible to have much larger ones. I am currently starting to read about the Arnoldi algorithm. I understand that it is related to Lanczsos.

EDIT2: If I have multiple matrices that I am "testing" and I know that there is a large submatrix that does not vary. Is it possible to ignore/discard it?

Improve this question
$\endgroup$
6
  • $\begingroup$ See my answer here: scicomp.stackexchange.com/a/1679/979 . This is a current research topic and current methods can do better than Lanczos. The problem of computing singular values is equivalent to the problem of computing eigenvalues. $\endgroup$
    – dranxo
    Mar 23, 2012 at 1:47
  • 4
    $\begingroup$ 400x400 matrix != large. Also what does largest mean if "There will probably be multiple eigenvalues of the same magnitude."? In numpy land: linalg.eig(random.normal(size=(400,400))) takes about half a second. Is this too slow? $\endgroup$
    – meawoppl
    Feb 21, 2013 at 19:08
  • $\begingroup$ @meawoppl yes half a second is too slow. This is because it is part of another algo that runs this calculation many times. $\endgroup$
    – power
    Feb 26, 2013 at 3:29
  • 2
    $\begingroup$ @power gotcah. Do you have an approximation to the eigenvector. i.e. is it likely similar to the last solution, or can you make a educated guess about its structure? $\endgroup$
    – meawoppl
    Feb 27, 2013 at 22:30
  • 1
    $\begingroup$ Related: mattermodeling.stackexchange.com/q/6563/5 $\endgroup$
    – user1271772
    Aug 13, 2021 at 18:03

6 Answers 6

Reset to default
17
$\begingroup$

It depends a lot on the size of your matrix, in the large-scale case also on whether it is sparse, and on the accuracy you want to achieve.

If your matrix is too large to allow a single factorization, and you need high accuracy, the Lanczsos algorithm is probably the fastest way. In the nonsymmetric case, the Arnoldi algorithm is needed, which is numerically unstable, so an implementation needs to address this (is somewhat awkward to cure).

If this is not the case in your problem, give more specific information in your question. Then add a comment to this answer, and I'll update it.

Edit: [This was for the old version of the question, asling for the largest eigenvalue.] As your matrix is small and apparently dense, I'd do Arnoldi iteration on B=(I-A)^{-1}, using an initial permuted triangular factorization of I-A to have cheap multiplication by B. (Or compute an explicit inverse, but this costs 3 times as much as the factorization.) You want to test whether B has a negative eigenvalue. Working with B in place of A, negative eigenvalues are much better separated, so if there is one, you should converge rapidly.

But I am curious about where your problem comes from. Nonsymmetric matrices usually have complex eigenvalues, so ''largest'' isn't even well-defined. Thus you must know more about your problem, which might help in suggesting how to solve it even faster and/or more reliably.

Edit2: It is difficult to get with Arnoldi a particular subset of interest. To get the absolutely largest eigenvalues reliably, you'd do subspace iteration using the original matrix, with a subspace size matching or exceeding the number of eigenvalues expected to be close to 1 or larger in magnitude. On small matrices, this will be slower than the QR algorithm but on large matrices it will be much faster.

Improve this answer
$\endgroup$
11
  • $\begingroup$ I need to test if the largest eigenvalue is greater than 1. Accuracy only needs to be to 1e-11. A "typical" matrix has so far been 386 x 386. I have been able to do a QR factorisation on it. However, it is also possible to have much larger ones. I am currently starting to read about the Arnoldi algorithm. I understand that it is related to Lanczsos. $\endgroup$
    – power
    Mar 21, 2012 at 0:21
  • $\begingroup$ This information belongs to your question - so please edit it, and also add more information (why are the eigenvalues real? or what does largest mean?) - see the edit of my answer. $\endgroup$
    – Arnold Neumaier
    Mar 21, 2012 at 10:48
  • $\begingroup$ sorry that I did not explain myself clearly. I also did not explain clearly that the eigenvalues are complex. I am testing if any eigenvalues have a magnitude of one or greater. $\endgroup$
    – power
    Mar 23, 2012 at 1:02
  • 1
    $\begingroup$ This makes more sense, but now my recipe with $(I-A)^{-1}$ works well only if the poor eigenvalue is indeed real >1. On the other hand, the new info probably implies that you have little choice but computing all eigenvalues. - Please upfdate your question to convey the extra info! $\endgroup$
    – Arnold Neumaier
    Mar 23, 2012 at 11:06
  • 1
    $\begingroup$ see edit 2 in my answer $\endgroup$
    – Arnold Neumaier
    Mar 30, 2012 at 8:38
10
$\begingroup$

The Power Iteration (or Power Method), e.g. what Dan is describing, should always converge, albeit at the rate $\left|\lambda_{n-1}/\lambda_{n}\right|$.

If $\lambda_{n-1}$ is close to $\lambda_n$, it will be slow, but you can use extrapolation to get around that. It may seem complicated, but an implementation in pseudo-code is given in the paper.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ what if |λ(n−1)| = |λ(n)| ? $\endgroup$
    – power
    Mar 20, 2012 at 3:24
  • $\begingroup$ @power, then the regular Power Iteration won't converge. I don't know how well the extrapolation methods will distinguish between the different eigenvalues, you'll have to read the paper for that. $\endgroup$
    – Pedro
    Mar 20, 2012 at 10:21
  • 3
    $\begingroup$ @power: All things reconsidered, if $|\lambda_{n-1}|=|\lambda_n|$, then the power iteration will still converge to the correct eigenvalue. The resulting eigenvector, which you don't seem to be interested in anyway, will be a linear combination of the eigenvectors corresponding to $\lambda_n$ and $\lambda_{n-1}$. $\endgroup$
    – Pedro
    Mar 30, 2012 at 8:42
  • $\begingroup$ do you have a reference to an academic paper or book that supports this? Also, what if \lambda_{n} is complex? $\endgroup$
    – power
    Apr 2, 2012 at 2:34
  • 7
    $\begingroup$ If there are several different eigenvalues of maximal modulus, the power iteration converges only under exceptional circumstances. It generally oscillates in a somewhat unprdictable manner. $\endgroup$
    – Arnold Neumaier
    Apr 2, 2012 at 6:43
7
$\begingroup$

There has been some good research on this recently. The new approaches use "randomized algorithms" which only require a few reads of your matrix to get good accuracy on the largest eigenvalues. This is in contrast to power iterations which require several matrix-vector multiplications to reach high accuracy.

You can read more about the new research here:

http://math.berkeley.edu/~strain/273.F10/martinsson.tygert.rokhlin.randomized.decomposition.pdf

http://arxiv.org/abs/0909.4061

This code will do it for you:

http://cims.nyu.edu/~tygert/software.html

https://bitbucket.org/rcompton/pca_hgdp/raw/be45a1d9a7077b60219f7017af0130c7f43d7b52/pca.m

http://code.google.com/p/redsvd/

https://cwiki.apache.org/MAHOUT/stochastic-singular-value-decomposition.html

If your language of choice isn't in there you can roll your own randomized SVD pretty easily; it only requires a matrix vector multiplication followed by a call to an off-the-shelf SVD.

Improve this answer
$\endgroup$
1
5
$\begingroup$

Here you will find an algorithmic introduction to the Jacobi-Davidson (JD) algorithm, which computes the maximum eigenvalue (by modulus/absolute value).

In this paper the mathematical aspects are explored. JD allows general (real or complex) matrices and can be used to compute ranges of eigenvalues.

Here you can find various library implementations JDQR and JDQZ (including a C interface, which you should be able to link to from R).

Improve this answer
$\endgroup$
8
  • $\begingroup$ I have not been able to find any literature that explicitly states that the Jacobi-Davidson method works for a real, general matrix. $\endgroup$
    – power
    Mar 20, 2012 at 2:15
  • $\begingroup$ Unless every article explicitly states a restriction and the convergence argument relies on the restriction that doesn't matter. $\endgroup$
    – Deathbreath
    Mar 20, 2012 at 12:38
  • $\begingroup$ Here is another explanation of JD. The considered matrices are completely general. No special structure is exploited and results specific to Hermitian matrices are compared and contrasted, e.g., convergence for general matrices is quadratic, but cubic for Hermitian matrices. $\endgroup$
    – Deathbreath
    Mar 20, 2012 at 12:44
  • $\begingroup$ thanks for this. I not find any C code for a general matrix, so I will have to write my own. The links to the algorithms seem to be only for Hermetian matrices. $\endgroup$
    – power
    Mar 27, 2012 at 9:10
  • 1
    $\begingroup$ @power you will also not find in literature a result that states that the standard QR implementations converge for a real, general matrix -- that is an open problem, and indeed not long ago a counterexample was found for the QR code in LAPACK. $\endgroup$
    – Federico Poloni
    Oct 22, 2014 at 9:07
3
$\begingroup$

In your original post, you say:

"I have also tried using the ARPACK interface in igraph R package. However, it gave an error for one of my matrices."

I would be interested to know more about the error. If you can make this matrix publicly available somewhere, I'd be interested in trying ARPACK on it.

Based on what I've read above, I would expect ARPACK would do a very good job of extracting the largest (or a few of the largest) eigenvalues of a sparse matrix. To be more specific, I would expect Arnoldi methods to work well for this case and that, of course, is what ARPACK is based on.

The slow convergence of the power method when there are closely-spaced eigenvalues in the region of interest was mentioned above. Arnoldi improves this by iterating with several vectors instead of the one in power method.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I will see if I can find my work from back then. I worked on this one year ago. $\endgroup$
    – power
    Mar 8, 2013 at 2:46
2
$\begingroup$

It's not the fastest way, but a reasonably quick way is to just hit an (initially random) vector with the matrix repeatedly, and then normalize every few steps. Eventually it will converge to the largest eigenvector, and the gain in norm for a single step is the associated eigenvalue.

This works best when the largest eigenvalue is substantially larger than any other eigenvalue. If another eigenvalue is close in magnitude to the largest, this will take a while to converge, and it may be difficult to determine if it has converged.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Thanks Dan, however: In my matrices, some of the other eigenvalues will be have a similar (if not the same) magnitude as the largest one. Is your method similar to Power Iteration and Rayleigh Quotient Iteration? Batterson and Smillie (1990) write that for some non-symmmetric matrices, Rayleigh Quotient Iteration will not converge. Batterson, S., Smillie, J (1990) "Rayleigh Quotient Iteration for Nonsymmetric Matrices", Mathematics of Computation, vol 55, num 191, P 169 - 178 $\endgroup$
    – power
    Mar 19, 2012 at 6:59
  • $\begingroup$ If other eigenvalues have the same magnitude as the largest one... then aren't those values also "the largest one" too? $\endgroup$
    – ely
    Mar 19, 2012 at 18:13
  • $\begingroup$ @EMS: They would still be "largest eigenvalues" but the presence of more than one largest would still kill convergence. $\endgroup$
    – Dan
    Mar 19, 2012 at 18:29
  • $\begingroup$ I'm just wondering which eigenvalue you want it to converge to. Things like Rayleigh quotient/Power method are meant when there is a distinct largest eigenvalue. Your question asks to find the largest eigenvalue, but then it sounds like this isn't actually well defined for your problem. I'm just misled by the title of the post. $\endgroup$
    – ely
    Mar 19, 2012 at 18:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.