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

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.

• What specifically do you mean by their distribution? – Geoffrey Irving Jan 18 '13 at 23:47
• @GeoffreyIrving Like a histogram. How many eigenvalues are there in a given interval $[a,b)$, or more practically, how many are there between [0,1), between [1,2), etc. Does this answer your question? – Szabolcs Jan 19 '13 at 0:13
• Does the term "approximate" refer to the binning of the eigenvalues (the larger the bins, the coarser the approx of the true distribution), or are you also looking to just an estimate/bound of the count of the eigenvalues in a bin? – Stefano M Jan 20 '13 at 16:18
• @Stefano The second one. I'm also looking for an estimate if the estimate is considerably easier (= faster) to compute for very large sparse matrices. – Szabolcs Jan 20 '13 at 16:54

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

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.

• Thanks for the pointers! I also found a reference to Sturm sequences (the one I linked in the question). Do you know of any software that already implements these things in a way that's usable for a sparse matrix of size larger than 50000? What I was really hoping for was a ready made software so I do not need to completely understand implement the method myself. If there's no existing public software, of course I'll have to do that, but it would be really good to be able to save some time ... – Szabolcs Jan 20 '13 at 16:55
• @Szabolcs I'm no expert in this area, but "Sturm seq. check" is common in large scale eigenvalue solvers for structural mechanics problems, so you can try to extrapolate code from there. You can also search for libraries for computing "modal density", which is what you are looking for in structural mechanics lingo. – Stefano M Jan 20 '13 at 21:10
• Thanks for the edit, I just noticed it. There's a lot of useful information here, I'll need some time to go through it. – Szabolcs Jan 23 '13 at 0:08

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/

• Computing pseudospectra is more expensive than computing eigenvalues. There are good reasons to compute pseudospectra, but it's different information from a "histogram" and not fast. – Jed Brown Jan 19 '13 at 7:43
• The term 'pseudospectrum' is new to me, but what I found about it online seems to say that they're used to deal more easily with non-Hermitian matrices. All my matrices are Hermitian here. Can this still be of use? – Szabolcs Jan 19 '13 at 15:59
• If it's applicable to a general class of matrices (non-Hermitian) then it's of course also applicable to a subset of those (Hermitian). @JedBrown is probably right that it's expensive to compute -- I was simply pointing out that there may be information in the pseudospectrum that, possibly, could be used to evaluate the number of eigenvalues in a subset of the complex plane (or real axis) without actually having to compute the eigenvalues themselves, in much the same way as a contour integral in the complex plane informs about singularities of a function without having to know their locations. – Wolfgang Bangerth Jan 20 '13 at 17:02
• 2-norm $\epsilon$-pseudospectra of Hermitian matrices tend to be a little boring: disks of radius $\epsilon$ (in the complex plane) centered around the eigenvalues (aligned on the real axis). – Stefano M Jan 20 '13 at 23:33