Tolerance seems to be ignored and there's no absolute tolerance?

[ChrisRackauckas]

N = 100
A = rand(Float64, N, N)
b = rand(Float64, N)
x = rand(Float64, N)
gmres!(x,A,b,tol=Float64(1e-60))
A*x - b
mean(abs.(A*x - b)) # 0.49568591270147505

gmres!(x,A,b,tol=Float64(1e5))
A*x - b
mean(abs.(A*x - b)) # 0.495852450033349

[haampie]

Thanks for opening the issue.

There is some issues with the example, but we can discuss how to improve the user experience here.

About the example: The expected value of A is fill(0.5, N, N) -- a rank one matrix with N - 1 eigenvalues at zero and 1 eigenvalue at N/2. The matrix A itself is some perturbation of this matrix. GMRES can't do much with this.

using Plots, LinearAlgebra
λs = eigvals(rand(100, 100))
scatter(real(λs), imag(λs), xlims = (-5, 52), ylims = (-25, 25), label = "eigvals")

At step k the residual vector of GMRES is rₖ = pₖ(A)r₀ where r₀ is the initial residual and pₖ is a kth order polynomial satisfying pₖ(0) = 1. To minimize the residual norm this polynomial should have zeros close to all eigenvalues. To converge quickly it should be a low-order polynomial. Judging from the distribution of the eigenvalues in the plot combined with the condition that pₖ(0) = 1, it is obvious that this is pretty hard. Also because GMRES restarts to ensure fixed memory demands it only builds a polynomial of a fixed order (by default 20 I think); it's just impossible to converge. In fact you really need a polynomial of order 100 to reduce the residual norm, but that defies the purpose of an iterative method completely:

> x, hist = gmres(rand(100,100), rand(100), restart=100, tol=1e-7, log=true)
> hist
Converged after 100 iterations.
> plot(hist, yscale = :log10, mark = :+, legend = :bottom, xlabel = "Iteration number")

If you try restart = 20 and maxiter = 10000 you'll see that it will just stagnate.

Better example: if you take A = rand(N, N) + N*I, then N-1 eigenvalues are clustered around N and the last eigenvalue is close to 1.5N. Now it's super easy to construct a polynomial that's 1 in the origin and zero near all eigenvalues around N and 1.5N :)

> x, hist = gmres(rand(100,100)+100I, rand(100), tol=1e-7, log=true);
> hist
Converged after 6 iterations.
> plot(hist, yscale = :log10, mark = :+, legend = :bottom, xlabel = "Iteration number")

---

there's no absolute tolerance

Yes, there's only a relative tolerance, maybe we should have atol and rtol everywhere.

Tolerance seems to be ignored

This is not true, it's just that the specified tolerance can never be attained.

---

So maybe the best thing for the user experience is to put some @warnings if the method has not converged to the specified tolerance and set log=true by default.

[ChrisRackauckas]

I see, great analysis. I guess this is a good example where simply using a rand input isn't a good check 👍

[jiahao]

Random matrices have special structure - it happens more often than not that they are terrible test candidates.

[ChrisRackauckas]

I'll have to find the time to sit in on a random matrix lecture some day.