Generalized minimal residual method

 ( Numerical Linear Algebra ) 

In mathematics, the generalized minimal residual method (GMRES) is an iterative method for the numerical solution of an indefinite nonsymmetric system of linear equations. The method approximates the solution by the vector in a Krylov subspace with minimal residual. The Arnoldi iteration is used to find this vector.

The GMRES method was developed by Yousef Saad and Martin H. Schultz in 1986. It is a generalization and improvement of the MINRES method due to Paige and Saunders in 1975.Paige and Saunders, "Solution of Sparse Indefinite Systems of Linear Equations" /ref> The MINRES method requires that the matrix is symmetric, but has the advantage that it only requires handling of three vectors. GMRES is a special case of the DIIS method developed by Peter Pulay in 1980. DIIS is applicable to non-linear systems.

The n-th Krylov sequence for this problem is $$K\_n\; =\; K\_n(A,r\_0)\; =\; \backslash operatorname\{span\}\; \backslash ,\; \backslash \{\; r\_0,\; Ar\_0,\; A^2r\_0,\; \backslash ldots,\; A^\{n-1\}r\_0\; \backslash \}.\; \backslash ,$$ where $r\_0\; =\; b\; -\; A\; x\_0$ is the initial residual given an initial guess $x\_0\; \backslash ne\; 0$. Clearly $r\_0\; =\; b$ if $x\_0\; =\; 0$.

GMRES approximates the exact solution of $Ax\; =\; b$ by the vector $x\_n\; \backslash in\; x\_0\; +\; K\_n$ that minimizes the Euclidean norm of the residual $r\_n=\; b\; -\; Ax\_n$.

The vectors $r\_0,Ar\_0,\backslash ldots\; A^\{n-1\}r\_0$ might be close to linearly dependent, so instead of this basis, the Arnoldi iteration is used to find orthonormal vectors $q\_1,\; q\_2,\; \backslash ldots,\; q\_n\; \backslash ,$ which form a basis for $K\_n$. In particular, $q\_1=\backslash |r\_0\backslash |\_2^\{-1\}r\_0$.

Therefore, the vector $x\_n\; \backslash in\; x\_0\; +\; K\_n$ can be written as $x\_n\; =\; x\_0\; +\; Q\_n\; y\_n$ with $y\_n\; \backslash in\; \backslash mathbb\{R\}^n$, where $Q\_n$ is the m-by- n matrix formed by $q\_1,\backslash ldots,q\_n$. In other words, finding the n-th approximation of the solution (i.e., $x\_n$) is reduced to finding the vector $y\_n$, which is determined via minimizing the residue as described below.

The Arnoldi process also constructs $\backslash tilde\{H\}\_n$, an ($n+1$)-by-$n$ upper Hessenberg matrix which satisfies $$AQ\_n\; =\; Q\_\{n+1\}\; \backslash tilde\{H\}\_n\; \backslash ,$$ an equality which is used to simplify the calculation of $y\_n$ (see ). Note that, for symmetric matrices, a symmetric tri-diagonal matrix is actually achieved, resulting in the MINRES method.

Because columns of $Q\_n$ are orthonormal, we have where $$e\_1\; =\; (1,0,0,\backslash ldots,0)^T\; \backslash ,$$ is the first vector in the standard basis of $\backslash mathbb\{R\}^\{n+1\}$, and$$\backslash beta\; =\; \backslash |r\_0\backslash |\; \backslash ,\; ,$$ $r\_0$ being the first trial residual vector (usually $b$). Hence, $x\_n$ can be found by minimizing the Euclidean norm of the residual$$r\_n\; =\; \backslash tilde\{H\}\_n\; y\_n\; -\; \backslash beta\; e\_1.$$ This is a linear least squares problem of size n.

This yields the GMRES method. On the $n$-th iteration:

1. calculate $q\_n$ with the Arnoldi method;
2. find the $y\_n$ which minimizes $\backslash |r\_n\backslash |$;
3. compute $x\_n\; =\; x\_0\; +\; Q\_n\; y\_n$;
4. repeat if the residual is not yet small enough.

This does not happen in general. Indeed, a theorem of Greenbaum, Pták and Strakoš states that for every nonincreasing sequence a₁, ..., a _m−1, a _m = 0, one can find a matrix A such that the = a _n for all n, where r _n is the residual defined above. In particular, it is possible to find a matrix for which the residual stays constant for m − 1 iterations, and only drops to zero at the last iteration.

In practice, though, GMRES often performs well. This can be proven in specific situations. If the symmetric part of A, that is $(A^T\; +\; A)/2$, is positive definite, then $$\backslash |r\_n\backslash |\; \backslash leq\; \backslash left(\; 1-\backslash frac\{\backslash lambda\_\{\backslash min\}^2(1/2(A^T\; +\; A))\}\{\; \backslash lambda\_\{\backslash max\}(A^T\; A)\}\; \backslash right)^\{n/2\}\; \backslash |r\_0\backslash |,$$ where $\backslash lambda\_\{\backslash mathrm\{min\}\}(M)$ and $\backslash lambda\_\{\backslash mathrm\{max\}\}(M)$ denote the smallest and largest eigenvalue of the matrix $M$, respectively.. NB all results for GCR also hold for GMRES, cf.

If A is symmetric matrix and positive definite, then we even have $$\backslash |r\_n\backslash |\; \backslash leq\; \backslash left(\; \backslash frac\{\backslash kappa\_2(A)^2-1\}\{\backslash kappa\_2(A)^2\}\; \backslash right)^\{n/2\}\; \backslash |r\_0\backslash |.$$ where $\backslash kappa\_2(A)$ denotes the condition number of A in the Euclidean norm.

In the general case, where A is not positive definite, we have $$\backslash frac\{\backslash |r\_n\backslash |\}\{\backslash |b\backslash |\}\; \backslash le\; \backslash inf\_\{p\; \backslash in\; P\_n\}\; \backslash |p(A)\backslash |\; \backslash le\; \backslash kappa\_2(V)\; \backslash inf\_\{p\; \backslash in\; P\_n\}\; \backslash max\_\{\backslash lambda\; \backslash in\; \backslash sigma(A)\}\; |p(\backslash lambda)|,\; \backslash ,$$ where P _n denotes the set of polynomials of degree at most n with p(0) = 1, V is the matrix appearing in the spectral decomposition of A, and σ( A) is the spectrum of A. Roughly speaking, this says that fast convergence occurs when the eigenvalues of A are clustered away from the origin and A is not too far from normal matrix.

All these inequalities bound only the residuals instead of the actual error, that is, the distance between the current iterate x _n and the exact solution.

The cost of the iterations grow as O( n²), where n is the iteration number. Therefore, the method is sometimes restarted after a number, say k, of iterations, with x _k as initial guess. The resulting method is called GMRES( k) or Restarted GMRES. For non-positive definite matrices, this method may suffer from stagnation in convergence as the restarted subspace is often close to the earlier subspace.

The shortcomings of GMRES and restarted GMRES are addressed by the recycling of Krylov subspace in the GCRO type methods such as GCROT and GCRODR. Recycling of Krylov subspaces in GMRES can also speed up convergence when sequences of linear systems need to be solved.

Another class of methods builds on the unsymmetric Lanczos iteration, in particular the BiCG method. These use a three-term recurrence relation, but they do not attain the minimum residual, and hence the residual does not decrease monotonically for these methods. Convergence is not even guaranteed.

The third class is formed by methods like CGS and BiCGSTAB. These also work with a three-term recurrence relation (hence, without optimality) and they can even terminate prematurely without achieving convergence. The idea behind these methods is to choose the generating polynomials of the iteration sequence suitably.

None of these three classes is the best for all matrices; there are always examples in which one class outperforms the other. Therefore, multiple solvers are tried in practice to see which one is the best for a given problem.

The minimum can be computed using a QR decomposition: find an ( n + 1)-by-( n + 1) orthogonal matrix Ω _n and an ( n + 1)-by- n upper triangular matrix $\backslash tilde\{R\}\_n$ such that $$\backslash Omega\_n\; \backslash tilde\{H\}\_n\; =\; \backslash tilde\{R\}\_n.$$ The triangular matrix has one more row than it has columns, so its bottom row consists of zero. Hence, it can be decomposed as $$\backslash tilde\{R\}\_n\; =\; \backslash begin\{bmatrix\}\; R\_n\; \backslash \backslash \; 0\; \backslash end\{bmatrix\},$$ where $R\_n$ is an n-by- n (thus square) triangular matrix.

The QR decomposition can be updated cheaply from one iteration to the next, because the Hessenberg matrices differ only by a row of zeros and a column: where h _n+1 = ( h_{1, n+1}, ..., h _{n+1, n+1})^T. This implies that premultiplying the Hessenberg matrix with Ω _n, augmented with zeroes and a row with multiplicative identity, yields almost a triangular matrix: This would be triangular if σ is zero. To remedy this, one needs the Givens rotation where $$c\_n\; =\; \backslash frac\{\backslash rho\}\{\backslash sqrt\{\backslash rho^2+\backslash sigma^2\}\}\; \backslash quad\backslash text\{and\}\backslash quad\; s\_n\; =\; \backslash frac\{\backslash sigma\}\{\backslash sqrt\{\backslash rho^2+\backslash sigma^2\}\}.$$ With this Givens rotation, we form Indeed, is a triangular matrix with $r\_\{n+1,n+1\}\; =\; \backslash sqrt\{\backslash rho^2+\backslash sigma^2\}$.

Given the QR decomposition, the minimization problem is easily solved by noting that Denoting the vector $\backslash beta\backslash Omega\_ne\_1$ by $$\backslash tilde\{g\}\_n\; =\; \backslash begin\{bmatrix\}\; g\_n\; \backslash \backslash \; \backslash gamma\_n\; \backslash end\{bmatrix\}$$ with g _n ∈ R ⁿ and γ _n ∈ R, this is The vector y that minimizes this expression is given by $$y\_n\; =\; R\_n^\{-1\}\; g\_n.$$ Again, the vectors $g\_n$ are easy to update.

 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;
 % Note: this is not the beta scalar in section "The method" above but
 % the beta scalar multiplied by e1
 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
 % 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;
     

 q = A*Q(:,k);   % Krylov Vector
 for i = 1:k     % Modified Gram-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);
     

 % 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;
     

   t = sqrt(v1^2 + v2^2);
     

   cs = v1 / t;  % see http://www.netlib.org/eispack/comqr.f
   sn = v2 / t;
     

• Biconjugate gradient method

• (2025). 9783528131357, Vieweg. ISBN 9783528131357
• Y. Saad (2025). 9780898715347, SIAM. ISBN 9780898715347
• Dongarra et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition, SIAM, Philadelphia, 1994