Finding eigenvalues of a complex symmetric tridiagonal matrix

Question

I am trying to find specific eigenvalues and -vectors of a large complex symmetric tridiagonal matrix (at least 10000x10000, and ideally larger). I know roughly which eigenvalues I am looking for, so I've been using scipy.linalg.sparse.eigs (really ARPACK) to find them. Unfortunately, this has turned out to be slow for large matrices - solving a 7000x7000 matrix took approximately 30 minutes.

I feel that there must be some way of using the symmetry properties of the matrix I am dealing with - it's not just any sparse complex matrix. As a comparison, when I use hermitian matrices of the same size, linalg.sparse.eigsh allows me to reduce the solution time by approximately a factor of 2.

Are there any freely available eigensolvers out there which implement efficient algorithms for diagonalising the kinds of matrices I am dealing with? Do such algorithms even exist?

Answer 1

You can apply a banded cholesky-like algorithm to generate L*L.', although the stability is questionable because you cannot bound the pivot growth (unlike the LL' positive definite case, which is well behaved). So perhaps a shift-invert strategy (since you know which eigenvalues you want) might be able to give you an answer. In the event you have to slog through a cluster to find a particular eigenvector nearby your chosen shift, you can exploit the fact that distinct eigenvectors (vi,vj) of a complex-symmetric matrix are orthogonal under a quirky non-hermitian "inner product", vi.'*vj = 0. You can use that fact to setup an outer loop that deflates your arpack residual against converged but unwanted eigenvectors on every arnoldi step, via gram-schmidt. Note this part is dodgy too, because that "inner product" $<vi,vi>$ can be negative/zero, so your deflation/projection/normalization gram-schmidt coefficients can go haywire, too. Ultimately, complex symmetry is a nice storage-saving property but does not really impart the same robustness as real/hermitian symmetry.

So lots of caveats, but FWIW I have written working software based on this - wouldn't suggest it otherwise.

Answer 2

LAPACK seems to have several methods for this case. You would have to wrap the fortran code with f2py however. Either way, the algorithms certainly exist, and the LAPACK documentation provides references to the relevant articles.

On a slightly different note, make sure you're using a new version of gfortran or (better yet) Intels fortran compiler ifort. It might save you 20% or more (on runtime).