Finding the distribution (histogram) of eigenvalues for large sparse matrices

Question

Are there any existing programs that are able to compute the (approximate) distribution of eigenvalues for very large (symmetric) sparse matrices?

Note that I do not need the eigenvalues themselves, only their distribution (finding all the eigenvalues is a more difficult problem), and I am primarily looking for existing software, not only the description of algorithms that I need to implement myself. My matrices are adjacency matrices of undirected graphs (or closely related matrices).

I found some papers discussing how this can be done in practice, so I was hoping that a working implementation exists somewhere.

Answer 1

In structural mechanics the number of eigenvalues of a matrix $K$ in a given range $(\alpha,\beta)$ is computed via the "Sturm sequence check", i. e. computing the $LDL^T$ factorizations of $K-\alpha I$ and $K-\beta I$ and counting the difference in the number of negative pivots.

If you have reasonably large bins, can be applied to your problem, and should be pretty straightforward to implement.

(A search on Lanczos shifted block algorithms should give more info, since this technique is often used in that context to check for missed eigenvalue/eigenvector pairs.).

This count is exact, and expensive for large $K$, so your request of an "estimate" or approximate count is still open. Please post any finding.

Edit/update

The authors of "A Padé-based factorization-free algorithm for identifying the eigenvalues missed by a generalized symmetric eigensolver" claim

To the best of the authors’ knowledge, no post-processing technique that does not require the factorization of a matrix related to $K$ and/or $M$ is currently available for checking whether an eigensolver applied to the solution of problem (1) has missed some eigenvalues in an arbitrary range of interest $[ \sigma_L , \sigma_R ]$.

This means that for obtaining an exact count of the eigenvalues in a given bin you have to use the classical Sturm sequence method (see also Sylvester inertia law).

General advice for implementing this approach in your case cannot be given, without an analysis of the properties of your adjacency matrices (dimensions, number on non zero entries, fill-ins after reordering, condition number of principal minors...).

Nevertheless I would suggest starting with a simple no-brainer approach, and see if you experience breakdowns of the implementation (assuming that this computation is not mission critical). I suggest to use the wonderful SuiteSparse by Tim Davis.

1. reorder your matrix to reduce fill-in, (e.g. calling SYMAMD from COLAMD on $A - \alpha_0 I$).
2. compute $L_i D_i L_i^T = A - \alpha_i I$, where $\alpha_i$, $i=0\dots m$ are the boundaries of your bins. (Try first a simple implementation like LDL without pivoting, and go for a pivoting approach only if you experience numerical difficulties.) (Note that without pivoting the symbolic factorization step can be recycled for all $i$.)
3. the count of negative diagonal terms in $D$ gives you the number of eigenvalues $\lambda < \alpha_i$.

This approach is effective only if $m$ is small or the factorization time is negligible with respect to eigendecomposition time: you have to perform some experiments to find out. Good luck.

Answer 2

I would imagine that something of this kind can be inferred from the pseudospectrum function. See the work by Marc Embree, for example: http://www.caam.rice.edu/~embree/