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.
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.)
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.
