Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sep 2016
Lecture topic: Iterative Methods for Linear Systems

Gaussian Elimination Modules & a Vector Norm

Gaussian Elimination with No Output Module

Here is a routine with no intermediate output:

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Return[A];  ] ;

Matrix Inversion Module

MatInv[A0_]:=Module[{A1,A2,A3,n=Length[A0]},A1=Join[A0,IdentityMatrix[n],2];
A2=GaussElim3[A1];
A3=Take[A2,{1,n},{n+1,2n}];
Return[A3]];

Vector Norm

norm[v_,p_]:=((((∑_(i=1)^(Length[v]))((Abs[v[[i]]])^p)))^(1/p))

Jacobi Iterative Method for Linear Systems

A Small System

Here is a simple linear system we want to solve:

n=4
A={{10,-1,2,0},{-1,11,-1,3},{2,-1,10,-1},{0,3,-1,8}}
b={6,25,-11,15}

Jacobi Iteration

Here is a module that solves the system using Jacobi Iteration:

jacobiiterative[A_,b_,xi_]:=Module[{n=Length[b],x,x0,i,j,k,maxit=1000,tol=(10^(-8))},
x0=N[xi];
x=x0;
For[ k=1, k≤maxit, k=k+1,
 For[i=1,i≤ n,i=i+1,
x[[i]]=(1/(A[[i,i]]))*(b[[i]]-(∑_(j=1)^(i-1))A[[i,j]]x0[[j]]-(∑_(j=i+1)^n)A[[i,j]]x0[[j]]);
];
If[norm[x-x0,2]<tol,Goto[end]];
x0=x;
];
Print["Algorithm did not converge to required tolerence"];
Label[end];Print[" npts = ",k];Return[x];
]

jacobiiterative[A,b,{0,0,0,0}]

Mathematica

Here is Mathematica's solution, using the inverse matrix, and another built-in method using LinearSolve.

mathsol=Inverse[A].b
LinearSolve[A,b]

Gaussian Elimination and Matrix Inversion

Here are a couple of other solutions that don't rely on Mathematica's built-in functions. The first uses our matrix inversion module:

MatInv[A].b

And the following uses Guassian Elimination directly:

matrix1=ConstantArray[0,{n,n}];
Table[matrix1[[i]]=Append[A[[i]],b[[i]]],{i,1,n}];
MatrixForm[matrix1]
sol=Array[f,n];
Do[sol[[i]]=GaussElim3[matrix1][[i,n+1]],{i,1,n}];
sol

For small systems, these both work very well.

Gauss-Seidel Iterative Method for Linear Systems

Small System

Here is a simple system we want to solve:

n=4
A={{10,-1,2,0},{-1,11,-1,3},{2,-1,10,-1},{0,3,-1,8}}
b={6,25,-11,15}

Here is a module that solves the system using Guass-Seidel Iteration. The only change from Jacobi Iteration is in red (previously computed (Subscript[x, i]) values are used immediately).

gaussseideliterative[A_,b_,xi_]:=Module[{n=Length[b],x,x0,i,j,k,maxit=1000,tol=(10^(-8))},
x0=N[xi];
x=x0;
For[ k=1, k≤maxit, k=k+1,
 For[i=1,i≤ n,i=i+1,
x[[i]]=(1/(A[[i,i]]))*(b[[i]]-(∑_(j=1)^(i-1))A[[i,j]]x[[j]]-(∑_(j=i+1)^n)A[[i,j]]x0[[j]]);
];
If[norm[x-x0,2]<tol,Goto[end]];
x0=x;
];
Print["Algorithm did not converge to required tolerence"];
Label[end];Print[" npts = ",k];Return[x];
]

gaussseideliterative[A,b,{0.,0.,0.,0.}]

We see that Gauss-Seidel speeds convergence in this case, since it only took 10 iterations (Jacobi took 24).

Larger System

Here is a larger system we want to solve:

n=40;
SeedRandom[2]
A=Table[RandomReal[{-1,1}],{i,1,n},{j,1,n}];
b=Table[RandomInteger[{-10,10}],{i,1,n}]
x0=Table[RandomReal[{0,1}],{i,1,n}]

Here is Mathematica's solution:

LinearSolve[A,b]

And ours:

gaussseideliterative[A,b,x0]

What went wrong? The Jacobi and Gauss-Seidel techniques are only guaranteed to find the solution if the matrix A is strictly diagonally dominant (for each row in the matrix, the main diagonal entry has a magnitude greater than the sum of the magnitudes of all the other entries in that row).

|(Subscript[A, i,i])|>(∑_(j=1,j≠i)^n)|(Subscript[A, i,j])|

Let’s check this for our matrix:

Table[Abs[A[[i,i]]]>(∑_(j=1)^(i-1))Abs[A[[i,j]]]+(∑_(j=i+1)^(Length[A]))Abs[A[[i,j]]],{i,1,Length[A]}]

Gee, no wonder our method failed!

Essentially, what is happening if the matrix is not diaginally dominant is that the sequence of matrices we construct does not converge.

Larger System that is diagonally dominant

Here is another larger system we want to solve:

n=40;
SeedRandom[2]
A=Table[RandomReal[{-1,1}],{i,1,n},{j,1,n}];
Do[A[[i,i]]=A[[i,i]]+27 i,{i,1,n}];
b=Table[RandomInteger[{-10,10}],{i,1,n}]
x0=Table[RandomReal[{0,1}],{i,1,n}]

Let's check for diagonal dominance:

Table[Abs[A[[i,i]]]>(∑_(j=1)^(i-1))Abs[A[[i,j]]]+(∑_(j=i+1)^(Length[A]))Abs[A[[i,j]]],{i,1,Length[A]}]

So our iterative methods should work

Here is Mathematica's solution:

LinearSolve[A,b]

And ours:

gaussseideliterative[A,b,x0]

These iterative methods work well if

the matrix is large,

has a diagonal without zero elements (or can be manipulated into such form), and

sparse (many entries are zero). This last condition helps with the diagonal dominance aspect.

There are many application for which the off diagonal matrix elements are small compared to the diagonal.

Expressing a matrix as A=D+L+U: 
Diagonal Dominance of A related to the Norm of (Subscript[T, j])=-(D^(-1))(L+U) and (Subscript[T, g])=-(((D+L))^(-1))U

What we are interested in when we do this is the matrix (Subscript[T, j])=-(D^(-1))(L+U) and (Subscript[T, g])=-(((D+L))^(-1))U. Here is an example on a small matrix

Clear[a]
n=4;
SeedRandom[2]
A=Array[a,{n,n}];

L=U=DD=ConstantArray[0,{n,n}];
Do[
{If[i>j,L[[i,j]]=A[[i,j]]];
If[i<j,U[[i,j]]=A[[i,j]]];
If[i==j,DD[[i,i]]=A[[i,i]]];}
,{i,1,n},{j,1,n}]
(Subscript[T, j])=-MatInv[DD].(L+U);
(Subscript[T, g])=-MatInv[DD+L].U;

A//MatrixForm

L//MatrixForm
U//MatrixForm
DD//MatrixForm

Simplify[(Subscript[T, j])]//MatrixForm
Simplify[(Subscript[T, g])]//MatrixForm

What is the matrix norm of T on the larger matrix we looked at earlier that our method failed on?

n=40;
SeedRandom[2]
A=Table[RandomReal[{-1,1}],{i,1,n},{j,1,n}];
L=U=DD=ConstantArray[0,{n,n}];
Do[
{If[i>j,L[[i,j]]=A[[i,j]]];
If[i<j,U[[i,j]]=A[[i,j]]];
If[i==j,DD[[i,i]]=A[[i,i]]];}
,{i,1,n},{j,1,n}]
(Subscript[T, j])=-MatInv[DD].(L+U);
(Subscript[T, g])=-MatInv[DD+L].U;
Print["(Subscript[T, j])"]
Norm[(Subscript[T, j]),1]
Norm[(Subscript[T, j]),2]
Norm[(Subscript[T, j]),∞]
Print["(Subscript[T, g])"]
Norm[(Subscript[T, g]),1]
Norm[(Subscript[T, g]),2]
Norm[(Subscript[T, g]),∞]

This norm is larger than one (any norm we choose to work with), and so it is not surprising that our method did not converge. Here is the result for the larger system which did work:

n=40;
SeedRandom[2]
A=Table[RandomReal[{-1,1}],{i,1,n},{j,1,n}];
Do[A[[i,i]]=A[[i,i]]+27 i,{i,1,n}];
L=U=DD=ConstantArray[0,{n,n}];
Do[
{If[i>j,L[[i,j]]=A[[i,j]]];
If[i<j,U[[i,j]]=A[[i,j]]];
If[i==j,DD[[i,i]]=A[[i,i]]];}
,{i,1,n},{j,1,n}]
(Subscript[T, j])=-MatInv[DD].(L+U);
(Subscript[T, g])=-MatInv[DD+L].U;
Print["(Subscript[T, j])"]
Norm[(Subscript[T, j]),1]
Norm[(Subscript[T, j]),2]
Norm[(Subscript[T, j]),∞]
Print["(Subscript[T, g])"]
Norm[(Subscript[T, g]),1]
Norm[(Subscript[T, g]),2]
Norm[(Subscript[T, g]),∞]

This matrix has norms less than one, which is why the iterative techniques converged for that matrix.

The key result is that a diagonally dominant matrix A=D+L+U will have Jacobi iteration matrix  (Subscript[T, j])=-(D^(-1))(L+U) which has a norm less than 1 (any norm, as we can see!). If this is true, then the Jacobi method will converge when applied to the matrix A. Also, if (Subscript[T, g])=(((D+L))^(-1))U has a norm less than one, then the Gauss-Seidel method will converge when applied to A.

For the example that converges, we saw that the Gauss-Seidel method converged more quickly than the Jacobi method.

Here is the definition of the spectral radius, and we can calculate the spectral radius for (Subscript[T, j]) and (Subscript[T, g]):

ρ[A0_]:=Module[{a1,λ},λ=Eigenvalues[A0];a1=Max[Abs[λ]] ;Return[a1]]

ρ[(Subscript[T, j])]
ρ[(Subscript[T, g])]

We see that the spectral radius is smaller for (Subscript[T, g]), so we would expect Gauss-Seidel to converge more quickly than the Jacobi method. Remember, the spectral radius is a lower bound for any of the natural norms.

For My Reference
Old Notation: Expressing a matrix as A=D-L-U: 
Diagonal Dominance of A related to the Norm of (Subscript[T, j])=(D^(-1))(L+U) and (Subscript[T, g])=(((D-L))^(-1))U

What we are interested in when we do this is the matrix T=(D^(-1))(L+U). Here is an example on a small matrix

Clear[a]
n=4;
SeedRandom[2]
A=Array[a,{n,n}];

L=U=DD=ConstantArray[0,{n,n}];
Do[
{If[i>j,L[[i,j]]=-A[[i,j]]];
If[i<j,U[[i,j]]=-A[[i,j]]];
If[i==j,DD[[i,i]]=A[[i,i]]];}
,{i,1,n},{j,1,n}]
(Subscript[T, j])=MatInv[DD].(L+U);
(Subscript[T, g])=MatInv[DD-L].U;

A//MatrixForm

L//MatrixForm
U//MatrixForm
DD//MatrixForm

(Subscript[T, j])//MatrixForm
(Subscript[T, g])//MatrixForm

What is the matrix norm of T on the larger matrix we looked at earlier that our method failed on?

n=40;
SeedRandom[2]
A=Table[RandomReal[{-1,1}],{i,1,n},{j,1,n}];
L=U=DD=ConstantArray[0,{n,n}];
Do[{If[i>j,L[[i,j]]=-A[[i,j]]];If[i<j,U[[i,j]]=-A[[i,j]]];If[i==j,DD[[i,i]]=A[[i,i]]];},{i,1,n},{j,1,n}]
(Subscript[T, j])=MatInv[DD].(L+U);
(Subscript[T, g])=MatInv[DD-L].U;
Print["(Subscript[T, j])"]
Norm[(Subscript[T, j]),1]
Norm[(Subscript[T, j]),2]
Norm[(Subscript[T, j]),∞]
Print["(Subscript[T, g])"]
Norm[(Subscript[T, g]),1]
Norm[(Subscript[T, g]),2]
Norm[(Subscript[T, g]),∞]

This norm is larger than one (any norm we choose to work with), and so it is not surprising that our method did not converge. Here is the result for the larger system which did work:

n=40;
SeedRandom[2]
A=Table[RandomReal[{-1,1}],{i,1,n},{j,1,n}];
Do[A[[i,i]]=A[[i,i]]+27 i,{i,1,n}];
L=U=DD=ConstantArray[0,{n,n}];
Do[{If[i>j,L[[i,j]]=-A[[i,j]]];If[i<j,U[[i,j]]=-A[[i,j]]];If[i==j,DD[[i,i]]=A[[i,i]]];},{i,1,n},{j,1,n}]
(Subscript[T, j])=MatInv[DD].(L+U);
(Subscript[T, g])=MatInv[DD-L].U;
Print["(Subscript[T, j])"]
Norm[(Subscript[T, j]),1]
Norm[(Subscript[T, j]),2]
Norm[(Subscript[T, j]),∞]
Print["(Subscript[T, g])"]
Norm[(Subscript[T, g]),1]
Norm[(Subscript[T, g]),2]
Norm[(Subscript[T, g]),∞]

The key result is that a diagonally dominant matrix A=D-L-U will have Jacobi iteration matrix  (Subscript[T, j])=(D^(-1))(L+U) which has a norm less than 1 (any norm, as we can see!). If this is true, then the Jacobi method will converge when applied to the matrix A. Also, if (Subscript[T, g])=(((D-L))^(-1))U has a norm less than one, then the Gauss-Seidel method will converge when applied to A.

For the example that converges, we saw that the Gauss-Seidel method converged more quickly than the Jacobi method.

Here is the definition of the spectral radius, and we can calculate the spectral radius for (Subscript[T, j]) and (Subscript[T, g]):

ρ[A0_]:=Module[{a1,λ},λ=Eigenvalues[A0];a1=Max[Abs[λ]] ;Return[a1]]

ρ[(Subscript[T, j])]
ρ[(Subscript[T, g])]

We see that the spectral radius is smaller for (Subscript[T, g]), so we would expect Gauss-Seidel to converge more quickly than the Jacobi method. Remember, the spectral radius is a lower bound for any of the natural norms (Theorem 7.15). See page 444 for a slightly more expanded discussion on choosing an iterative method for a given matrix A.