gmres

Solve a square linear system using gmres with default settings, and then adjust the tolerance and number of iterations used in the solution process.

Create a random sparse matrix A with 50% density and nonzeros on the main diagonal. Also create a random vector b for the right-hand side of $\mathrm{Ax}=\mathit{b}$.

rng default
A = sprandn(400,400,0.5) + 12*speye(400);
b = rand(400,1);

Solve $\mathrm{Ax}=\mathit{b}$ using gmres. The output display includes the value of the relative residual error $\frac{\Vert \mathit{b}-\mathrm{Ax}\Vert }{\Vert \mathit{b}\Vert }$.

x = gmres(A,b);

gmres stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.35.

By default gmres uses 10 iterations and a tolerance of 1e-6, and the algorithm is unable to converge in those 10 iterations for this matrix. Since the residual is still large, it is a good indicator that more iterations (or a preconditioner matrix) are needed. You also can use a larger tolerance to make it easier for the algorithm to converge.

Solve the system again using a tolerance of 1e-4 and 100 iterations.

tol = 1e-4;
maxit = 100;
x = gmres(A,b,[],tol,maxit);

gmres stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.0045.

Even with a looser tolerance and more iterations the residual error does not improve enough for convergence. When an iterative algorithm stalls in this manner it is a good indication that a preconditioner matrix is needed. However, gmres also has an input that controls the number of inner iterations. By specifying a value for the inner iterations, gmres does more work per outer iteration.

Solve the system again using a restart value of 100 and a maxit value of 20. Rather than doing 100 iterations once, gmres performs 100 iterations between restarts and repeats this 20 times.

restart = 100;
maxit = 20;
x = gmres(A,b,restart,tol,maxit);

gmres(100) converged at outer iteration 2 (inner iteration 75) to a solution with relative residual 9.3e-05.

In this case specifying a large restart value for gmres enables it to converge to a solution within the allowed number of iterations. However, large restart values can consume a lot of memory when A is also large.

Examine the effect of using a preconditioner matrix with non-restarted gmres to solve a linear system.

Load west0479, a real 479-by-479 nonsymmetric sparse matrix.

load west0479
A = west0479;

Define b so that the true solution to $\mathrm{Ax}=\mathit{b}$ is a vector of all ones.

b = sum(A,2);

Set the tolerance and maximum number of iterations.

tol = 1e-12;
maxit = 20;

Use gmres to find a solution at the requested tolerance and number of iterations. Specify five outputs to return information about the solution process:

[x,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0

fl0 = 
1

rr0

rr0 = 
0.7603

it0

it0 = 1×2

     1    20

fl0 is 1 because gmres does not converge to the requested tolerance 1e-12 within the requested 20 iterations. The best approximate solution that gmres returns is the last one (as indicated by it0(2) = 20). MATLAB® stores the residual history in rv0.

To aid with the slow convergence, you can specify a preconditioner matrix. Since A is nonsymmetric, use ilu to generate the preconditioner $\mathit{M}=\mathit{L}\mathit{U}$. Specify a drop tolerance to ignore nondiagonal entries with values smaller than 1e-6. Solve the preconditioned system ${\mathit{M}}^{-1}\mathit{A}\mathit{x}={\mathit{M}}^{-1}\mathit{b}$ by specifying L and U as inputs to gmres. Note that when you specify a preconditioner, gmres calculates the residual norm of the preconditioned system for the outputs rr1 and rv1.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1

fl1 = 
0

rr1

rr1 = 
1.1870e-13

it1

it1 = 1×2

     1     6

The use of an ilu preconditioner produces a relative residual less than the prescribed tolerance of 1e-12 at the sixth iteration. The first residual rv1(1) is norm(U\(L\b)), where M = L*U. The last residual rv1(end) is norm(U\(L\(b-A*x1))).

You can follow the progress of gmres by plotting the relative residuals at each iteration. Plot the residual history of each solution with a line for the specified tolerance.

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Using a preconditioner with restarted gmres.

Load west0479, a real 479-by-479 nonsymmetric sparse matrix.

load west0479
A = west0479;

Define b so that the true solution to $\mathrm{Ax}=\mathit{b}$ is a vector of all ones.

b = sum(A,2);

Construct an incomplete LU preconditioner with a drop tolerance of 1e-6.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

The benefit to using restarted gmres is to limit the amount of memory required to execute the method. Without restart, gmres requires maxit vectors of storage to keep the basis of the Krylov subspace. Also, gmres must orthogonalize against all of the previous vectors at each step. Restarting limits the amount of workspace used and the amount of work done per outer iteration.

Execute gmres(3), gmres(4), and gmres(5) using the incomplete LU factors as preconditioners. Use a tolerance of 1e-12 and a maximum of 20 outer iterations.

tol = 1e-12;
maxit = 20;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,3,tol,maxit,L,U);
[x4,fl4,rr4,it4,rv4] = gmres(A,b,4,tol,maxit,L,U);
[x5,fl5,rr5,it5,rv5] = gmres(A,b,5,tol,maxit,L,U);
fl3

fl3 = 
0

fl4

fl4 = 
0

fl5

fl5 = 
0

fl3, fl4, and fl5 are all 0 because in each case restarted gmres drives the relative residual to less than the prescribed tolerance of 1e-12.

The following plot shows the convergence history of each restarted gmres method. gmres(3) converges at outer iteration 5, inner iteration 3 (it3 = [5, 3]) which would be the same as outer iteration 6, inner iteration 0, hence the marking of 6 on the final tick mark.

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

In general the larger the number of inner iterations, the more work gmres does per outer iteration and the faster it can converge.

Examine the effect of supplying gmres with an initial guess of the solution.

Create a tridiagonal sparse matrix. Use the sum of each row as the vector for the right-hand side of $\mathrm{Ax}=\mathit{b}$ so that the expected solution for $\mathit{x}$ is a vector of ones.

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

Use gmres to solve $\mathrm{Ax}=\mathit{b}$ twice: one time with the default initial guess, and one time with a good initial guess of the solution. Use 200 iterations and the default tolerance for both solutions. Specify the initial guess in the second solution as a vector with all elements equal to 0.99.

maxit = 200;
x1 = gmres(A,b,[],[],maxit);

gmres converged at iteration 27 to a solution with relative residual 9.5e-07.

x0 = 0.99*e;
x2 = gmres(A,b,[],[],maxit,[],[],x0);

gmres converged at iteration 7 to a solution with relative residual 6.7e-07.

In this case supplying an initial guess enables gmres to converge more quickly.

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = gmres(A,b,[],tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

Solve a linear system by providing gmres with a function handle that computes A*x in place of the coefficient matrix A.

One of the Wilkinson test matrices generated by gallery is a 21-by-21 tridiagonal matrix. Preview the matrix.

A = gallery('wilk',21)

A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

The Wilkinson matrix has a special structure, so you can represent the operation A*x with a function handle. When A multiplies a vector, most of the elements in the resulting vector are zeros. The nonzero elements in the result correspond with the nonzero tridiagonal elements of A. Moreover, only the main diagonal has nonzeros that are not equal to 1.

The expression $\mathrm{Ax}$ becomes:

$\mathrm{Ax}=\left[\begin{array}{cccccccccc}10& 1& 0& \cdots & & \cdots & & \cdots & 0& 0\\ 1& 9& 1& 0& & & & & & 0\\ 0& 1& 8& 1& 0& & & & & ⋮\\ ⋮& 0& 1& 7& 1& 0& & & & \\ & & 0& 1& 6& 1& 0& & & ⋮\\ ⋮& & & 0& 1& 5& 1& 0& & \\ & & & & 0& 1& 4& 1& 0& ⋮\\ ⋮& & & & & 0& 1& 3& \ddots & 0\\ 0& & & & & & 0& \ddots & \ddots & 1\\ 0& 0& \cdots & & \cdots & & \cdots & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ {\mathit{x}}_4\\ {\mathit{x}}_5\\ ⋮\\ \\ ⋮\\ \\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

The resulting vector can be written as the sum of three vectors:

$\mathrm{Ax}=\left[\begin{array}{c}0+10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}+0\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_1\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}10{\mathit{x}}_1\\ 9{\mathit{x}}_2\\ ⋮\\ 10{\mathit{x}}_{21}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

In MATLAB®, write a function that creates these vectors and adds them together, thus giving the value of A*x:

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

(This function is saved as a local function at the end of the example.)

Now, solve the linear system $\mathrm{Ax}=\mathit{b}$ by providing gmres with the function handle that calculates A*x. Use a tolerance of 1e-12, 15 outer iterations, and 10 inner iterations before restart.

b = ones(21,1);
tol = 1e-12;  
maxit = 15;
restart = 10;
x1 = gmres(@afun,b,restart,tol,maxit)

gmres(10) converged at outer iteration 5 (inner iteration 10) to a solution with relative residual 5.3e-13.

x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

Check that afun(x1) produces a vector of ones.

afun(x1)

ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end