search

JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
407 stars 106 forks source link

Stability of minimal residual part in BiCGStab(l) #190

Open haampie opened 6 years ago

haampie commented 6 years ago

The minimal residual part of BiCGStab(l) is unstable since the normal equations M'M x = M'b are solved by explicitly forming M'M. It seems that JuliaCI fails as a result of it for it, even though Travis does not (which is odd).

The original paper "BiCGStab(l) for linear equations involving unsymmetric matrices with complex spectrum" proposes Gram-Schmidt for orthogonalization similar as in GMRES, but section 3.4 actually advocates working with M'M because the dimension is small (l = 2 or 4 typically) and it saves a bunch of inner products.

I have not tried yet if replacing it by QR'ing first actually helps, but I figure it will. Also I'm not sure what the "standard" implementation of BiCGStab(l) is. Probably there is a whole zoo of variants.

haampie commented 6 years ago

Actually the paper "BiCGstab(l) and other hybrid Bi-CG methods" by Sleijpen, Fokkema and Van der Vorst discusses the accuracy vs speed trade-off, but concludes that it is hard to know what is best (obviously).

It proposes Bunch-Kaufman as a more stable alternative to LU / Cholesky when inverting M'M. Maybe worth a shot as well.