Computing eigendecomposition of extremely large (though sparse band) matrix

Question

I have a very large graph (about over a billion nodes) and I need to compute all the eigenvectors and eigenvalues of the graph (i.e., of the Laplacian of the adjacency matrix) for downstream analysis. Naively, this could take centuries to run and require petabytes of storage space. However, the graph is extremely locally connected (and symmetric). That is, I know that nodes are only connected to the few nodes around it. For example, node 5,000 would never be connected to node 10,000. Not only are the extreme vast majority of the values in the matrix zero, but all of the non-zero values are close to the diagonal (local connectivity).

How can I compute the eigendecomposition in a way that doesn’t take a very long time?

Related: Compute all eigenvalues of a very big and very sparse adjacency matrix

Potentially useful: DSBEVD, DSBTRD

I believe this is a "band matrix."

Here's about what the matrix would look like, except much larger:

Note: I am wondering about an algorithm or method that will accomplish this, though a language-specific implementation is still fine, since I can at least read how it is implemented.

Answer 1

The LAPACK function dsbevd does this. I will describe the general logic of what it is doing.

The band matrix is transformed into a tridiagonal matrix, which is one that has non-zero elements only on the diagonal and the first off-diagonals, making it much simpler to handle computationally.

Then, it uses a series of orthogonal transformations (Givens rotations or Householder reflections) to zero out elements outside the tridiagonal band in a way that preserves the matrix’s symmetry, so the that the eigenvalues remain consistent with the original matrix.

If only eigenvalues are needed, the algorithm uses an iterative method, (QL or QR algorithm), to find the eigenvalues of the tridiagonal matrix. These methods work by iteratively refining the matrix and ensuring that off-diagonal elements tend toward zero, leaving the diagonal elements as the eigenvalues.

If both eigenvalues and eigenvectors are required, a more advanced method called divide-and-conquer is used, which splits the matrix into smaller, independent submatrices (based on where off-diagonal elements are small or negligible) and solves each subproblem separately. The results from these subproblems are then combined to give the full set of eigenvalues and eigenvectors.

If eigenvectors are computed, they are initially found for the tridiagonal matrix. However, since the original matrix was transformed during the reduction to tridiagonal form, the eigenvectors must be transformed back to match the original matrix. This is done by applying the same orthogonal transformations used during the reduction step in reverse order.