Test of iterative methods - Trefethen’s example

Taken from Trefethen, Bau: Numerical Linear Algebra, SIAM 1997, p. 315

Preliminaries

Clear all symbols from previous evaluations to avoid conflicts

Clear["Global`*"]

Problem

Solve Ax = b for a very sparse array

n=1000;
A=SparseArray[{{i_,i_}->0.5+Sqrt[i],{i_,j_}/;Abs[i-j]==1->1.0,{i_,j_}/;Abs[i-j]==100->1.0},{n,n}]
b=ConstantArray[1,n];
xDirect=LinearSolve[A,b];

Norm[A.xDirect-b]

The matrix is symmetric positive definite, all eigenvalues are positive. Thus all methods should converge.

eigenA=Sort[Eigenvalues[A]];
Print["κ(A) = ",Max[eigenA]/Min[eigenA]]
ListPlot[{eigenA},PlotRange->All,AxesLabel->{"n","eigenvalue"}]

Basic iterative methods

In general, we solve a system of linear equations

A x=b

and use a standard decomposition

A=D-(OverTilde[L])-(OverTilde[U])=D(I-L-U)

where D  is the diagonal of A,-(OverTilde[L]) and -(OverTilde[U]) are the strictly lower and upper triangular parts respectively and  L=(D^(-1))(OverTilde[L]), U=(D^(-1))(OverTilde[U]),
More generally will write A=M-K.

One step of the general iterative method can be written as

(Subscript[x, m+1])=R (Subscript[x, m])+c  where  R=(M^(-1))K, c=(M^(-1))b

and the method is convergent if and only if the spectral radius of the matrix R

ρ(R) =max|(Subscript[λ, i])|   where  (Subscript[λ, i]) are eigenvalues of R

satisfies

ρ(R) <1.

Diag=DiagonalMatrix[Diagonal[A]];
D1=Inverse[Diag];
L=-D1.LowerTriangularize[A,-1];
U=-D1.UpperTriangularize[A,1];
Id=IdentityMatrix[n];
(* MatrixForm[U] *)

An initial guess and the maximum number of iterations will be the same for all methods.

x0=ConstantArray[0.0,n];
niter=200;

And here is a common part of all basic iterative algorithms. It is supposed that (Subscript[n, iter]), (Subscript[x, 0]) and (Subscript[x, direct]) are already assigned.

basicIterativeMethod[matrixR_,vectorC_]:=
Module[{x,e},
x=x0;
e=ConstantArray[1.0 (10^(-20)),niter+1];
e[[1]]=Max[Abs[x-xDirect]];
Do[
x=matrixR.x+vectorC;
e[[i+1]]=Max[Abs[x-xDirect]],
{i,1,niter}
];
{x,e}
]

Jacobi method

One step of the Jacobi method is

(Subscript[x, m+1,j])=(1/(Subscript[a, jj]))((Subscript[b, j])-(∑_(k≠j))(Subscript[a, jk])(Subscript[x, m,k])) ,j=1,...,n

or in a matrix form

(Subscript[x, m+1])=(Subscript[R, Jac])(Subscript[x, m])+(Subscript[c, Jac])   where (Subscript[R, Jac])=(D^(-1))((OverTilde[L])+(OverTilde[U]))=L+U  and  (Subscript[c, Jac])=(D^(-1))b

R=L+U;
c=D1.b;
(* MatrixForm[R]*)
eigen=Eigenvalues[R];
ρJac=Max[Abs[eigen]];
Print["ρ((Subscript[R, Jac])) = ",ρJac]
ListPlot[Sort[eigen],PlotRange->All,PlotLegends->{"Eigenvalues of (Subscript[R, Jac])"},AxesLabel->{"n","eigenvalue"}]
{xIter,errorJac}=basicIterativeMethod[R,c];
ListLogPlot[{errorJac},PlotRange->{10^1,10^(-18)},PlotStyle->{Blue},PlotLegends->{"Jacobi"},AxesLabel->{"niter","error"}]

Gauss-Seidel method

One step of the Gauss-Seidel method is

(Subscript[x, m+1,j])=(1/(Subscript[a, jj]))((Subscript[b, j])-(∑_(k=1)^(j-1))(Subscript[a, jk])(Subscript[x, m+1,k])-(∑_(k=j+1)^n)(Subscript[a, jk])(Subscript[x, m,k])) ,j=1,...,n

or in a matrix form

(Subscript[x, m+1])=(Subscript[R, GS])(Subscript[x, m])+(Subscript[c, GS])   where (Subscript[R, GS])=(((D-(OverTilde[L])))^(-1))(OverTilde[U])=(((I-L))^(-1))U  and  (Subscript[c, GS])=(((D-(OverTilde[L])))^(-1))b=(((I-L))^(-1))(D^(-1))b

ImL1=Inverse[Id-L];
R=ImL1.U;
c=ImL1.D1.b;
(*MatrixForm[R]*)
eigen=Eigenvalues[R];
ρGS=Max[Abs[eigen]];
Print["ρ((Subscript[R, GS])) = ",ρGS]
ListPlot[Sort[Re[eigen]],PlotRange->All,PlotLegends->{"Real part of eigenvalues of (Subscript[R, GS])"},AxesLabel->{"n","eigenvalue"}]
{xIter,errorGS}=basicIterativeMethod[R,c];
ListLogPlot[{errorJac,errorGS},PlotRange->{10^1,10^(-18)},PlotStyle->{Blue,Red},PlotLegends->{"Jacobi","Gauss-Seidel"},AxesLabel->{"niter","error"}]

Successive overrelaxation method with optimal ω

One step of the SOR(ω) method is

(Subscript[x, m+1,j])=(1-ω)(Subscript[x, m,j])+(ω/(Subscript[a, jj]))((Subscript[b, j])-(∑_(k=1)^(j-1))(Subscript[a, jk])(Subscript[x, m+1,k])-(∑_(k=j+1)^n)(Subscript[a, jk])(Subscript[x, m,k])) ,j=1,...,n

or in a matrix form

(Subscript[x, m+1])=(Subscript[R, SOR])(Subscript[x, m])+(Subscript[c, SOR])   where (Subscript[R, SOR])=(((D-ω(OverTilde[L])))^(-1))[(1-ω)D+ω(OverTilde[U])]=(((I-ω L))^(-1))[(1-ω)I+ω U] 
and  (Subscript[c, GS])=((ω(D-ω(OverTilde[L])))^(-1))b=((ω(I-ω L))^(-1))(D^(-1))b

ω=2/(1+Sqrt[1-((ρJac)^2)]);
Print["ω = ",ω]
ImL1=Inverse[Id-ω L];
R=ImL1.((1-ω)Id+ω U);
c=ω ImL1.D1.b;
(*MatrixForm[R]*)
eigen=Eigenvalues[R];
ρSOR=Max[Abs[eigen]];
Print["ρ((Subscript[R, SOR])) = ",ρSOR]
ListPlot[Sort[Re[eigen]],PlotLegends->{"Real part of eigenvalues of (Subscript[R, SOR])"},AxesLabel->{"n","eigenvalue"}]
{xIter,errorSORopt}=basicIterativeMethod[R,c];
ListLogPlot[{errorJac,errorGS,errorSORopt},PlotRange->All,PlotStyle->{Blue,Red,Green},PlotLegends->{"Jacobi","Gauss-Seidel","optimal SOR"},AxesLabel->{"niter","error"}]

Successive overrelaxation method with estimated ω

Print["ρ((Subscript[R, Jac])) = ",ρJac]
(* Determination of ρ((Subscript[R, Jac])) from the error *)
nit1=5;nit2=10;
ρJacEst=Exp[(Log[errorJac[[nit2]]]-Log[errorJac[[nit1]]])/(nit2-nit1)];
Print["estimated ρ((Subscript[R, Jac])) = ",ρJacEst]
(* Use this value in SOR algorithm *)
ω=2/(1+Sqrt[1-((ρJacEst)^2)]);
Print["ω = ",ω]
ImL1=Inverse[Id-ω L];
R=ImL1.((1-ω)Id+ω U);
c=ω ImL1.D1.b;
(*MatrixForm[R]*)
eigen=Eigenvalues[R];
ρSORest=Max[Abs[eigen]];
Print["estimated ρ((Subscript[R, SOR])) = ",ρSORest, " (Optimal ρ((Subscript[R, SOR])) = ",ρSOR,")"]
ListPlot[Sort[Re[eigen]],PlotLegends->{"Real part of eigenvalues of (Subscript[R, SOR])"},AxesLabel->{"n","eigenvalue"}]
{xIter,errorSOR}=basicIterativeMethod[R,c];
ListLogPlot[{errorJac,errorGS,errorSORopt,errorSOR},PlotRange->All,PlotStyle->{Blue,Red,Green,Orange},PlotLegends->{"Jacobi","Gauss-Seidel","SOR with optimal ω","SOR with estimated ω"},AxesLabel->{"niter","error"}]

Gradient iterative methods

Solving a system of linear equations

A x=b

with a symmetric positive definite matrix A by searching a minimum of the function

ϕ(x)=(1/2)(x^T)A x-(x^T) b

An initial guess and the maximum number of iterations will be the same as in basic iterative method.

Method of the steepest descent

xIter=x0;
r=b-A.x0;
error=ConstantArray[1.0 (10^(-50)),niter+1];
error[[1]]=Max[Abs[xIter-xDirect]];
Do[
w=A.r;
α=r.r/r.w;
xIter=xIter+α r;
r=r-α w;
error[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
errorMSD=error;
ListLogPlot[{errorJac,errorGS,errorSORopt,errorMSD},PlotRange->All,PlotStyle->{Blue,Red,Green,Brown},PlotLegends->{"Jacobi","Gauss-Seidel","SOR with optimal ω","Method of steepest descent"},AxesLabel->{"niter","error"}]

Conjugate Gradient Method

xIter=x0;
r=b-A.x0;
p=r;
γ=r.r;
error=ConstantArray[1.0 (10^(-50)),niter+1];
error[[1]]=Max[Abs[xIter-xDirect]];
Do[
w=A.p;
α=γ/p.w;
xIter=xIter+α p;
r=r-α w;
β=1/γ;
γ=r.r;
β=β γ;
p=r+β p;
error[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
errorCG=error;
ListLogPlot[{errorJac,errorGS,errorSORopt,errorCG},PlotRange->All,PlotStyle->{Blue,Red,Green,Black},PlotLegends->{"Jacobi","Gauss-Seidel","SOR with optimal ω","Conjugated gradients"},AxesLabel->{"niter","error"}]

Conjugate Gradient Method with Preconditioning

As a preconditioner we use a diagonal of A

M=DiagonalMatrix[Diagonal[A]];
xIter=x0;
r=b-A.x0;
error=ConstantArray[1.0 (10^(-50)),niter+1];
error[[1]]=Max[Abs[xIter-xDirect]];
Do[
z=LinearSolve[M,r];
If[i==1,
p=z;γ=r.z,
β=1/γ;γ=r.z;β=β γ;p=z+β p
];
w=A.p;
α=γ/p.w;
xIter=xIter+α p;
r=r-α w;
error[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
(* Save results for later comparison with other methods *)
errorCGP=error;
ListLogPlot[{errorJac,errorGS,errorSORopt,errorCG,errorCGP},PlotRange->All,PlotStyle->{Blue,Red,Green,Brown,Purple},PlotLegends->{"Jacobi","Gauss-Seidel","SOR with optimal ω","Conjugated gradients","Conjugated gradients with preconditioning"},AxesLabel->{"niter","error"}]

The condition number of A was

eigenA=Sort[Eigenvalues[A]];
Print["κ(A) = ",Max[eigenA]/Min[eigenA]]

The condition number of (M^(-1))A is

eigenMA=Sort[Eigenvalues[Inverse[M].A]];
Print["κ((M^(-1))A) = ",Max[eigenMA]/Min[eigenMA]]

Other Krylov space methods

MINRES equivalent to CR algorithm

Special case of GMRES for symmetric, positive definite matrix A called MINRES = minimal residuals.
The Arnoldi algorithm is replaced by the Lanczos algorithm (only three-term recurrence) and the whole algorithm produces the same approximation as Stiefel’s Conjugate Residual (CR) method which can be written in a very similar way as the CG method (see e.g. report https://web.stanford.edu/group/SOL/reports/SOL-2011-2R.pdf):

xIter=x0;
r=b-A.x0;
s=A.r;
γ=r.s;
p=r;
q=s;
error=ConstantArray[1.0 (10^(-50)),niter+1];
error[[1]]=Max[Abs[xIter-xDirect]];
Do[
α=γ/(q.q);
xIter=xIter+α p;
r=r-α q;
s=A.r;
β=1/γ;
γ=r.s;
β=β γ;
p=r+β p;
q=s+β q;
error[[i+1]]=Max[Abs[xIter-xDirect]],
{i,1,niter}
]
(* Save results for later comparison with other methods *)
errorCR=error;
ListLogPlot[{errorGS,errorSORopt,errorCG,errorCGP,errorCR},PlotRange->All,PlotStyle->{Red,Green,Brown,Purple,Blue},PlotLegends->{"Gauss-Seidel","SOR with optimal ω","Conjugated gradients","Conjugated gradients with preconditioning","Conjugate Residual = MINRES"},AxesLabel->{"niter","error"}]

GMRES

This algorithm should produce the same results as CR algorithm if applied to a symmetric, positive definite system.

xIter=x0;
r=b-A.x0;
rr=Norm[r];
error=ConstantArray[1.0 (10^(-50)),niter];
error[[1]]=Max[Abs[xIter-xDirect]];
q=ConstantArray[0.0,{n,niter}];
H=ConstantArray[0.0,{niter+1,niter}];
q[[All,1]]=r/rr;
Do[
q[[All,i]]=A.q[[All,i-1]];
Do[
H[[j,i-1]]=Conjugate[q[[All,j]]].q[[All,i]];
q[[All,i]]=q[[All,i]]-H[[j,i-1]]*q[[All,j]],
{j,1,i-1}];
H[[i,i-1]]=Norm[q[[All,i]]];
q[[All,i]]=q[[All,i]]/H[[i,i-1]];
y=LeastSquares[H[[1;;i,1;;i-1]],rr*UnitVector[i,1]];
xIter=x0+q[[All,1;;i-1]].y;
error[[i]]=Max[Abs[xIter-xDirect]],
{i,2,niter}
]
(* Save results for later comparison with other methods *)
errorGMRES=error;
ListLogPlot[{errorGS,errorSORopt,errorCG,errorCGP,errorCR,errorGMRES},PlotRange->All,PlotStyle->{Blue,Red,Green,Brown,Purple,Black},PlotLegends->{"Gauss-Seidel","SOR with optimal ω","Conjugated gradients","Conjugated gradients with preconditioning","Conjugate Residual = MINRES","GMRES = MINRES"},AxesLabel->{"niter","error"}]

GMRES with preconditioning

Instead of A we use (M^(-1))A where M is the diagonal

M=DiagonalMatrix[Diagonal[A]];
MA=Inverse[M].A;
xIter=x0;
r=Inverse[M].b-MA.x0;
rr=Norm[r];
error=ConstantArray[1.0 (10^(-50)),niter];
error[[1]]=Max[Abs[xIter-xDirect]];
q=ConstantArray[0.0,{n,niter}];
H=ConstantArray[0.0,{niter+1,niter}];
q[[All,1]]=r/rr;
Do[
q[[All,i]]=MA.q[[All,i-1]];
Do[
H[[j,i-1]]=Conjugate[q[[All,j]]].q[[All,i]];
q[[All,i]]=q[[All,i]]-H[[j,i-1]]*q[[All,j]],
{j,1,i-1}];
H[[i,i-1]]=Norm[q[[All,i]]];
q[[All,i]]=q[[All,i]]/H[[i,i-1]];
y=LeastSquares[H[[1;;i,1;;i-1]],rr*UnitVector[i,1]];
xIter=x0+q[[All,1;;i-1]].y;
error[[i]]=Max[Abs[xIter-xDirect]],
{i,2,niter}
]
(* Save results for later comparison with other methods *)
errorGMRES=error;
ListLogPlot[{errorGS,errorSORopt,errorCG,errorCGP,errorCR,errorGMRES},PlotRange->All,PlotStyle->{Blue,Red,Green,Brown,Purple,Black},PlotLegends->{"Gauss-Seidel","SOR with optimal ω","Conjugated gradients","Conjugated gradients with preconditioning","Conjugate Residual = MINRES","GMRES with preconditioning"},AxesLabel->{"niter","error"}]