# Finding ALL Eingenvalues of a Sparse Integer Matrix

I would like to find ALL Eingenvalues of a huge, very sparse integer matrix. This matrix has a lot of known properties, e.g. that it is symmetric and nearly tridiagonal, with very few (max. ca. 4 per row for max. ca. 10% of the rows) entries off the super- and diagonal. The matrix is very large (for scientific reasons I intend to go to matrix sizes limited by the memory available for computations; let's say 1 TB for the dense matrix), and contains only integers.

Concrete example to have numbers: a matrix with 214'000 rows, of which the diagonal, the superdiagonal and 4 additional entries for every 10th row is filled (ca. 0.001% of entries are non-zero, 0.0002% off the tridiagonal).

With all this knowledge, I would expect to find optimised algorithms to find these eigenvalues. However, so far, the answers I find either make use of the sparsity and implement an algorithm from e.g. the Spectra library, which in turn should not be used for finding all Eigenvalues, or then making the matrix dense and using directly an algorithm from e.g. LAPACK. There, algorithms exist that deal with tridiagonal matrices, but no sparse method to tridiagonalize a matrix.

I feel like the answer to this question will be "your intuition is wrong", but I am genuinely puzzled why there doesn't seem to be any algorithms that make use of either the fact that Eigenvalues of an integer matrix must be algebraic integers, or a matrix being so close to tridiagonal, given the performance of finding the Eigenvalues of a Hessenberg matrix.

So, please, either educate me on which algorithm I have not found during my research, or help me understand why "nearly tridiagonal" is in no way better than "dense".

New contributor
BernhardWebstudio is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• What do you mean by integers here? The usual $\mathbb{Z}$? Algebraic integers? Because "eigenvalues of an integer matrix must be integers" surely doesn't hold for $\mathbb{Z}$. Commented Jul 9 at 16:50
• "eigenvalues of an integer matrix must be integers" does not hold - counter example is a symmetric circulant matrix with 1 on the diagonal, and only 2 other non-zero entries on each row/column which can also be set to 1. Since for an order N matrix the evals are related to the nth roots of unity, they clearly are not integers. Commented Jul 9 at 17:04
• What do you mean by "very large"? Rough order of magnitude? Commented Jul 9 at 17:05
• The Wikipedia article mentions that the eigenvalues are algebraic integers. In general, they are complex numbers. Commented Jul 10 at 12:13
• I don't think it is realistic to compute all eigenvalues of matrices of this size. If you want all eigenvalues, the sparsity of the original matrix does not buy you very much because you will have to build eigenspaces of size equal to the matrix dimension, and the representations of the eigenspace is dense. The fact that the matrix has integer entries also does not buy you anything because the eigenvalues and eigenvectors are not integers -- all you know is that they are algebraic numbers, but in general not even rational. Commented Jul 10 at 19:36

## 2 Answers

For something like an eigensolver, a 1% fill isn't really sparse. During the tridiagonalization step, the solver has to fill in some of the zeros each time it eliminates one of the off-digonal zeros. This fill-in will in turn introduce more fill in, until your trailing matrix is 100% dense.

If you want to test it for yourself, you can tridiagonalizing your matrices with dsytrd and look at how many zeros are in your Householder reflectors. The first ones will be sparse, but they'll become increasingly dense as the factorization progresses.

You can maybe do something complicated where you use a sparse datastructure for the first few hundred columns/rows. But, I'm skeptical that it would provide a benefit given the complexity it would require. Writing good sparse factorization routines is fairly hard.

My suggestion would be to try to get access to a cluster or supercomputer and just use a dense library. n=250,000 is doable with a few dozen GPUs, and n=2,000,000 should be doable on a large GPU cluster or supercomputer. I know there are quantum chemistry apps that routinely run million-by-million Hermitian eigensolvers on the DOE supercomputers. I've run a dense SVD with singular vectors of size n=100,000 on 8 nodes of the Summit supercomputer (48 V100 GPUs) in a few minutes, and a values-only symmetric eigensolve will be faster.

### Update for a sparsity of 4 non-tridiagonal elements per 10 rows

This is the level of sparsity that I would expect a sparse tridiagonalization routine to be beneficial. You'll still have a large amount of fill-in, and you'll reach a point in which it's cheaper to switch to a dense routine. Although, the exact performance will depend heavily on the exact fill structure.

Doing some quick estimation, I think the first few thousand columns would have a manageable amount of fill. But, I'd guess you'd reach a basically dense matrix somewhere between a quarter to two-thirds of the way through the matrix (but it depends a lot on the matrix's structure).

I would suggest running dsytrd on one of your matrices, and (assuming you're using lower storage) plot the number of nonzeros below the tridiagonal in each column. That won't tell the whole story, but that will tell you the maximum number of elements that each column can fill in the trailing columns. I'm guessing the optimal point to switch to a dense solver would be somewhere between 1% and 20%.

To actually answer the question you asked: I'm guessing such a routine isn't common because it's a niche usecase. If the matrix isn't sparse enough, fill-in makes the matrix dense anyways. If you only need a few eigenvalues, iterative methods are faster. And if you need the eigenvectors, those will be dense anyways.

I did realize that Matlab actually provides a sparse eigensolver. It can't compute eigenvectors, it only supports real symmetric matrices, and I can't remember how optimized it is. But, if you can use Matlab, that might be useful.

• Thank you for your answer. I have edited my question with more precise information about the sparsity; it is actually more in the range of 0.001% (including the entries on the super- and diagonal) for the system with 250'000 rows, but I suspect your answer stays the same. Commented 2 days ago
• Actually, at that sparsity, I suspect there's likely be a benefit for a sparse tridiagonalization routine. I'll update my answer. Commented 2 days ago

You are right that there is nothing off-the-shelf for this job. Computing all eigenvalues of a large matrix while keeping sparsity is not easy, and usually one can get away with computing the first eigenvalues, the last, or those closest to a specific point (is this your case, maybe?).

Similarly, I know of no algorithms that use integrality and algebraic integers: working with characteristic polynomials is ill-conditioned, so it isn't a common choice.

If you just want to get the job done once, getting access to HPC cluster as suggested in another answer might be the fastest solution.

If you want to meddle directly with numerical code, I suggest trying a divide-and-conquer algorithm. Your matrix is tridiagonal-plus-rank $$\frac{n}{10}$$, so you can compute the eigenvalues of the tridiagonal part, and then treat the low rank part as the composition of many successive rank-1 modifications of this matrix. Rootfinding on the secular equation can be used to compute the eigenvalues of these rank-1 modifications. The total time cost is cubic, $$n^2 \times \frac{n}{10}$$, not parallelizable since you have to solve $$\frac{n}{10}$$ secular equation in sequence, but the space cost is just $$O(n)$$, so this algorithm is very cache-friendly and I would imagine it to be faster than the classical dense algorithms.

• Given that D&C is sometimes used for finding the eigenvalues of the tridiagonal part, one could maybe just do D&C on the whole matrix. That might provide some parallelism, depending on the exact structure. Commented 2 days ago
• You might be able to do a structure preserving implicit QR. I think polynomial roots are often found using a unitary plus rank-1 representation of the companion matrix. Tridiagonal plus rank-k might just work. Commented yesterday