Linear Systems Iterative Solutions

1. Linear SystemsIterative 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 PDE’s • 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 PDE’s • No changes over time (steady state): • Laplace’s 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 PDE’s 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 Laplace’s Equation:

8. Linear System n2 n

9. 3D Simulations n3 n2 n

10. Gaussian Elimination • What happens to these banded matrices when Gaussian Elimination is applied? • Matrix only has about 7n3non-zero elements. • Matrix size = N2, where N=n3 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 n=300 • Memory cost: • 189,000,000 = 189*106 elements • Floats => 756MB • Doubles => 1.4GB • Full matrix would be: 7.29*1014! • Gaussian Elimination • Floats => 1.9*1013MB • With n=300, simulating weather for the state of Ohio would have samples > 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 A=D-M, where D is the diagonal elements of A. • Ax=b • Dx-Mx=b • Dx = Mx+b 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. • Let’s 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 i=0; i<N; i++) { // For each equation double sum = b[i]; for (int j=0; j<N; j++) {// Compute new xi if( i <> j ) sum += -A(i,j)x(j); } temp[i] = sum / A[i,i]; } // Test for convergence … x = temp; } Complexity: Each Iteration: O(N2) Total: O(MN2)

22. Jacobi Method - Complexity while( !converged ) { for (int i=0; i<N; i++) { // For each equation double sum = b[i]; foreach (double element in nonZeroElements[i]) { // Compute new xi if( i <> j ) sum += -A(i,j)*x(j); } temp[i] = sum / A[i,i]; } // Test for convergence … x = temp; } Complexity: Each Iteration: O(pN) Total: O(MpN) p= # non-zero elements For our 2D Laplacian Equation, p=4 N=n2 with n=300 => N=90,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 A=D+L+U, where D is the diagonal elements of A, L is the lower triangular part of A and U is the upper part. • Ax=b • Dx+Lx+Ux=b • (D+L)x = b-Ux (D+L)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 i=0; i<N; i++) { // For each equation double sum = b[i]; foreach (double element in nonZeroElements[i]) { if( i <> j ) sum += -A(i,j)x(j); } x[i] = sum / A[i,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? • 10x+y=12 => high-slope (more vertical) • x+10y=21 => 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