Linear Systems Iterative Solutions

1
Linear Systems Iterative Solutions

• CSE 541
• Roger Crawfis

2
Sparse Linear Systems

• Computational Science deals with the simulation
  of natural phenomenon, such as weather, blood
  flow, impact collision, earthquake response, etc.
• To simulate these issues such as heat transfer,
  electromagnetic radiation, fluid flow and shock
  wave propagation need to be taken into account.
• Combining initial conditions with general laws of
  physics (conservation of energy and mass), a
  model of these may involve a Partial Differential
  Equation (PDE).

3
Example PDEs

• The Wave equation
• 1D ? 2? /? t2 -c2 ? 2? /? x2
• 3D ? 2? /? t2 -c2 ? 2?
• Note ? 2? ? 2? /? x2 ? 2? /? y2 ? 2? /? z2
• ?(x,y,z,t) is some continuous function of space
  and time (e.g., temperature).

4
Example PDEs

• No changes over time (steady state)
• Laplaces Equation
• This can be solved for very simple geometric
  configurations and initial conditions.
• In general, we need to use the computer to solve
  it.

5
Example PDEs
?Up

• Second derivatives

? middle
?Left
? right
?Down
? 2 ? (? Left ? Right ? Up ? Down 4 ?
Middle ) / h2 0
6
Finite Differences

• Fundamentally we are taking derivatives
• Grid spacing or step size of h.
• Finite-Difference method uses a regular grid.

7
Finite Differences

• A very simple problem
• Find the electrostatic potential inside a box
  whose sides are at a given potential
• Set up a n by n grid on which the potential is
  defined and satisfies Laplaces Equation

8
Linear System
n2
n
9
3D Simulations
n3
n
n2
10
Gaussian Elimination

• What happens to these banded matrices when
  Gaussian Elimination is applied?
• Matrix only has about 7n3 non-zero elements.
• Matrix size N2, where Nn3 or n6 elements
• Gaussian Elimination on these suffers from fill.
• The forward elimination phase will produce n2
  non-zero elements per row, or n5 non-zero
  elements.

11
Memory Costs

• Example n300
• Memory cost
• 189,000,000 189106 elements
• Floats gt 756MB
• Doubles gt 1.4GB
• Full matrix would be 7.291014!
• Gaussian Elimination
• Floats gt 1.91013MB
• With n300, simulating weather for the state of
  Ohio would have samples gt 1km apart.
• Remember, this is h in central differences.

12
Solutions for Sparse Matrices

• Need to keep memory (and computation) low.
• These types of problems motivate the Iterative
  Solutions for Linear Systems.
• Iterate until convergence.

13
Jacobi Iteration

• One of the easiest splitting of the matrix A is
  AD-M, where D is the diagonal elements of A.
• Axb
• Dx-Mxb
• Dx Mxb
• x(k)D-1Mx(k-1)D-1b
• Trivial to compute D-1.

14
Jacobi Iteration

• Another way to understand this is to treat each
  equation separately
• Given the ith equation, solve for xi.
• Assume you know the other variables.
• Use the current guess for the other variables.

15
Jacobi iteration
16
Jacobi Iteration

• Cute, but will it work?
• Algorithms, even mathematical ones, need a
  mathematical framework or analysis.
• Lets first look at a simple example.

17
Jacobi Iteration - Example

• Example system
• Initial guess
• Algorithm

18
Jacobi Iteration - Example
1st iteration 2nd iteration
19
Jacobi Iteration - Example

• x(3) 0.427734375
• y(3) 1.177734375
• z(3) 2.876953125
• x(4) -0.351
• y(4) 0.620
• z(4) 1.814

Actual Solution x 0 y 1 z 2
20
Jacobi Iteration

• Questions
• How many iterations do we need?
• What is our stopping criteria?
• Is it faster than applying Gaussian Elimination?
• Are there round-off errors or other precision and
  robustness issues?

21
Jacobi Method - Implementation

• while( !converged )
• for (int i0 iltN i) // For each equation
• double sum bi
• for (int j0 jltN j) // Compute new xi
• if( i ltgt j )
• sum -A(i,j)x(j)
• tempi sum / Ai,i
• // Test for convergence
• x temp

Complexity Each Iteration O(N2) Total O(MN2)
22
Jacobi Method - Complexity

• while( !converged )
• for (int i0 iltN i) // For each equation
• double sum bi
• foreach (double element in nonZeroElementsi)
• // Compute new xi
• if( i ltgt j )
• sum -A(i,j)x(j)
• tempi sum / Ai,i
• // Test for convergence
• x temp

Complexity Each Iteration O(pN) Total
O(MpN) p non-zero elements
For our 2D Laplacian Equation, p4 Nn2 with
n300 gt N90,000
23
Jacobi Iteration

• Cute, but does it work for all matrices?
• Does it work for all initial guesses?
• Algorithms, even mathematical ones, need a
  mathematical framework or analysis.
• We still do not have this.

24
Gauss-Seidel Iteration

• Split the matrix A into three parts ADLU,
  where D is the diagonal elements of A, L is the
  lower triangular part of A and U is the upper
  part.
• Axb
• DxLxUxb
• (DL)x b-Ux

(DL)x(k)b-Ux(k-1)
25
Gauss-Seidel Iteration

• Another way to understand this is to again treat
  each equation separately
• Given the ith equation, solve for xi.
• Assume you know the other variables.
• Use the most current guess for the other
  variables.

26
Gauss-Seidel Iteration

• Looking at it more simply

This iteration
Last iteration
27
Gauss-Seidel Iteration

• Questions
• How many iterations do we need?
• What is our stopping criteria?
• Is it faster than applying Gaussian Elimination?
• Are there round-off errors or other precision and
  robustness issues?

28
Gauss-Seidel - Implementation

• while( !converged )
• for (int i0 iltN i) // For each equation
• double sum bi
• foreach (double element in nonZeroElementsi)
• if( i ltgt j )
• sum -A(i,j)x(j)
• xi sum / Ai,i
• // Test for convergence
• temp x

Complexity Each Iteration O(pN) Total
O(MpN) p non-zero elements
Differences from Jacobi
29
Convergence

• Jacobi Iteration can be shown to converge from
  any initial guess if A is strictly diagonally
  dominant.
• Diagonally dominant
• Strictly Diagonally dominant

30
Convergence

• Gauss-Seidel can be shown to converge is A is
  symmetric positive definite.

31
Convergence - Jacobi

• Consider the convergence graphically for a 2D
  system

Equation 1
Equation 2
y
x
y
x
Initial guess
32
Convergence - Jacobi

• What if we swap the order of the equations?

Equation 2
Equation 1
x

• Not diagonally dominant
• Same set of equations!

Initial guess
33
Diagonally Dominant

• What does diagonally dominant mean for a 2D
  system?
• 10xy12 gt high-slope (more vertical)
• x10y21 gt low-slope (more horizontal)
• Identity matrix (or any diagonal matrix) would
  have the intersection of a vertical and a
  horizontal line. The b vector controls the
  location of the lines.

34
Convergence Gauss-Seidel
Equation 1
y
Equation 2
x
y
x
Initial guess
35
Convergence - SOR

• Successive Over-Relaxation (SOR) just adds an
  extrapolation step.
• w 1.3 impliesgo an extra30

Equation 1
Extrapolation
y
x
Equation 2
Extrapolation
y
x
Initial guess
This would Extrapolation at the very end (mix of
Jacobi and Gauss-Seidel.
36
Convergence - SOR

• SOR with Gauss-Seidel

Equation 1
x
Equation 2
y
x
Initial guess
Extrapolation in bold