[Math] Fast algorithm for approximating eigenvalue distribution of large sparse matrix

eigenvalues-eigenvectorsMATLABmatricesnumerical linear algebrasparse matrices

I am interested in the eigenvalue distribution of a huge $2^{16} \times 2^{16}$ Hermitian sparse matrix with spectrum contained in $[-1,1]$. That is I don't need to know all eigenvalues exactly, but rather the approximate number of eigenvalues in, say, the intervals $[-1,-0.99],\dots,[0.99,1]$.

The Matlab command eig fails since it doesn't accept sparse matrices and the matrix is too big for being stored as a normal matrix. The command eigs doesn't help me very much since it only gives me the $k$ biggest eigenvalues and takes forever.

Are the fast approximate algorithms for approximating the spectral density?

Best Answer

You can apply a common approach known as the spectrum slicing. Assume that $M$ is a Hermitian matrix factored as $$\tag{1}M=LDL^*$$ (a complex variant of the $LDL^T$ factorisation), where $L$ is unit lower triangular and $D$ diagonal. The inertia (the number of negative, zero, and positive eigenvalues) of $M$ and $D$ is the same and since $D$ is diagonal, you can get the inertia of $M$ pretty easily.

Now fix a $\sigma\in\mathbb{R}$. The eigenvalues of $M-\sigma I$ (where $I$ is the identity matrix) are the eigenvalues of $M$ minus $\sigma$. Consider the factorisation $$ M-\sigma I=L_{\sigma}D_{\sigma}L_{\sigma}^*. $$ Both $M-\sigma I$ and $D_{\sigma}$ have again the same inertia so, e.g., the number of negative diagonal entries of $D_{\sigma}$ determines the number of eigenvalues of $M$ smaller than $\sigma$.

In your case, since you know that the spectrum of $M$ is contained in the interval $[-1,1]$, you need to perform only two factorisations to determine how many eigenvalues are contained in the intervals $[-1,-\alpha)$ and $(\alpha,1]$ for some $\alpha\in(0,1)$. In particular, $$ M+\alpha I=L_+D_+L_+^*, \quad\text{and}\quad M-\alpha I=L_-D_-L_-^*. $$ Then, the number of negative diagonal entries in $D_+$ gives you the number of eigenvalues of $M$ smaller than $-\alpha$ {that is, in the interval $[-1,-\alpha)$} and the number of positive diagonal entries of $D_-$ gives the number of eigenvalues of $M$ greater than $\alpha$ {that is, in the interval $(\alpha,1]$}.

In MATLAB, the factorisation you look for is implemented in the function ldl.

Actually, you could also use eigs to compute, say, $k$ eigenvalues largest in magnitude (that is, close to the boundaries of the interval $[-1,1]$) and check if the result gave you an eigenvalue with a magnitude smaller than $\alpha$ with both signs). If it is so, then the result of eigs contains all eigenvalues in $[-1,-\alpha)$ and $(\alpha,1]$ plus some eigenvalues outside of these intervals. If this criterion is not satisfies, you need to increase $k$. Of course, this could be rather costly especially if the eigenvalues would be clustered around $-1$ and $1$ which would require to use relatively high $k$.

P.S.: The matrix $D$ is not generally diagonal, but block diagonal with $1\times 1$ and $2\times 2$ blocks. Still, the eigenvalues of $D$ (and hence their signs which matter here) can be easily obtained.

Related Question