How are preconditioned GMRES relative residual computed?

I cannot reproduce the computed relative residual using preconditioned GMRES method. Here is the minimal reproducible example:

rng(0);
A = rand(1000,1000);
x = rand(1000,1);

Afun = @(x)(A*x);
b = Afun(x);

Ainv = inv(A+rand(1000,1000).*1e-4);
precond = @(x)(Ainv*x);

[Y,flag,relres,piter,resvec] = gmres(Afun, b, 32,1e-12,32, precond);

fprintf('  reported relres: %10.4e\n',relres)

% manually compute relres
manual_relres = norm(b - Afun(Y))/norm(b);
fprintf('  manual relres: %10.4e\n',manual_relres)

% manually compute preconditioned relres
manual_precond_relres = norm(precond(b - Afun(Y)))/norm(precond(b));
fprintf('  manual precond relres: %10.4e\n',manual_precond_relres)

The output of my result is:

reported relres: 3.2622e-13

manual relres: 3.3596e-13

manual precond relres: 3.3844e-11

Answer 1

[Heiko Weichelt]

Thanks for bringing this to our attention. I can reproduce the issue you mentioned @响 李.

We discussed this internally and found that there are mismatchs between the documentation of gmres, the help text in the file, and the algorithm itself.

The main issue is that we avoid un-necessary computations, if we converge. You can see this, that if you don't converge, like in @Abhinav Aravindan's answers, that are unfortunately using the preconditioner wrongly, the values actually match.

If you converge, however, we don't match and your relres seems to return , as we avoid applying the preconditioner to the current residual again.

Notice that the algorithm uses the provided stopping critera correctly, to compare against the relative residual (without the influence of the preconditioner).

We will discuss this internally further and will try to improve the involved pieces in the future so that documentation and values returned match better. Please let me know of any other questions you may have.

[SL: typo fix]

Answer 2

[Abhinav Aravindan]

Hi @响 李,

The inconsistency in the results arises from how the preconditioner matrix is defined. The gmres function specifies a preconditioner matrix "M" and computes "x" by solving the system "". The preconditioner matrix is to be defined as "M" and not "" which is likely causing the discrepancies in your results. You may modify the preconditioner matrix without the inversion in Line 8 and 9 as follows to obtain the expected results.

M = A+rand(1000,1000).*1e-4;
precond = @(x)(M*x);

Please refer to the below example for further clarity:

% Define Parameters
load west0479
A = west0479;
b = sum(A,2);
tol = 1e-12;
maxit = 20;

% Solution without Preconditioner Matrix
[x,fl0,rr0,it0,rv0] = gmres(A,b,32,tol,maxit);

% Calculate Relative Residual Error
manual_relres = norm((b-A*x))/norm(b);
fprintf('Solution without Preconditioner\n');
fprintf('  reported relres: %10.4e\n',rr0);
fprintf('  manual relres: %10.4e\n\n',manual_relres);

% Solution with Preconditioner Matrix
[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
M = L*U;
[x1,fl1,rr1,it1,rv1] = gmres(A,b,32,tol,maxit,M);

% Calculate Relative Residual Error
manual_precond_relres = norm(M\(b-A*x1))/norm(M\b);
fprintf('Solution with Preconditioner\n');
fprintf('  reported relres: %10.4e\n',rr1);
fprintf('  manual precond relres: %10.4e\n',manual_precond_relres);

Output:

Here is the documentation of gmres for your reference:

• https://www.mathworks.com/help/releases/R2023b/matlab/ref/gmres.html

I hope this helps resolve your query!

[响]

Thank you for your quick reply! But I think I am referring to passing a "function handle" as "M" rather than a matrix.

[Abhinav Aravindan]

Hi @响 李,

The gmres function seems to works as expected even while using a function handle as the preconditioner matrix, accomodating the change mentioned in my answer:

rng(0);
A = rand(1000,1000);
x = rand(1000,1);

Afun = @(x)(A*x);
b = Afun(x);

% Replaced Ainv with variable M (not using inverse)
M = A+rand(1000,1000).*1e-4;
precond = @(x)(M*x);

[Y,flag,relres,piter,resvec] = gmres(Afun, b, 32,1e-12,32);

fprintf('  reported relres without precond: %10.4e\n',relres)

% manually compute relres
manual_relres = norm(b - Afun(Y))/norm(b);
fprintf('  manual relres: %10.4e\n',manual_relres)

[Y,flag,relres,piter,resvec] = gmres(Afun, b, 32,1e-12,32, precond);

fprintf('  reported relres with precond: %10.4e\n',relres)

% manually compute preconditioned relres
manual_precond_relres = norm(precond(b - Afun(Y)))/norm(precond(b));
fprintf('  manual precond relres: %10.4e\n',manual_precond_relres)

Output:

[响]

Thanks for the help again, but your preconditioned GMRES is not converging to the required accuracy (1e-12) in this case. So you might be not able to reproduce my result (inconsistency between reported and manual)

My example provides a good preconditioner that should be able to converge in very few steps.