1 / 36

Linear Systems Iterative Solutions

Linear Systems Iterative Solutions. CSE 541 Roger Crawfis. Sparse Linear Systems. Computational Science deals with the simulation of natural phenomenon, such as weather, blood flow, impact collision, earthquake response, etc.

fraley
Download Presentation

Linear Systems Iterative Solutions

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  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

More Related