I want to run a block Lanczos block-tridiagonalization on a hermitian, sparse matrix (of relatively small size $\sim 10^2 \times 10^2$). However the matrix typically has many eigenvalues that are highly degenerate, and the eigenvalues are closely spaced. I only care about the invariant subspace of the initial vectors, so it is desired for the procedure to terminate when the elements of the subdiagonal become close to zero.

However, when I try a block Lanczos tridiagonalization in this way, it is numerically unstable because of the above mentioned eigenvalue structure (as e.g. discussed in the answer here). I am already re-orthogonalising all vectors at each step of the Lanczos procedure.

Are there any techniques with which to make the block Lanczos iteration more stable under these circumstances? Is this already implemented in a code?

  • 5
    $\begingroup$ At that size you might as well switch to a dense algorithm (readily available in LAPACK), if Lanczos is giving you trouble. $\endgroup$ Feb 12, 2023 at 16:55
  • $\begingroup$ Thank you @rchilton1980! Yes indeed, the matrix size is not much of an issue - I'm using block Lanczos, because I require the Krylov basis / the block tridiagonalisation similarity matrix of the Krylov basis (this stems from a physics problem, specifically this Block-Lanczos Density Matrix Renormalisation Group technique) - my apologies for the naive question, but were you referring to an alternative, dense algorithm giving this as well? Or were you thinking of a entirely different algorithm, for e.g. just eigenvalues? $\endgroup$
    – lm1909
    Feb 13, 2023 at 11:40

1 Answer 1


I'll expand into a full answer to get more space.

The algorithms from LAPACK can provide both the eigenvectors and the eigenvalues. On symmetric/hermitian systems, like those for which Lanczos is applicable, these algorithms are very robust. I'd suggest dsyev.f or cheev.f as good points of entry for symmetric/hermitian systems, respectively. There might be banded ones too, if you poke around, though I'd bet this is not worth the effort unless you start looking at larger instances.

The cited paper is a lot to parse, so I may be understating the problem. However it's hard to imagine a situation where the Krylov basis would somehow be more useful than the eigenvectors themselves. If you want to compute the invariant subspace associated with some arbitrary (initial) vector, you can always just project it onto the eigenvectors, they're all orthogonal.

EDIT: You might also consider toying with this quickly in matlab, to see if you can postprocess your needed subspace information from the eigenvectors themselves. Here the routine of interest is [V,D] = eig(A)

  • $\begingroup$ +1 TIL defective symmetric/Hermitian matrices still have orthogonal eigenvectors (spectral theorem here: en.wikipedia.org/wiki/… ) $\endgroup$
    – user20857
    Feb 14, 2023 at 5:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.