I want to compute the spectrum (all the eigenvalues) of a large sparse matrix (hundreds of thousands of rows). This is hard.

I am willing to settle for an approximation. Are there approximation methods to do this?

While I hope for a general answer to this question I would also be satisfied with an answer in the following specific case. My matrix is a Normalized Laplacian of a large graph. Eigenvalues will be between 0 and 2 with a large number of them clustered around 1.

  • $\begingroup$ Is the matrix sparse or dense? $\endgroup$ – Aron Ahmadia Jun 20 '12 at 7:31
  • $\begingroup$ The matrix is sparse. I've edited the question to reflect this. $\endgroup$ – MRocklin Jun 20 '12 at 13:11
  • $\begingroup$ Why do you want all the eigenvalues? This is almost universally a bad thing to do when you have a sparse or structured matrix, thus it's important to know how you plan to use it. $\endgroup$ – Jed Brown Jun 22 '12 at 5:49
  • $\begingroup$ The spectrum of a graph laplacian carries some important information that I'd like to inspect. I don't need them all, I just need to know roughly where they are. $\endgroup$ – MRocklin Jun 22 '12 at 14:38

If your graph is undirected (as I suspect), the matrix is symmetric, and you cannot do anything better than the Lanczsos algorithm (with selective reorthogonalization if necessary for stability). As the full spectrum consists of 100000 numbers, I giess you are mainly interested in the spectral density.

To get an approximate spectral density, take the spectrum of the leading Krylov subspace of dimension 100 or so, and replace its discrete density by a smoothed version.

The leading Krylov spectrum will have nearly resolved well-isolated eigenvalues (should these exist), approximates the eigenvalues at the end of the nonisolates spectrum, and is somewhat random in-between, with a distribution whose cumulative distribution function resembles that of the true spectrum. It would converge to it in exact arithmetic if the dimension grows. (If your operator were infinite-dimensional, this would still be the case, and you'd get the integral of the true spectral density function on the continuous spectrum.)

  • $\begingroup$ Won't the spectrum of the leading Krylov subspace just be the 100 largest eigenvalues? I'm also interested in the distribution of the moderate and smallest eigenvalues. $\endgroup$ – MRocklin Jun 20 '12 at 13:15
  • 1
    $\begingroup$ @MRocklin: No. I augmented my answer to give more details. $\endgroup$ – Arnold Neumaier Jun 21 '12 at 8:20

Arnold Neumaier's answer is discussed in more detail in section 3.2 of the paper "Approximating Spectral Densities of Large Matrices" by Lin Lin, Yousef Saad and Chao Yang (2016).

Some other methods are also discussed but the numerical analysis at the end of the paper shows that the Lanczos method outperforms these alternatives.


If you're ok with thinking about things that are not eigenvalues but functions that in some sense still tell you something about the spectrum, then I think you should check out some of the work by Mark Embree at Rice University.


Here's yet another way to characterize the spectrum.

Given an eigenvalue problem $\mathbf{A} v_k = \lambda_k v_k$ (assume real symmetric $\mathbf{A}$ and separated eigenvalues; although the latter is probably not necessary), let's attempt to estimate the smeared spectral density $$ S(\omega) = \sum_k \frac{\pi^{-1}\sigma}{\sigma^2 + (\lambda_k - \omega)^2} = \frac{\sigma}{\pi} \mathop{\mathrm{Tr}} [\sigma^2 + (\omega - \mathbf{A})^2]^{-1} $$ After hitting e.g. http://dx.doi.org/10.1016/0377-0427(96)00018-0 in a literature search, we know that an unbiased Monte Carlo estimator to the trace is $$ S(\omega) = \frac{\sigma}{\pi}\langle z^T [\sigma^2 + (\omega - \mathbf{A})^2]^{-1} z \rangle $$ where each entry of the random vector $z$ contains either $+1$ or $-1$ with probability 0.5 for each. For given $\sigma$ and $\omega$, the inverse product $[\sigma^2 + (\omega - \mathbf{A})^2]^{-1} z$ can be computed for example with the conjugate gradient method, or sparse LU on $[\omega + i \sigma - \mathbf{A}]^{-1} [\omega - i \sigma - \mathbf{A}]^{-1}$ to minimize fill-in. This allows estimation of $S(\omega)$ also for large matrices. In practice, it seems the CG solution doesn't need to be very accurate, and neither are many vectors necessary in computing the average. This may depend on the problem.

The above appears to weigh parts of the spectrum more evenly than a similarly smeared Krylov spectral density --- try diag(linspace(0, 1, 150000)) --- although maybe there is a way to correct for this?. This is somewhat similar to the pseudospectral approach, but the result indicates the (smeared) number of eigenvalues in the vicinity to point $\omega$, rather than the inverse distance to the nearest eigenvalue.

EDIT: A better performing alternative for computing the above quantity is to compute Chebyshev moments (via similar stochastic evaluation as above) and then reconstruct the spectral density from them. This requires neither matrix inversions nor separate computations for each $\omega$. See http://theorie2.physik.uni-greifswald.de/downloads/publications/LNP_chapter19.pdf and references therein.


See the paper "On Sampling-based Approximate Spectral Decomposition" by Sanjiv Kumar, Mehryar Mohri & Ameet Talwalkar (ICML 2009.). It uses sampling of columns of your matrix.

Since your matrix is symmetric you should do the following:

Let A be your n*n matrix. You want to reduce the computation of the eigenvalues of an n*n matrix to the computation of the eigenvalues of an k*k matrix. First choose your value of k. Let's say you choose k=500, since you can easily compute the eigenvalues of a 500*500 matrix. Then, randomly choose k columns of the matrix A. Contruct the matrix B that keeps only these columns, and the corresponding rows.

B = A(x,x) for a random set of k indexes x

B is now a k*k matrix. Compute the eigenvalues of B, and multiply them by (n/k). You now have k values which are approximately distributed like the n eigenvalues of A. Note that you get only k values, not n, but their distribution will be correct (up to the fact that they are an approximation).


You can always use the Gershgorin circle Theorem bounds to approximate the eigenvalues.

If the off-diagonal terms are small, the diagonal itself is a good approximation of the spectrum. Otherwise if you end up with an approximation of the eigenspace (by other methods) you could try to express the diagonal entries in this system. This will lead to a matrix with smaller off-diagonal terms and the new diagonal will be a better approximation of the spectrum.

  • $\begingroup$ Gerschgoring gives no aprroximations but error bounds, so is irrelevant here. Moreover, using your method on a sparse matrix would require a dense eigenvector matrix, which is impossible to store for the OPs problem. $\endgroup$ – Arnold Neumaier Jun 22 '12 at 12:17
  • $\begingroup$ As I said, the diagonal is itself an approximation of the spectrum with the error bounds given by the Gershgorin circle theorem, of course Gershgorin error bounds are not approximations. The diagonal will be a good approximation if the off-diagonal terms are small, wich I believe is the case since OP said that matrix is sparse. $\endgroup$ – FKaria Jun 22 '12 at 17:50
  • 5
    $\begingroup$ Most sparse matrices arising in practice have some significant off-diagonal elements in each row and column, which makes the diagonal very poor approximations (e.g., for a Laplacian of a regular graph the diagonal is constant), and the error bounds useless. $\endgroup$ – Arnold Neumaier Jun 22 '12 at 19:38

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.