Efficient way to diagonalize a sparse, block tridiagonal matrix

Question

Is there an efficient way to diagonalize a block tridiagonal NxN matrix? The matrix is made up of blocks of size K with N = 2 * K * K.

I need an algorithm that scales way better than the standard O(N^3) scaling, since I am usually interested in N ca. 10.000 and K ca. 100 and need to diagonalize a lot of them. I need both the eigenvalues AND the eigenvectors, if possible to full precision.

Further information, if useful: The matrix I am considering is hermitian, the off-diagonal matrices are just the unit matrix and the matrices on the main diagonal are pentagonal. I implemented the exact form of my matrix in this minimal reproducable example (python):

import scipy.linalg as linalg
import numpy as np

#N is usually in this form and I want to be able to have a code that is still fast, if  I set the 5 to a 50.
N = 2* 5**2

#Creating Matrix in my form...

N_band_max = int(np.sqrt(N/2))*2

#setting bands...
Main = np.random.random(N)
Diag1 = np.random.random(N-1) + 1j*np.random.random(N-1)
Diag2 = np.ones(N-2)
Diag3 = np.ones(N-N_band_max) 

#This is for the block tridiagonal structure

for i in range(len(Diag2)):
    if(i%N_band_max == N_band_max-2 or i%N_band_max == N_band_max-1):
        Diag2[i] = 0.
        
for i in range(len(Diag1)):
    if(i%2 == 1 ):
        Diag1[i] = 0.


H = np.diag(Main) + np.diag(Diag1,1) + np.diag(Diag1.conjugate(),-1) + np.diag(Diag2,2) + np.diag(Diag2.conjugate(),-2) + np.diag(Diag3,N_band_max) + np.diag(Diag3.conjugate(),-N_band_max)

#Up to here, I just created a matrix in the desired form, for reproducibility.
#This is the part I want to be FASTER. Here, any speedup would be highly appreciated
EigVals,EigVecs = linalg.eigh(H)

I'd be grateful for any help. I've also read this paper:

https://arxiv.org/abs/1306.0217

but while it in principle provides an algorithm for my problem, I was not able to implement it efficiently (in python) and the authors did not publish their code.

(The same question from a mathematical standpoint and with a nice latex format is here: https://mathoverflow.net/questions/427654/diagonalizing-a-block-tridiagonal-matrix)

Efficient way to diagonalize a sparse, block tridiagonal matrix

0 Answers