pcg

Solve system of linear equations — preconditioned conjugate gradients method

Syntax

x = pcg(A,b)

x = pcg(A,b,tol)

x = pcg(A,b,tol,maxit)

x = pcg(A,b,tol,maxit,M)

x = pcg(A,b,tol,maxit,M1,M2)

x = pcg(A,b,tol,maxit,M1,M2,x0)

[x,flag] = pcg(___)

[x,flag,relres] = pcg(___)

[x,flag,relres,iter] = pcg(___)

[x,flag,relres,iter,resvec] = pcg(___)

Description

x = pcg(A,b) attempts to solve the system of linear equations A*x = b for x using the Preconditioned Conjugate Gradients Method. When the attempt is successful, pcg displays a message to confirm convergence. If pcg 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.

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

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

x = pcg(A,b,tol,maxit,M) specifies a preconditioner matrix M and computes x by effectively solving the system $$H^{-1}A{H}^{-T}y=H^{-1}b$$ for y, where $$y=H^{T}x$$ and $$H=M^{1/2}={\left(M_1{M}_2\right)}^{1/2}$$. The algorithm does not form H explicitly. Using a preconditioner matrix can improve the numerical properties of the problem and the efficiency of the calculation.

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

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

[x,flag] = pcg(___) 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, pcg does not display any diagnostic messages.

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

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

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

Examples

Iterative Solution to Linear System

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

Create a random symmetric sparse matrix A. Also create a vector b of the row sums of A for the right-hand side of $\mathrm{Ax}=\mathit{b}$ so that the true solution $\mathit{x}$ is a vector of ones.

rng default
A = sprand(400,400,.5);
A = A'*A;
b = sum(A,2);

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

x = pcg(A,b);

pcg 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 3.6e-06.

By default pcg uses 20 iterations and a tolerance of 1e-6, and the algorithm is unable to converge in those 20 iterations for this matrix. However, the residual is close to the tolerance, so the algorithm likely just needs more iterations to converge.

Solve the system again using a tolerance of 1e-7 and 150 iterations.

x = pcg(A,b,1e-7,150);

pcg converged at iteration 129 to a solution with relative residual 1e-07.

Using pcg with Preconditioner

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

Create a symmetric positive definite, banded coefficient matrix.

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

Define b for the right-hand side of the linear equation $\mathrm{Ax}=\mathit{b}$.

b = ones(size(A,1),1);

Set the tolerance and maximum number of iterations.

tol = 1e-8;
maxit = 100;

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

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

• fl0 is a flag indicating whether the algorithm converged.

• rr0 is the relative residual of the computed answer x.

• it0 is the iteration number when x was computed.

• rv0 is a vector of the residual history for $\Vert \mathit{b}-\mathrm{Ax}\Vert$.

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

fl0 = 
1

rr0

rr0 = 
0.0131

it0

it0 = 
100

fl0 is 1 because pcg does not converge to the requested tolerance of 1e-8 within the requested 100 iterations.

To aid with the slow 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}}$. Solve the preconditioned system by specifying L and L' as inputs to pcg.

L = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L,L');
fl1

fl1 = 
0

rr1

rr1 = 
8.0992e-09

it1

it1 = 
79

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

Now, use the michol option to create a modified incomplete Cholesky preconditioner.

L = ichol(A,struct('michol','on'));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L,L');
fl2

fl2 = 
0

rr2

rr2 = 
9.9614e-09

it2

it2 = 
47

This preconditioner is better than the one produced by the incomplete Cholesky factorization with zero fill for the coefficient matrix in this example, so pcg is able to converge even quicker.

You can see how the preconditioners affect the rate of convergence of pcg by plotting each of the residual histories starting from the initial estimate (iterate number 0). Add 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')
semilogy(0:length(rv2)-1,rv2/norm(b),'-o')
yline(tol,'r--');
legend('No Preconditioner','Default ICHOL','Modified ICHOL','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Supplying Initial Guess

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

pcg converged at iteration 35 to a solution with relative residual 9.5e-07.

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

pcg converged at iteration 7 to a solution with relative residual 8.7e-07.

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

Returning Intermediate Results

You also can use the initial guess to get intermediate results by calling pcg 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] = pcg(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.

Using Function Handle Instead of Numeric Matrix

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

Use gallery to generate a 20-by-20 positive definite tridiagonal matrix. The super- and subdiagonals have ones, while the main diagonal elements count down from 20 to 1. Preview the matrix.

n = 20;
A = gallery('tridiag',ones(n-1,1),n:-1:1,ones(n-1,1));
full(A)

ans = 20×20

    20     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1    19     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1    18     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1    17     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1    16     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1    15     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1    14     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1    13     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1    12     1     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1    11     1     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. 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:

$\left[\begin{array}{cccccccccc}20& 1& 0& \cdots & & \cdots & & \cdots & 0& 0\\ 1& 19& 1& 0& & & & & & 0\\ 0& 1& 18& 1& 0& & & & & ⋮\\ ⋮& 0& 1& 17& 1& 0& & & & \\ & & 0& 1& 16& 1& 0& & & ⋮\\ ⋮& & & 0& 1& 15& 1& 0& & \\ & & & & 0& 1& 14& 1& 0& ⋮\\ ⋮& & & & & 0& 1& 13& \ddots & 0\\ 0& & & & & & 0& \ddots & \ddots & 1\\ 0& 0& \cdots & & \cdots & & \cdots & 0& 1& 1\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ {\mathit{x}}_4\\ {\mathit{x}}_5\\ ⋮\\ \\ ⋮\\ \\ {\mathit{x}}_{20}\end{array}\right]=\left[\begin{array}{c}2{0\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+19{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+18{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{18}+2{\mathit{x}}_{19}+{\mathit{x}}_{20}\\ {\mathit{x}}_{19}+{\mathit{x}}_{20}\end{array}\right]$.

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

$\left[\begin{array}{c}2{0\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+19{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+18{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{18}+2{\mathit{x}}_{19}+{\mathit{x}}_{20}\\ {\mathit{x}}_{19}+{\mathit{x}}_{20}\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_1\\ ⋮\\ {\mathit{x}}_{19}\end{array}\right]+\left[\begin{array}{c}20{\mathit{x}}_1\\ 19{\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{20}\\ 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:19)] + ...
    [(20:-1:1)'].*x + ...
    [x(2:20); 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 pcg with the function handle that calculates A*x. Use a tolerance of 1e-12 and 50 iterations.

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

pcg converged at iteration 20 to a solution with relative residual 4.4e-16.

x1 = 20×1

    0.0476
    0.0475
    0.0500
    0.0526
    0.0555
    0.0588
    0.0625
    0.0666
    0.0714
    0.0769
      ⋮

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

afun(x1)

ans = 20×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)
y = [0; x(1:19)] + ...
    [(20:-1:1)'].*x + ...
    [x(2:20); 0];
end

Input Arguments

A — Coefficient matrix

matrix | function handle

Coefficient matrix, specified as a symmetric positive definite 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. See Determine Whether Matrix Is Symmetric Positive Definite for info on how to confirm that A is symmetric positive definite.

Specifying A as a Function Handle

You can optionally specify the coefficient matrix as a function handle instead of a matrix. 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). Parameterizing Functions explains how to provide additional parameters to the function afun, if necessary. The function call afun(x) must return the value of A*x.

Data Types: double | function_handle
Complex Number Support: Yes

b — Right-hand side of linear equation

column vector

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

Data Types: double
Complex Number Support: Yes

tol — Method tolerance

[] or 1e-6 (default) | positive scalar

Method tolerance, specified as a positive scalar. Use this input to trade-off accuracy and runtime in the calculation. pcg 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

maxit — Maximum number of iterations

[] or min(size(A,1),20) (default) | positive scalar integer

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

M, M1, M2 — Preconditioner matrices (as separate arguments)

eye(size(A)) (default) | matrices | function handles

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 pcg 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.

pcg treats unspecified preconditioners as identity matrices.

Specifying M as a Function Handle

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

To use a function handle, use the function signature function y = mfun(x). Parameterizing Functions explains how to provide additional parameters to the function mfun, if necessary. The function call mfun(x) must return the value of M\x or M2\(M1\x).

Data Types: double | function_handle
Complex Number Support: Yes

x0 — Initial guess

[] or a column vector of zeros (default) | column vector

Initial guess, specified as a column vector with length equal to size(A,2). If you can provide pcg 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

x — Linear system solution

column vector

Linear system solution, returned as a column 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 pcg is the one with minimal residual norm computed over all the iterations.

flag — Convergence flag

scalar

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 — pcg converged to the desired tolerance tol within maxit iterations.
1	Failure — pcg iterated maxit iterations but did not converge.
2	Failure — The preconditioner matrix M or M = M1*M2 is ill conditioned.
3	Failure — pcg stagnated after two consecutive iterations were the same.
4	Failure — One of the scalar quantities calculated by the pcg algorithm became too small or too large to continue computing.

relres — Relative residual error

scalar

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

iter — Iteration number

scalar

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

Data Types: double

resvec — Residual error

vector

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

Preconditioned Conjugate Gradients Method

The preconditioned conjugate gradients method (PCG) was developed to exploit the structure of symmetric positive definite matrices. Several other algorithms can operate on symmetric positive definite matrices, but PCG is the quickest and most reliable at solving those types of systems [1].

Tips

• Convergence of most iterative methods depends on the condition number of the coefficient matrix, cond(A). When A is square, you can use equilibrate to improve its condition number, 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.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• After execution, code generation for pcg does not display a diagnostic message.

• You can specify only input A as sparse. Other inputs (vector/matrix) must be nonsparse.

• You must enable support for variable-size arrays.

• The number of iterations iter required for the algorithm to converge in the generated code might be different from the iteration values in MATLAB^®.

Thread-Based Environment

Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

This function fully supports thread-based environments. For more information, see Run MATLAB Functions in Thread-Based Environment.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The pcg function supports GPU array input with these usage notes and limitations:

• When input A is a sparse matrix:

  • If you use two preconditioners, M1 and M2, then they must be lower triangular and upper triangular matrices, or both of them must be function handles. Using lower triangular and upper triangular preconditioner matrices instead of function handles can significantly improve computation speed.

  • For GPU arrays, pcg does not detect stagnation (Flag 3). Instead, it reports failure to converge (Flag 1).

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Distributed Arrays

Partition large arrays across the combined memory of your cluster using Parallel Computing Toolbox™.

Usage notes and limitations:

• If M1 is a function, then it is applied independently to each row.

For more information, see Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox).

Version History

Introduced before R2006a