Iterative Solution Methods

Iterative Solution Methods



Starts with an initial approximation for the
solution vector (x0)
At each iteration updates the x vector by using
the sytem Ax=b
During the iterations A, matrix is not changed
so sparcity is preserved

Each iteration involves a matrix-vector product

If A is sparse this product is efficiently done
1
Iterative solution procedure




Write the system Ax=b in an equivalent form
x=Ex+f (like x=g(x) for fixed-point iteration)
Starting with x0, generate a sequence of
approximations {xk} iteratively by
xk+1=Exk+f
Representation of E and f depends on the type
of the method used
But for every method E and f are obtained from
A and b, but in a different way
2
Convergence




As k, the sequence {xk} converges to the
solution vector under some conditions on E
matrix
This imposes different conditions on A matrix
for different methods
For the same A matrix, one method may
converge while the other may diverge
Therefore for each method the relation
between A and E should be found to decide on
the convergence
3
Different Iterative methods



Jacobi Iteration
Gauss-Seidel Iteration
Successive Over Relaxation (S.O.R)


SOR is a method used to accelerate the convergence
Gauss-Seidel Iteration is a special case of SOR
method
4
Jacobi iteration
a11 x1  a12 x2    a1n xn  b1
a21 x1  a22 x2    a2 n xn  b2

an1 x1  an 2 x2    ann xn  bn
 x10 
 0
x2 
0

x 

 0
 xn 
1
1 
(b1  a12 x20    a1n xn0 )
k 1
xi  bi 
a11
aii 
1
x12 
(b2  a21 x10  a23 x30    a2 n xn0 )
a22
1
x1n 
(bn  an1 x10  an 2 x20    ann1 xn01 )
ann
x11 

aij x   aij x 

j 1
j i 1

i 1
n
k
j
k
j
5
xk+1=Exk+f iteration for Jacobi method
A can be written as A=L+D+U (not decomposition)
0 0 a11 0
0  0 a12 a13 
 a11 a12 a13   0
a
  a
 0 a
  0 0 a 
a
a
0
0
0
22
23 
 21 22 23   21
 
 
a31 a32 a33  a31 a32 0  0 0 a33  0 0
0 
Ax=b  (L+D+U)x=b
1
k 1
xi 
aii
Dxk+1
i 1
n


k
k
bi   aij x j   aij x j 
j 1
j i 1


 
Lx k
Uxk
Dxk+1 =-(L+U)xk+b
xk+1=-D-1(L+U)xk+D-1b
E=-D-1(L+U)
f=D-1b
6
Gauss-Seidel (GS) iteration
Use the latest
update
a11 x1  a12 x2    a1n xn  b1
a21 x1  a22 x2    a2 n xn  b2

 x10 
 0
x2 
0

x 

 0
 xn 
an1 x1  an 2 x2    ann xn  bn
1
1 
(b1  a12 x20    a1n xn0 )
k 1
xi  bi 
a11
aii 
1
1
1
0
0
x2 
(b2  a21x1  a23 x3    a2 n xn )
a22
1
x1n 
(bn  an1 x11  an 2 x12    ann1 x1n 1 )
ann
x11 
i 1
a x
j 1
ij
k 1
j

  aij x 
j i 1

n
k
j
7
x(k+1)=Ex(k)+f iteration for Gauss-Seidel
Ax=b  (L+D+U)x=b
k 1
i
x
1

aii
Dxk+1
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i 1


 
Lx k 1
(D+L)xk+1 =-Uxk+b
Uxk
xk+1=-(D+L)-1Uxk+(D+L)-1b
E=-(D+L)-1U
f=-(D+L)-1b
8
Comparison



Gauss-Seidel iteration converges more rapidly
than the Jacobi iteration since it uses the latest
updates
But there are some cases that Jacobi iteration
does converge but Gauss-Seidel does not
To accelerate the Gauss-Seidel method even
further, successive over relaxation method can
be used
9
Successive Over Relaxation Method

GS iteration can be also written as follows
k 1
i
x
1
x 
aii
k
i
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i


xik 1  xik   ik
Correction term
 i2
xi3
xi2
xi1
x
0
i

 i1
2
i
 i1
 i0
Multiply with
 1
Faster
convergence
 i0
10
SOR
xik 1  xik   ik
k 1
i
x
1
 x 
aii
k
i
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i


1
k 1
k
xi  (1   ) xi  
aii
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i 1


1<<2 over relaxation (faster convergence)
0<<1 under relaxation (slower convergence)
There is an optimum value for 
Find it by trial and error (usually around 1.6)
11
x(k+1)=Ex(k)+f iteration for SOR
k 1
i
x
1
 (1   ) x  
aii
k
i
i 1
n


k 1
k
bi   aij x j   aij x j 
j 1
j i 1


Dxk+1=(1-)Dxk+b-Lxk+1-Uxk
(D+ L)xk+1=[(1-)D-U]xk+b
E=(D+ L)-1[(1-)D-U]
f= (D+ L)-1b
12
The Conjugate Gradient Method
d 0  r0  b  Ax0
T
ri ri
i  T
d i Adi
xi 1  xi   i Adi
• Converges if A is a
symmetric positive
definite matrix
• Convergence is
faster
T
i 1 i 1
T
i i
r r
i 1 
r r
d i 1  ri 1  i 1d i
13
Convergence of Iterative Methods
Define the solution vector as x̂
k
Define an error vector as e
x  e  xˆ
k
k
Substitute this into x
k 1
 Ex  f
k
e k 1  xˆ  E (e k  xˆ )  f  Exˆ  f  Eek
e
k 1
 Ee  EEe
k
k 1
 EEEe
k 2
E
( k 1) 0
e
14
Convergence of Iterative Methods
iteration
e
k 1
 E
( k 1) 0
e  E
( k 1)
e
0
power
The iterative method will converge for any initial
iteration vector if the following condition is satisfied
Convergence condition
Lim e k 1  0 if
k 
Lim E ( k 1)  0
k 
15
Norm of a vector
A vector norm should satisfy these conditions
x  0 for every nonzero vector x
x  0 iff x is a zero vector
αx   x
for scalar α
x y  x  y
Vector norms can be defined in different forms as
long as the norm definition satisfies these conditions
16
Commonly used vector norms
Sum norm or ℓ1 norm
x 1  x1  x2    xn
Euclidean norm or ℓ2 norm
x 2  x  x   x
2
1
2
2
2
n
Maximum norm or ℓ norm
x   max i xi
17
Norm of a matrix
A matrix norm should satisfy these conditions
A 0
A  0 iff A is a zero matrix
αA   A
for scalar α
A B  A  B
Important identitiy
Ax  A x
x is a vector
18
Commonly used matrix norms
Maximum column-sum norm or ℓ1 norm
m
A 1  max  aij
1 j  n
i 1
Spectral norm or ℓ2 norm
A 2  maximum eigenvalue of A A
T
Maximum row-sum norm or ℓ norm
n
A   max  aij
1i  m
j 1
19
Example

Compute the ℓ1 and ℓ norms of the matrix
3 9 5
7 2 4 


6 8 1
16
19
17  A 
13
15
10
 A1
20
Convergence condition
lim e
k 
k 1
 0 if
lim E
k 
( k 1)
0
Express E in terms of modal matrix P and 
:Diagonal matrix with eigenvalues of E on the diagonal
E  PP
1
E ( k 1)  PP 1 PP 1  PP 1
E ( k 1)  P( k 1) P 1
1k 1



k 1

2

k 1  




k 1 

n 


lim E ( k 1)  0  lim P( k 1) P 1  0  lim ( k 1)  0
k 
k 
k 
 lim ki 1  0   i  1 for i  1,2,...,n
k 
21
Sufficient condition for convergence
If the magnitude of all eigenvalues of iteration matrix
E is less than 1 than the iteration is convergent
It is easier to compute the norm of a matrix than
to compute its eigenvalues
Ex  x
Ex   x 
  x  E x    E
Ex  E x 
E  1 is a sufficient condition for convergence
22
Convergence of Jacobi iteration
E=-D-1(L+U)

 0

 a21
 a22
E 



 an1
 a
 nn
a12

a11
0



a23

a22







ann1

ann
a1n 

a11 

a2 n 

a22 
 
an 1n 


an 1n 1 

0 

23
Convergence of Jacobi iteration
Evaluate the infinity(maximum row sum) norm of E
E

n
aij
j 1
i j
aii
1 
 1 for i  1,2,...,n
n
 aii   aij
j 1
i j
Diagonally dominant matrix
If A is a diagonally dominant matrix, then Jacobi
iteration converges for any initial vector
24
Stopping Criteria

Ax=b

At any iteration k, the residual term is
rk=b-Axk

Check the norm of the residual term
||b-Axk||

If it is less than a threshold value stop
25
Example 1 (Jacobi Iteration)
 4  1 1  x1   7 
 4  8 1  x    21

 2  

 2 1 5  x3   15 
0 
x 0  0
0
b  Ax0
2
 26.7395
Diagonally dominant matrix
7  x20  x30
7
x 
  1.75
4
4
0
0
21

4
x

x
21
1
3
x12 
  2.625
8
8
0
0
15

2
x

x
15
1
2
x31 
  3 .0
5
5
1
1
b  Ax1  10.0452
2
26
Example 1 continued...
7  x12  x31
7  2.625  3
x 

 1.65625
4
4
21  4 x11  x31 21  4 1.75  3
2
x2 

 3.875
8
8
15  2 x11  x12 15  2 1.75  2.625
2
x3 

 4.225
5
5
2
1
7  3.875  4.225
 1.6625
4
21  4 1.65625  4.225
x23 
 3.98125
8
15  2 1.65625  3.875
x33 
 2.8875
5
b  Ax 2
2
 6.7413
x13 
b  Ax2
2
 1.9534
Matrix is diagonally dominant, Jacobi iterations are converging
27
Example 2
 2 1 5  x1   15 
 4  8 1  x    21

 2  

 4  1 1  x3   7 
0 
x 0  0
0
b  Ax0
2
 26.7395
The matrix is not diagonally dominant
 15  x20  5 x30  15
x 

  7 .5
2
2
21  4 x10  x30
21
1
x2 
  2.625
8
8
x31  7  4 x10  x20
 7 .0
1
1
b  Ax1  54.8546
2
28
Example 2 continued...
 15  2.625  5  7
 11.3125
2
21  4  7.5  7
x12 
 0.25
8
x31  7  4  7.5  2.625  39.625
x11 
b  Ax 2
2
 208.3761
The residual term is increasing at each iteration,
so the iterations are diverging.
Note that the matrix is not diagonally dominant
29
Convergence of Gauss-Seidel iteration



GS iteration converges for any initial vector if A
is a diagonally dominant matrix
GS iteration converges for any initial vector if A
is a symmetric and positive definite matrix
Matrix A is positive definite if
xTAx>0 for every nonzero x vector
30
Positive Definite Matrices



A matrix is positive definite if all its eigenvalues
are positive
A symmetric diagonally dominant matrix with
positive diagonal entries is positive definite
If a matrix is positive definite


All the diagonal entries are positive
The largest (in magnitude) element of the whole
matrix must lie on the diagonal
31
Positive Definitiness Check
20 12 25
12 15 2 


25 2 5 
Not positive definite
Largest element is not on the diagonal
5
20 12
12  15 2 


 5
2 25
Not positive definite
All diagonal entries are not positive
20 12 5  Positive definite
12 15 2 

 Symmetric, diagonally dominant, all
 5 2 25 diagonal entries are positive
32
Positive Definitiness Check
20 12 5 
12 15 2 


 8 2 25
A decision can not be made just by investigating
the matrix.
The matrix is diagonally dominant and all diagonal
entries are positive but it is not symmetric.
To decide, check if all the eigenvalues are positive
33
Example (Gauss-Seidel Iteration)
 4  1 1  x1   7 
 4  8 1  x    21

 2  

 2 1 5  x3   15 
0 
x 0  0
0
b  Ax0
2
 26.7395
Diagonally dominant matrix
7
7  x20  x30

 1.75
x 
4
4
21  4 x11  x30 21  4 1.75
1

 3.5
x2 
8
8
15  2 x11  x12 15  2 1.75  3.5
1

 3.0
x3 
5
5
1
1
b  Ax1  3.0414
2
b  Ax1  10.0452
2
Jacobi iteration
34
Example 1 continued...
7  x12  x31
7  3.5  3
x 

 1.875
4
4
21  4 x12  x31 21  4 1.875  3
2
x2 

 3.9375
8
8
15  2 x12  x22 15  2 1.875  3.9375
2
x3 

 2.9625
5
5
2
1
b  Ax 2
b  Ax 2
2
2
 0.4765
 6.7413
Jacobi iteration
When both Jacobi and Gauss-Seidel iterations
converge, Gauss-Seidel converges faster
35
Convergence of SOR method



If 0<<2, SOR method converges for any
initial vector if A matrix is symmetric and
positive definite
If >2, SOR method diverges
If 0<<1, SOR method converges but the
convergence rate is slower (deceleration) than
the Gauss-Seidel method.
36
Operation count




The operation count for Gaussian Elimination or
LU Decomposition was 0 (n3), order of n3.
For iterative methods, the number of scalar
multiplications is 0 (n2) at each iteration.
If the total number of iterations required for
convergence is much less than n, then iterative
methods are more efficient than direct
methods.
Also iterative methods are well suited for
sparse matrices
37