Make Lanczos algorithm more robust

[JohanSchott]

The Lanczos algorithm may suffer from numerical instability, see discussion in e.g. https://en.wikipedia.org/wiki/Lanczos_algorithm.
The generous unit-test tolerances indicate numerical round-off errors are indeed present (compared the same code on Mac OS X and Ubuntu).
Until now, the errors I have seen have been smaller than the physics of interest.
But it is good to be aware of this numerical instability issue, in particular if very subtle features are of interest.

In https://en.wikipedia.org/wiki/Lanczos_algorithm, a few ideas and references are provided of how to improve the numerical stability of the Lanczos algorithm, that might be relevant for this repository.

[johanjoensson]

I have done some investigation of a few methods.

1. Using blocks of vectors for the Lanczos method, instead of a single vector. This is useful if we want to calculate offdiagonal elements in the Greens function. Numerically, this increases stability a bit as well.
2. Partial reorthogonalization is probably the best compromise between computational cost and numerical accuracy, link to old paper
3. Full reorthogonalization is easy to implement, but veeeery computationally expensive.

In general, all methods that adds some orthogonalization-step to the method will require that we save all the Krylov vectors we generate. This increases the memory footprint of the calculation quite a bit. In my experience, the Greens function does not converge faster when we add orthogonalization. I have not actually looked at the numerical accuracy of the resulting Greens function though.

[JohanSchott]

Restarting the algorithm after a certain number of iterations seems to be a popular approach, e.g. Thick-Restart Lanczos method: https://zenodo.org/records/1236144

[JohanSchott]

Not sure how big this issue is.

[johanjoensson]

No, I am not sure that this will actually help with calculating the Greens function. The restarted methods work by focusing on certain eigenvalues/eigenvectors. If what you want is just certain eigenvalues/eigenvectors this is great, but for the Greens function to converge quickly, what eigenvalues should we focus on?

[JohanSchott]

I don't know. Perhaps it's not suitable for calculating Green's functions.

[johanjoensson]

The way we calculate Greens functions relies on the tridiagonalization of the Hamiltonian. With restarted methods I think we replace the tridiagonal structure with something different, but maybe there is a clever way to calculate the Greens function (we only need one element of the inverse of the Hamiltonian).

[JohanSchott]

When I quickly read about restarted methods I agree it does seem the matrix is not tridiagonal. So probably the three implementation ideas you outlined above are more promising.