How are preconditioned GMRES relative residual computed?
28 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Heiko Weichelt
on 18 Oct 2024
Edited: Steven Lord
on 18 Oct 2024
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]
0 Comments
More Answers (1)
Abhinav Aravindan
on 14 Oct 2024
Edited: Abhinav Aravindan
on 14 Oct 2024
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:
I hope this helps resolve your query!
3 Comments
Abhinav Aravindan
on 14 Oct 2024
Edited: Abhinav Aravindan
on 14 Oct 2024
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:
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!