search

SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
245 stars 52 forks source link

Every GMRES Algorithm Fails Silently From Stagnation #297

Closed packquickly closed 11 months ago

packquickly commented 1 year ago

Every implemented GMRES algorithm often fails to converge when restart is less than the problem dimension, but will return without indicating failure.

This failure (likely due to stagnation) seems to depend on the absolute difference between the problem dimension and the Krylov subspace dimension (gmres_restart.) The rate of failure is slightly different for each algorithm, but the lowest failure rate resulting from a 1-dimensional difference in the example case I provided is 30%, ie. problem_dim = 100 and gmres_restart = 99 will silently fail and return an incorrect solution ~30% of the time. As will problem_dim=25 and gmres_restart=24. If the difference is 5 dimensions, that failure rate is 90% for KrylovKit, and >99% for the other two algorithms.

See the issues raised here for reproducing code: https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/issues/338 https://github.com/Jutho/KrylovKit.jl/issues/70 https://github.com/JuliaSmoothOptimizers/Krylov.jl/issues/729

ChrisRackauckas commented 1 year ago

It sounds like these are upstream issues?

packquickly commented 1 year ago

Not entirely.

  1. We probably don't want users to need to understand that GMRES can fail silently sometimes regardless of upstream implementation.
  2. Krylov.jl returns a stats object which we can use to detect this failure mode. However, when used to solve a linear system in LinearSolve.jl this information is not returned and the failure becomes silent.
packquickly commented 1 year ago

Actually, on further inspection Krylovkit.jl catches this as well (albeit does not raise an error, grr) and this information seems lost in the call to LinearSolve.jl

ChrisRackauckas commented 11 months ago

This was fixed back in https://github.com/SciML/LinearSolve.jl/pull/334