Generalized minimal residual method (GRMRES)

The generalized minimal residual method (GMRES) is an iterative method for the numerical solution of a nonsymmetric system of linear equations. The method approximates the solution by the vector in a Krylov subspace with minimal residual. The algorithm builds an orthogonal basis of the Krylov subspaces and then solves the least square problem to find this vector. To built an orthogonal basis, one may use the Arnoldi iterations. After the orthogonalization, we solve the least square by the QR decomposition, where the Given rotation is used.

At every iteration, a matrix-vector product A q_n must be computed. This costs about 2m^2 floating-point operations for general dense matrices of size m, but the cost can decrease to O(m) for sparse matrices. In addition to the matrix-vector product, O(nm) floating-point operations must be computed at the n-th iteration.

Convergence

The n-th iterate minimizes the residual in the Krylov subspace K^n. Since every subspace is contained in the next subspace, the residual does not increase. After m iterations, where m is the size of the matrix A, the Krylov space K^m is the whole of R^m and hence the GMRES method arrives at the exact solution. However, the idea is that after a small number of iterations (relative to m), the vector x_n is already a good approximation to the exact solution.

A pseudocode of the GMRES algorithm (in MATLAB) is listed bellow.

    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(1) = 1;
      e=[error];
      r_norm=norm(r);
      Q(:,1) = r/r_norm;
      beta = r_norm*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
      
      %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);
  for i = 1:k
    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;
  end
end

Related program(s):

BEM4I miniApp

CalculiX solver