There are quite a few algorithms available to numerically compute eigenvalues of a matrix. Suppose one is not interested in computing the full spectrum of a general matrix, but only the largest eigenvalue in magnitude, and that we have additional knowledge on the matrix at hand namely:

  • it has real and non-positive entries
  • it is symmetric
  • it is definitely not sparse
  • the range of sizes that should be handled is in between 10000 ~ 50000

At first sight, power iteration methods should be relevant for me as they allow for computation of the eigenvalues one by one, in order of decreasing magnitude, so I could just stop after the first one. But it does not really use my assumptions so I figured that it might be possible to do better than that. As it turns out, the Wikipedia entry states the following :

For symmetric matrices, the power iteration method is rarely used, since its convergence speed can be easily increased without sacrificing the small cost per iteration; see, e.g., Lanczos iteration and LOBPCG.

I can also see that the scipy Python package implements the Implicitly Restarted Lanczos Method, for its scipy.sparse.linalg.eigsh function that "Find k eigenvalues and eigenvectors of the real symmetric square matrix or complex hermitian matrix A.".

So it is now my understanding that the power method is not the way to go, but is the Lanczos iteration method the best direction for my problem ? Or is there a way to use my second assumption that the matrix only contains real entries of the same sign to speed up the process ?

Some preliminary experimental results for a small (2000 rows and cols) matrix :

naive full spectrum computation with scipy.linalg.eig
--- 0.7484889030456543 seconds ---
rho: 36.7184965271

full spectrum computation with scipy.sparse.linalg.eigsh:
--- 0.7304418087005615 seconds ---
rho: 36.7184965271

eigsh stopped at k=1:
--- 0.06366991996765137 seconds ---
rho: 36.7184965271

cvxpy for computing spectral norm based on Rodrigo de Azevedo's hint:
--- 0.19952964782714844 seconds ---
rho: 36.718496527064126

Interestingly moving from scipy.linalg.eig to scipy.sparse.linalg.eigsh does not speed up things a lot when computing a full spectrum. I also compared eigsh to the alternative method proposed by Rodrigo de Azevedo.

  • $\begingroup$ If your matrix is symmetric, the spectral radius is equal to the spectral norm, which can be computed using semidefinite programming. Take a look at this toy example. $\endgroup$ – Rodrigo de Azevedo Sep 6 '18 at 19:11
  • $\begingroup$ I completely agree with what you say, but how would you argue that it would be better to compute a spectral norm compared to a spectral radius ? $\endgroup$ – ippiki-ookami Sep 9 '18 at 8:56
  • $\begingroup$ The spectral norm can be computed via convex optimization. Can one say the same about the spectral radius? $\endgroup$ – Rodrigo de Azevedo Sep 9 '18 at 9:02
  • $\begingroup$ Perhaps it is not, but I tried your idea with a randomly generated ~2000 rows toy matrix, and it seems that on my machine the Lanczos method stopped at the first eigenvalue outperforms computing the spectral norm through convex optimization. Is this a surprising result to you ? Should I prefer the latter for other reasons (amenability to parallel computation, etc..) ? $\endgroup$ – ippiki-ookami Sep 9 '18 at 9:23
  • $\begingroup$ Frankly, I would just compute the dominant eigenvalue of $A^\top A = A^2$. My first comment was not very smart; using SDP is interesting when one wants to find a matrix $A$ that has minimal spectral norm. You don't have a design problem, however. Your $A$ is given. $\endgroup$ – Rodrigo de Azevedo Sep 9 '18 at 9:41

When you call the Python routine it actually calls Fortran routines I believe that choose specific sub-routines. If you look at Netlib. It has many routines for symmetric eigenvalues problems which are all Lanczos algorithms. I am not sure about the only positive part. It will take into account the fact your matrix isn't sparse.

However, the issue is with the size of your matrix. A dense 10k x 10k matrix may actually work. A dense 50k x 50k matrix is more than 12 GB.

In terms of the fastest, there are likely libraries if you look them up for low-rank matrix approximation based on randomized algorithms.


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.