What and Why
As an advanced technique, we could implement the Lanczos algorithm originally scheduled for [#168]. There is a branch lanczos-attempt that holds a first implementation; however, there are some caveats that are not implemented but mandatory.
Note that the Lanczos algorithm is numerically unstable; in particular, errors w.r.t. the orthogonalization seem to accumulate rather quickly, reaching scalar products with values of 1e-5 for supposedly orthogonal vectors already with ~30 basis functions. This requires a careful implementation, but also impacts the usability: How to produce code that produces somewhat reliable results and exposes the issues without overwhelming the user?
- It seems rather crucial for convergence that you always orthogonalize first, then apply the Hamiltonian again. Wikipedia says here that this is the most stable approach; how true...
- It is mandatory that you exploit the tridiagonal form in the orthogonalization, otherwise errors accumulate (even) faster
- It may be necessary to improve the scalar product / normalization scheme for the tensor library. The current implementation is pretty direct, maybe one could specialize standard forms to corresponding BLAS functions?
- After some time, you will run into problems where the Krylov space does not increase in size any more. This must be handled gracefully (the first attempt threw exceptions for simplicity => not very user-friendly).