Residual in GMRES

[MichelKern]

I am confused by the computation of the GMRES residual, and the stopping test:

The variable e_all_iter is initialized as r_norm which is the actual norm of the initial residual (line 1273)

Then, during the iterations, one has (lines 1461 - 1465)

              combined_error_iter=fabs(Beta[k_counter+1])/r_norm;
              e_all_iter[k_counter+1]=(combined_error_iter);

so that e_all_iter is the decrease in the residual (but then the first value is wrong ?), and this is the value that gets printed.

Finally, the stopping test, and the value reported (lines 1471 - 1473) is combined_error_iter/e_all_iter[0], so there is one extra division by the initial residual.

I think what should be done is to set e_all_iter to the actual value of the residual, and when it is needed to divide by e_all_iter[0] (for deciding where to stop). Maybe also report the "target" residual when the iterations are started. Something like: "initial residual is " << r_norm << " and stopping criteria is residual < " << tolerance * r_norm

[mjayadharan]

I don't remember the exact variables and might need to take a close look at them. Let me get back to you on this. Manu.

On Tue, Sep 8, 2020, 8:07 AM Michel Kern notifications@github.com wrote:

I am confused by the computation of the GMRES residual, and the stopping test:

The variable e_all_iter is initialized as r_norm which is the actual norm of the initial residual (line 1273)

Then, during the iterations, one has (lines 1461 - 1465)

          combined_error_iter=fabs(Beta[k_counter+1])/r_norm;
          e_all_iter[k_counter+1]=(combined_error_iter);

so that e_all_iter is the decrease in the residual (but then the first value is wrong ?), and this is the value that gets printed.

Finally, the stopping test, and the value reported (lines 1471 - 1473) is combined_error_iter/e_all_iter[0], so there is one extra division by the initial residual.

I think what should be done is to set e_all_iter to the actual value of the residual, and when it is needed to divide by e_all_iter[0] (for deciding where to stop). Maybe also report the "target" residual when the iterations are started. Something like: "initial residual is " << r_norm << " and stopping criteria is residual < " << tolerance * r_norm

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/mjayadharan/MMMFE-ST-DD/issues/29, or unsubscribe https://github.com/notifications/unsubscribe-auth/AIR5RWLW54OTEXKUW4KRDVDSEYNBRANCNFSM4Q74IHHA .

[mjayadharan]

The local gmres code is a matrix free adaptation of the following matlab implementation which can be found here.

function [x, e] = gmres( A, b, x, max_iterations, threshold)
  n = length(A);
  m = max_iterations;

  % use x as the initial vector
  r = b - A * x;

  b_norm = norm(b);
  error = norm(r) / b_norm;

  % initialize the 1D vectors
  sn = zeros(m, 1);
  cs = zeros(m, 1);
  %e1 = zeros(n, 1);
  e1 = zeros(m+1, 1);
  e1(1) = 1;
  e = [error];
  r_norm = norm(r);
  Q(:,1) = r / r_norm;
  beta = r_norm * e1;     %Note: this is not the beta scalar in section "The method" above but the beta scalar multiplied by e1
  for k = 1:m

    % run arnoldi
    [H(1:k+1, k) Q(:, k+1)] = arnoldi(A, Q, k);

    % eliminate the last element in H ith row and update the rotation matrix
    [H(1:k+1, k) cs(k) sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k);

    % update the residual vector
    beta(k + 1) = -sn(k) * beta(k);
    beta(k)     = cs(k) * beta(k);
    error       = abs(beta(k + 1)) / b_norm;

    % save the error
    e = [e; error];

    if (error <= threshold)
      break;
    end
  end
  % if threshold is not reached, k = m at this point (and not m+1) 

  % calculate the result
  y = H(1:k, 1:k) \ beta(1:k);
  x = x + Q(:, 1:k) * y;
end

%----------------------------------------------------%
%                  Arnoldi Function                  %
%----------------------------------------------------%
function [h, q] = arnoldi(A, Q, k)
  q = A*Q(:,k);   % Krylov Vector
  for i = 1:k     % Modified Gramm-Schmidt, keeping the Hessenberg matrix
    h(i) = q' * Q(:, i);
    q = q - h(i) * Q(:, i);
  end
  h(k + 1) = norm(q);
  q = q / h(k + 1);
end

%---------------------------------------------------------------------%
%                  Applying Givens Rotation to H col                  %
%---------------------------------------------------------------------%
function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k)
  % apply for ith column
  for i = 1:k-1
    temp   =  cs(i) * h(i) + sn(i) * h(i + 1);
    h(i+1) = -sn(i) * h(i) + cs(i) * h(i + 1);
    h(i)   = temp;
  end

  % update the next sin cos values for rotation
  [cs_k sn_k] = givens_rotation(h(k), h(k + 1));

  % eliminate H(i + 1, i)
  h(k) = cs_k * h(k) + sn_k * h(k + 1);
  h(k + 1) = 0.0;
end

%%----Calculate the Given rotation matrix----%%
function [cs, sn] = givens_rotation(v1, v2)
%  if (v1 == 0)
%    cs = 0;
%    sn = 1;
%  else
    t = sqrt(v1^2 + v2^2);
%    cs = abs(v1) / t;
%    sn = cs * v2 / v1;
    cs = v1 / t;  % see http://www.netlib.org/eispack/comqr.f
    sn = v2 / t;
%  end
end

• It could be confusing that in _localgmres, there is division by r_norm instead of b_norm, compared to the snippet above. This is because our initial guess x_0 is zero and hence the residual norm, r_norm=norm(Ax_0-b), is same as b_norm.
• Regarding what is being printed in the output, you are welcome to change it and send me a pull request, I will merge the changes to the master. Let me know if you have any other questions before closing this issue.

[mjayadharan]

Also, if there is a possibility that division of error by r_norm will be confusing some time in the future, we could simply add a variable b_norm which will hold the place for error in first iteration. Of course this will be redundant, but will increase the readability of the algorithm.

[MichelKern]

It's fairly clear that r_norm is actually norm(b), though you're right that it could get confusing later on. You're probably better off by erring on the side of clarity vs redundancy.

But the main issue is not the printing, but the test: in the matlab code the test is

if (error <= threshold)

whereas the C++ code has

  if (combined_error_iter/e_all_iter[0] < tolerance)

but the division by e_all_iter[0] has already taken place (e_all_iter is the same thing as error, and both actually contain residual reduction )

I might still be mistaken, so I want to make sure we both agree before I make any changes.

(I also removed a comment above that was for the wrong issue)

[mjayadharan]

I now see what you are pointing to. Sorry I overlooked it the first time. There indeed are two divisions by b_norm. This seems be a bug, which might affect the convergence and number of iterations. Please go ahead and make the changes. Thanks.

[mjayadharan]

This issue has been resolved, using a pull request from @MichelKern .

[MichelKern]

Thanks for the merge. I still prefer to give you an opportunity (or to force you ?) to review the changes before merging into master.

----- Le 9 Sep 20, à 14:57, Manu Jayadharan notifications@github.com a écrit :

I now see what you are pointing to. Sorry I overlooked it the first time. There indeed are two divisions by b_norm . This seems be a bug, which might affect the convergence and number of iterations. Please go ahead and make the changes. Thanks.

— You are receiving this because you were assigned. Reply to this email directly, [ https://github.com/mjayadharan/MMMFE-ST-DD/issues/29#issuecomment-689545116 | view it on GitHub ] , or [ https://github.com/notifications/unsubscribe-auth/AA4OOBKVWHSTIW375NPOKWTSE53T5ANCNFSM4Q74IHHA | unsubscribe ] .

-- Michel Kern Centre de recherche Inria de Paris, 2 rue Simone Iff, bureau A413, 75012 Paris tel. : +33 (0)1 80 49 42 36,

[mjayadharan]

I have looked at the changes carefully before merging the request. I wasn't sure I was supposed to open a review formally before merging. May be I will do that from next time. Thanks for pointing out.

[MichelKern]

No, no I meant review in a loose sense, just to make sure you agree with the changes. A formal review is probably overkillEnvoyé depuis mon appareil mobile Samsung. -------- Message d'origine --------De : Manu Jayadharan notifications@github.com Date : 12/09/2020 20:52 (GMT+01:00) À : mjayadharan/MMMFE-ST-DD MMMFE-ST-DD@noreply.github.com Cc : Michel Kern michel.kern@inria.fr, Mention mention@noreply.github.com Objet : Re: [mjayadharan/MMMFE-ST-DD] Residual in GMRES  (#29) I have looked at the changes carefully before merging the request. I wasn't sure I was supposed to open a review formally before merging. May be I will do that from next time. Thanks for pointing out.

—You are receiving this because you were mentioned.Reply to this email directly, view it on GitHub, or unsubscribe.