9

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?

Rodrigo de Azevedo
  • 23,223
Dominik
  • 444

1 Answers1

6

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.

Algebraic Pavel
  • 23,818
  • Thank you for your answer! That is an interesting application of Sylvester's inertia theorem. However, I am not sure how fast this procedure is compared to eig. I don't get the ldL command in MATLAB, the $D$ it produces is not diagonal but block diagonal from which it isn't immediately possible to read of the positive eigenvalues. Secondly it does not work for sparse hermitian matrices, so the problem that the matrix is too big to be stored remains. What am I doing wrong? – Dominik Jun 02 '14 at 14:52
  • @Dominik Use [L,D,P]=ldl(M) or [L,D,p]=ldl(M,'vector') for sparse matrices, it also computes the permutation matrix from pivoting. For the blocks in $D$, it's true, I've made a small edit. Anyway, the signs of the eigenvalues of $D$ are easy to get the matrix has only small blocks up to size 2. It can still be much faster than eig as ldl uses a sparse elimination. BTW if you'd be interested only in the total number of eigenvalues in $[-1,-\alpha)\cup(\alpha,1]$, you can apply the similar approach to the matrix $M^TM$ (which, however, is generally less sparse than $M$). – Algebraic Pavel Jun 02 '14 at 15:27
  • Also, I'm not completely certain with that but I guess that the $2\times 2$ blocks actually have one positive and one negative eigenvalue. – Algebraic Pavel Jun 02 '14 at 15:29
  • I keep getting the error: Error using ldl: Complex sparse LDL is not yet available. I have the latest Matlab version installed. Thank you very much, anyway! I'll have a look into fast ldl algorithms. – Dominik Jun 03 '14 at 11:21
  • @Dominik Ok, I haven't noticed that the complex version is not yet implemented. What you can do is to write $M$ as $M_R+M_I$, where $M_R$ is the real part of $M$ and $M_I$ is the imaginary part. Since $M$ is Hermitian, $M_R$ is symmetric and $M_I$ is skew-symmetric. Then instead of $M$ consider the matrix $\underline{M}=\begin{bmatrix}M_R&M_I\-M_I&M_R\end{bmatrix}$. This matrix is symmetric (and sparse if $M$ is sparse) and real so you can use the Matlab's ldl function on it. Also, it has the same eigenvalues as $M$ except that the multiplicity of each eigenvalue of $M$ is doubled. – Algebraic Pavel Jun 03 '14 at 12:32
  • Of course, I forgot the imaginary unit. The splitting is actually $M=M_R+\iota M_I$, where $\iota^2=-1$ is the forgotten one. – Algebraic Pavel Jun 03 '14 at 12:59