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

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

n = 400; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = sum(A,2);

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

x = minres(A,b);

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

By default minres uses 20 iterations and a tolerance of 1e-6, and the algorithm is unable to converge in those 20 iterations for this matrix. Since the residual is on the order of 1e-2, it is a good indicator that more iterations 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 250 iterations.

x = minres(A,b,1e-4,250);

minres converged at iteration 200 to a solution with relative residual 7e-13.

Examine the effect of using a preconditioner matrix with minres to solve a linear system.

Create a symmetric positive definite, banded coefficient matrix.

A = delsq(numgrid('S',102));

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 = 100;

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

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

fl0 = 1

rr0

rr0 = 0.0013

it0

it0 = 100

fl0 is 1 because minres does not converge to the requested tolerance 1e-12 within the requested 100 iterations.

To aid with convergence, you can specify a preconditioner matrix. Since A is symmetric, use ichol to generate the preconditioner $\mathit{M}=\mathit{L}{\mathit{L}}^{\mathit{T}}$. Specify the 'ict' option to use incomplete Cholesky factorization with threshold dropping, and specify a diagonal shift value of 1e-6 to avoid nonpositive pivots. Solve the preconditioned system by specifying L and L' as inputs to minres.

setup = struct('type','ict','diagcomp',1e-6,'droptol',1e-14);
L = ichol(A,setup);
[x1,fl1,rr1,it1,rv1,rvcg1] = minres(A,b,tol,maxit,L,L');
fl1

fl1 = 0

rr1

rr1 = 2.6691e-15

it1

it1 = 4

The use of an ilu preconditioner produces a relative residual less than the prescribed tolerance of 1e-12 at the fourth iteration. The output rv1(1) is norm(b), and the output rv1(end) is norm(b-A*x1).

You can follow the progress of minres 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(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ICHOL preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Examine the effect of supplying minres 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 minres 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 = minres(A,b,[],maxit);

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

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

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

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

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

Solve a linear system by providing minres 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 minres with the function handle that calculates A*x. Use a tolerance of 1e-12 and 50 iterations.

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = minres(@afun,b,tol,maxit)

minres converged at iteration 11 to a solution with relative residual 4.1e-16.

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