Solving $Ax=b$: iterative

This section is appropriate if you have a background in mathematics and have encountered matrices before. If not then feel free to attempt this section however the material assumes a relatively high level of knowledge of linear algebra.

MATLAB has many inbuilt methods for solving $Ax=b$. Many of these are iterative and suitable for different kinds of matrices. If you have information about the structure of $A$ and know which iterative solver is the best to use (see Linear Algebra and its Applications by Gilbert Strang for details), you can specify the solver in MATLAB. The help files will specify the types of matrix to which they should be applied. We will now use two of these solvers based on the ‘conjugate gradient’ method.

Walkthrough

Create a Symmetric Positive Definite (SPD) matrix (and check that it is SPD) by using the following commands:

A=randn(50);

% Make sure A is SPD
A=A*A';

% Verify that this gives us a symmetric matrix
spy(A-A')
find((A-A')>1e-10)

% check positive definite: all eigenvalues > 0
e=eig(A);
min(e)

Create a suitable right-hand side $b$:

b=sum(A,2);
Question

With $A$ and $b$ defined as above, what is the exact solution to $Ax=b$?

Expand for solution
Solution

Every element of $x$ should be exactly 1. If this is not clear, make sure you understand what sum(A,2) is doing, and try writing out a 3x3 system by hand.

Calculate $x$ using the preconditioned Conjugate Gradient method (CG) (see Linear Algebra and its Applications by Gilbert Strang for details if you are interested). Note that $A$ must be SPD to use this solver. Try using it on a non-SPD matrix.

x1=pcg(A,b);
x2=pcg(A,b,1e-6,100);

What does the second command do differently? Look at doc pcg. Now look at the solutions:

format long
max(abs(x1-x2))
plot(x1-x2)

This gives the following figure for the error $x1-x2$:

pcg error

which shows that the solutions are not the same.

Question

Why is this?

Expand for solution
Solution

When running the command

x1=pcg(A,b);

you should have received the following output:

>> x1=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 19) has relative residual 0.0019.

Therefore the method did not converge in the given maximum number of steps. It defaults to 20. The command

x2=pcg(A,b,1e-6,100);

sets the maximum number of iterations at 100, which is enough for the method to converge.

Now clear the workspace

clear

The Biconjugate Gradient method can be used in a similar way. Note that $A$ does not have to be SPD in this case.

A=randn(50);
b=sum(A,2);
x=bicg(A,b,1e-6,100);

plot(x)

Now clear the workspace

clear

To learn more about these iterative methods, see:

doc pcg
doc bicg
doc bicgstab
doc gmres