How are preconditioned GMRES relative residual computed?

37 次查看（过去 30 天）
响 李
响 李 2024-10-14
编辑： Steven Lord 2024-10-18
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

请先登录，再回答此问题。

采纳的回答

Heiko Weichelt
Heiko Weichelt 2024-10-18
编辑：Steven Lord 2024-10-18
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 ||b-A*x||/||M\b||, 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]

更多回答（1 个）

Abhinav Aravindan
Abhinav Aravindan 2024-10-14
编辑：Abhinav Aravindan 2024-10-14
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:
I hope this helps resolve your query!
  3 个评论
显示 1更早的评论
Abhinav Aravindan
Abhinav Aravindan 2024-10-14
编辑：Abhinav Aravindan 2024-10-14
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:
响 李
响 李 2024-10-14
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.

请先登录，再进行评论。

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by