hermitian_lanczos and loss of orthogonality ?

[loiseaujc]

Hi,

I'm working on a package called LightKrylov. It is quite similar to Krylov.jl except that it is developed for Fortran. The philosophy is quite the same: provide some abstract implementations of standard Krylov processes and Krylov-based solvers which people can use simply by extending the abstract vector and linear operator types to accommodate for their own definition/data structure.

I'd like to say that Krylov.jl served as the inspiration, but to fairly honest I found about it quite recently. Nonetheless, I'd like to congratulates you guys. It looks like an awesome package! In comparison, LightKrylov is nowhere as feature-complete yet. For the moment, we focused on the most widely used techniques which, in my opinion are:

• Krylov processes : Hermitian Lanczos, Arnoldi and Golub-Kahan.
• Krylov solver : GMRES (and FGMRES) for nonsymmetric systems and Conjugate Gradient for sym. pos. def. ones.

When looking at your code base, I couldn't help but notice that you are not doing any kind of re-orthogonalization (whether full or partial) of the Krylov bases for the symmetric and non-symmetric Lanczos factorizations. I was wondering if there are any particular reasons for that? Don't your implementations suffer from a gradual loss of (bi-) orthogonality as the Krylov process proceeds?

[amontoison]

Hi @loiseaujc!
Thanks for your kind message.

The Julia functions that implement Krylov processes, like hermitian_lanczos, arnoldi, or golub_kahan, are not directly used in the Krylov methods. We use them to easily verify properties or prototype new Krylov methods. However, I could easily add options to perform partial re-orthogonalization because I store all the Krylov basis vectors.

For the Krylov methods, it will be more complex because the storage requirements in the workspaces (KrylovSolver) will differ depending on this option.

The good news is that we have Krylov implementations based on the incomplete Arnoldi process (DIOM and DQGMRES), and we can say that these methods are similar to CG and MINRES with partial re-orthogonalization.
DIOM is an interpolation between CG and FOM, while DQGMRES interpolates between MINRES and GMRES.

Regarding bi-orthogonality, the non-Hermitian Lanczos process is inherently unstable, and I have never investigated whether partial re-orthogonalization could help. Adding a partial_reorthogonalization option in the Krylov process nonhermitian_lanczos could provide more insight into its impact.

Partial/full re-orthogonalization comes at a cost in terms of storage and computation, and it is not always necessary. It depends on the problem. In general, if we need it, it likely means that a preconditioner was not used or the problem was not properly scaled.

Update: I added local reorthogonalization (special case of partial reorthogonalization) in the PR ##957 for hermitian_lanczos.
Partial reorthogonalization is complex to support because it changes the structure of the projection in the Krylov basis (not anymore a triadiagonal).