qmr

Solve system of linear equations — quasi-minimal residual method

Description

example

x = qmr(A,b) attempts to solve the system of linear equations A*x = b for x using the Quasi-minimal Residual Method. When the attempt is successful, qmr displays a message to confirm convergence. If qmr fails to converge after the maximum number of iterations or halts for any reason, it displays a diagnostic message that includes the relative residual norm(b-A*x)/norm(b) and the iteration number at which the method stopped.

example

x = qmr(A,b,tol) specifies a tolerance for the method. The default tolerance is 1e-6.

example

x = qmr(A,b,tol,maxit) specifies the maximum number of iterations to use. qmr displays a diagnostic message if it fails to converge within maxit iterations.

example

x = qmr(A,b,tol,maxit,M) specifies a preconditioner matrix M and computes x by effectively solving the system M1Ax=M1b for x. Using a preconditioner matrix can improve the numerical properties of the problem and the efficiency of the calculation.

example

x = qmr(A,b,tol,maxit,M1,M2) specifies factors of the preconditioner matrix M such that M = M1*M2.

example

x = qmr(A,b,tol,maxit,M1,M2,x0) specifies an initial guess for the solution vector x. The default is a vector of zeros.

example

[x,flag] = qmr(___) returns a flag that specifies whether the algorithm successfully converged. When flag = 0, convergence was successful. You can use this output syntax with any of the previous input argument combinations. When you specify the flag output, qmr does not display any diagnostic messages.

example

[x,flag,relres] = qmr(___) also returns the relative residual norm(b-A*x)/norm(b). If flag is 0, then relres <= tol.

example

[x,flag,relres,iter] = qmr(___) also returns the iteration number iter at which x was computed.

example

[x,flag,relres,iter,resvec] = qmr(___) also returns a vector of the residual norm at each iteration, including the first residual norm(b-A*x0).

Examples

collapse all

Solve a square linear system using qmr 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. Also create a random vector b for the right-hand side of Ax=b.

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

Solve Ax=b using qmr. The output display includes the value of the relative residual error Ax-bb.

x = qmr(A,b);
qmr 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.12.

By default qmr uses 20 iterations and a tolerance of 1e-6, and the algorithm is unable to converge in those 40 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 reduce the tolerance to make it easier for the algorithm to converge.

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

x = qmr(A,b,1e-4,100);
qmr 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.063.

Even with a looser tolerance and more iterations the residual error does not improve much. When an iterative algorithm stalls in this manner it is a good indication that a preconditioner matrix is needed.

Calculate the incomplete Cholesky factorization of A, and use the L' factor as a preconditioner input to qmr.

L = ichol(A);
x = qmr(A,b,1e-4,100,L');
qmr converged at iteration 57 to a solution with relative residual 6.1e-05.

Using a preconditioner improves the numerical properties of the problem enough that qmr is able to converge.

Examine the effect of using a preconditioner matrix with qmr 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 is a vector of all ones.

b = sum(A,2);

Set the tolerance and maximum number of iterations.

tol = 1e-12;
maxit = 20;

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

  • x0 is the computed solution to A*x0 = b.

  • fl0 is a flag indicating whether the algorithm converged.

  • rr0 is the residual of the computed answer x0.

  • it0 is the iteration number when x0 was computed.

  • rv0 is a vector of the residual history for Ax-b.

[x0,fl0,rr0,it0,rv0] = qmr(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.7984
it0
it0 = 17

fl0 is 1 because qmr does not converge to the requested tolerance 1e-12 within the requested 20 iterations. The 17th iterate is the best approximate solution and is the one returned as indicated by it0 = 17.

To aid with the slow convergence, you can specify a preconditioner matrix. Since A is nonsymmetric, use ilu to generate the preconditioner M=LU. Specify a drop tolerance to ignore nondiagonal entries with values smaller than 1e-6. Solve the preconditioned system M-1Ax=M-1b by specifying L and U as inputs to qmr.

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

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

You can follow the progress of qmr 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','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Examine the effect of supplying qmr 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 Ax=b so that the expected solution for 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 qmr to solve Ax=b twice: one time with the default initial guess, and one time with a good initial guess of the solution. Use 200 iterations for both solutions and specify the initial guess as a vector with all elements equal to 0.99.

maxit = 200;
x1 = qmr(A,b,[],maxit);
qmr converged at iteration 27 to a solution with relative residual 9.5e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = qmr(A,b,[],maxit,[],[],x0);
qmr converged at iteration 7 to a solution with relative residual 6.7e-07.

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

Returning Intermediate Results

You also can use the initial guess to get intermediate results by calling qmr in a for-loop. Each call to the solver performs a few iterations and stores the calculated solution. Then you use that solution as the initial vector for the next batch of iterations.

For example, this code performs 100 iterations four times and stores the solution vector after each pass in the for-loop:

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

X(:,k) is the solution vector computed at iteration k of the for-loop, and R(k) is the relative residual of that solution.

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

Create an nonsymmetric tridiagonal matrix. Preview the matrix.

A = gallery('wilk',21) + diag(ones(20,1),1)
A = 21×21

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

Since this tridiagonal matrix has a special structure, you can represent the operation A*x with a function handle. As each row of A multiplies the elements in x, only a few of the results are nonzero (corresponding to the nonzeros on the tridiagonals).

The expression Ax becomes:

Ax=[1020019200120010010200110][x1x2x3x21]=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21].

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

Ax=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21]=[0x1x2x20]+[10x19x29x2010x21]+2[x2x3x210].

Likewise, the expression for ATx becomes:

ATx=[1010029100210020010100210][x1x2x3x21]=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21].

ATx=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21]=2[0x1x2x20]+[10x19x29x2010x21]+[x2x3x210].

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

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

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

Now, solve the linear system Ax=b by providing qmr with the function handle that calculates A*x and A'*x. Use a tolerance of 1e-6 and 25 iterations. Specify b as the row sums of A so that the true solution for x is a vector of ones.

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = qmr(@afun,b,tol,maxit)
qmr converged at iteration 19 to a solution with relative residual 4.7e-07.
x1 = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

Local Functions

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

Input Arguments

collapse all

Coefficient matrix, specified as a square matrix or function handle. This matrix is the coefficient matrix in the linear system A*x = b. Generally, A is a large sparse matrix or a function handle that returns the product of a large sparse matrix and column vector.

Specifying A as a Function Handle

You can specify the coefficient matrix as a function handle instead of a matrix to save memory in the calculation. The function handle returns matrix-vector products instead of forming the entire coefficient matrix, making the calculation more efficient.

To use a function handle, use the function signature function y = afun(x,opt). Parameterizing Functions explains how to provide additional parameters to the function afun, if necessary. The function afun must satisfy these conditions:

  • afun(x,'notransp') returns the product A*x.

  • afun(x,'transp') returns the product A'*x.

An example of an acceptable function is:

function y = afun(x,opt,B,C,n)
if strcmp(opt,'notransp')
    y = [B*x(n+1:end); C*x(1:n)];
else
    y = [C'*x(n+1:end); B'*x(1:n)];
end
The function afun uses B and C to compute either A*x or A'*x (depending on the specified flag) without actually forming the entire sparse matrix A = [zeros(n) B; C zeros(n)]. This exploits the sparsity pattern of the matrix to save memory in the computation of A*x and A'*x.

Data Types: double | function_handle
Complex Number Support: Yes

Right-hand side of linear equation, specified as a column vector. b must be a column vector with length equal to size(A,1).

Data Types: double
Complex Number Support: Yes

Method tolerance, specified as a positive scalar. Use this input to trade-off accuracy and runtime in the calculation. qmr must meet the tolerance within the number of allowed iterations to be successful. A smaller value of tol means the answer must be more precise for the calculation to be successful.

Data Types: double

Maximum number of iterations, specified as a positive scalar integer. Increase the value of maxit to allow more iterations for qmr to meet the tolerance tol. Generally, a smaller value of tol means more iterations are required to successfully complete the calculation.

Preconditioner matrices, specified as separate arguments of matrices or function handles. You can specify a preconditioner matrix M or its matrix factors M = M1*M2 to improve the numerical aspects of the linear system and make it easier for qmr to converge quickly. You can use the incomplete matrix factorization functions ilu and ichol to generate preconditioner matrices. You also can use equilibrate prior to factorization to improve the condition number of the coefficient matrix. For more information on preconditioners, see Iterative Methods for Linear Systems.

qmr treats unspecified preconditioners as identity matrices.

Specifying M as a Function Handle

You can specify any of M, M1, or M2 as function handles instead of matrices to save memory in the calculation. The function handle performs matrix-vector operations instead of forming the entire preconditioner matrix, making the calculation more efficient.

To use a function handle, first create a function with the signature function y = mfun(x,opt). Parameterizing Functions explains how to provide additional parameters to the function mfun, if necessary. The function mfun must satisfy these conditions:

  • mfun(x,'notransp') returns the value of M\x or M2\(M1\x).

  • mfun(x,'transp') returns the value of M'\x or M1'\(M2'\x).

An example of an acceptable function is:

function y = mfun(x,opt,a,b)  
if strcmp(opt,'notransp')
    y = x.*a;
else
    y = x.*b;
end
end
In this example the function mfun uses a and b to compute either M\x = x*a or M'\x = x*b (depending on the specified flag) without actually forming the entire sparse matrix M.

Data Types: double | function_handle
Complex Number Support: Yes

Initial guess, specified as a column vector with length equal to size(A,2). If you can provide qmr with a more reasonable initial guess x0 than the default vector of zeros, then it can save computation time and help the algorithm converge faster.

Data Types: double
Complex Number Support: Yes

Output Arguments

collapse all

Linear system solution, returned as a vector. This output gives the approximate solution to the linear system A*x = b. If the calculation is successful (flag = 0), then relres is less than or equal to tol.

Whenever the calculation is not successful (flag ~= 0), the solution x returned by qmr is the one with minimal residual norm computed over all the iterations.

Convergence flag, returned as one of the scalar values in this table. The convergence flag indicates whether the calculation was successful and differentiates between several different forms of failure.

Flag Value

Convergence

0

Success — qmr converged to the desired tolerance tol within maxit iterations.

1

Failure — qmr iterated maxit iterations but did not converge.

2

Failure — The preconditioner matrix M or M = M1*M2 is ill conditioned.

3

Failure — qmr stagnated after two consecutive iterations were the same.

4

Failure — One of the scalar quantities calculated by the qmr algorithm became too small or too large to continue computing.

Relative residual error, returned as a scalar. The relative residual error relres = norm(b-A*x)/norm(b) is an indication of how accurate the answer is. If the calculation converges to the tolerance tol within maxit iterations, then relres <= tol.

Data Types: double

Iteration number, returned as a scalar. This output indicates the iteration number at which the computed answer for x was calculated.

Data Types: double

Residual error, returned as a vector. The residual error norm(b-A*x) reveals how close the algorithm is to converging for a given value of x. The number of elements in resvec is equal to the number of iterations. You can examine the contents of resvec to help decide whether to change the values of tol or maxit.

Data Types: double

More About

collapse all

Quasi-minimal Residual Method

The QMR algorithm was developed as an improvement to BiCG. While GMRES uses an orthogonal basis for the Krylov subspace and computes a minimum residual solution, QMR uses a bi-orthogonal basis and therefore computes only a quasi-minimal residual solution.

QMR typically converges much more smoothly than BiCG, and it also uses look-ahead techniques to avoid breakdowns in almost all cases. The computational cost of QMR is only slightly higher than for BiCG [1].

Tips

  • Convergence of most iterative methods depends on the condition number of the coefficient matrix, cond(A). You can use equilibrate to improve the condition number of A, and on its own this makes it easier for most iterative solvers to converge. However, using equilibrate also leads to better quality preconditioner matrices when you subsequently factor the equilibrated matrix B = R*P*A*C.

  • You can use matrix reordering functions such as dissect and symrcm to permute the rows and columns of the coefficient matrix and minimize the number of nonzeros when the coefficient matrix is factored to generate a preconditioner. This can reduce the memory and time required to subsequently solve the preconditioned linear system.

References

[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

[2] Freund, Roland W. and Nöel M. Nachtigal, “QMR: A quasi-minimal residual method for non-Hermitian linear systems,” SIAM Journal: Numer. Math. 60, 1991, pp. 315–339.

Extended Capabilities

Introduced before R2006a