Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GMRES Fails Silently From Stagnation #338

Open
packquickly opened this issue May 3, 2023 · 0 comments
Open

GMRES Fails Silently From Stagnation #338

packquickly opened this issue May 3, 2023 · 0 comments

Comments

@packquickly
Copy link

packquickly commented May 3, 2023

GMRES 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 (restart.) ie. roughly 90% of the time GMRES will fail silently on a problem of dimension 25 if restart = 24, and will still fail 90% of the time with a problem of dimension 100 if restart = 99. For a difference in dimension of 5, the test below returns an incorrect value greater than 99% of the time.

Reproducing code:

import IterativeSolvers: gmres
import LinearAlgebra: cond, norm

for problem_dim in [25, 100]
    for excess_dim in [0, 1, 5]
        succeeded = 0
        restart = problem_dim - excess_dim
        for i = 1:100 
            condition_number = Inf
            matrix = nothing
            while condition_number > 1000
                matrix = randn(Float64, (problem_dim, problem_dim))
                condition_number = cond(matrix)
            end

            true_vec = randn(Float64, (problem_dim,))
            b = matrix * true_vec

            julia_soln = gmres(
                matrix, b, restart=restart, log=false, abstol=1e-7
            )

            residual_norm = norm(julia_soln - true_vec)

            succeeded += (residual_norm < 1)  # very loose tolerance!
        end
        println("problem dim: $problem_dim. restart dim: $restart. succeeded: $succeeded")
    end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant