c6.2. ITERATIVE METHODS �2006 Gilbert Strang 6.2 Iterative Methods New solution methods are needed when a problem Ax = b is too large and expensive for ordinary elimination. We are thinking of sparse matrices A, so that multiplications Ax are relatively cheap. If A has at most p nonzeros in every row, then Ax needs at most pn multiplications. Typical applications are to large finite difference or finite element equations, where we often write A = K. We are turning from elimination to look at iterative methods. There are really two big decisions, the preconditioner P and the choice of the method itself: 1. A good preconditioner P is close to A but much simpler to work with. 2. Options include pure iterations (6.2), multigrid (6.3), and Krylov methods (6.4), including the conjugate gradient method. Pure iterations compute each new xk+1 from xk − P −1(Axk − b). This is called “stationary ” because every step is the same. Convergence to x� = A−1b, studied below, will be fast when all eigenvalues of M = I − P −1A are small. It is easy to suggest a preconditioner, but not so easy to suggest an excellent P (Incomplete LU is a success). The older iterations of Jacobi and Gauss-Seidel are less favored (but they are still important, you will see good points and bad points). Multigrid begins with Jacobi or Gauss-Seidel iterations, for the one job that they do well. They remove high frequency components (rapidly oscillating parts) to leave a smooth error. Then the central idea is to move to a coarser grid —where the rest of the error can be destroyed. Multigrid is often dramatically successful. Krylov spaces contain all combinations of b, Ab, A2b, . . . and Krylov methods look for the best combination. Combined with preconditioning, the result is terrific. When the growing subspaces reach the whole space Rn , those methods give the exact solution A−1b. But in reality we stop much earlier, long before n steps are complete. The conjugate gradient method (for positive definite A, and with a good preconditioner ) has become truly important. The goal of numerical linear algebra is clear: Find a fast stable algorithm that uses the special properties of the matrix. We meet matrices that are symmetric or triangular or orthogonal or tridiagonal or Hessenberg or Givens or Householder. Those are at the core of matrix computations. The algorithm doesn’t need details of the entries (which come from the specific application). By concentrating on the matrix structure, numerical linear algebra offers major help. Overall, elimination with good numbering is the first choice ! But storage and CPU time can become excessive, especially in three dimensions. At that point we turn from elimination to iterative methods, which require more expertise than K\F . The next pages aim to help the reader at this frontier of scientific computing.c�2006 Gilbert Strang Stationary Iterations We begin with old-style pure stationary iteration. The letter K will be reserved for “Krylov” so we leave behind the notation KU = F . The linear system becomes Ax = b. The large sparse matrix A is not necessarily symmetric or positive definite: Linear system Ax = b Residual rk = b − Axk Preconditioner P � A The preconditioner P attempts to be close to A and still allow fast iterations. The Jacobi choice P = diagonal of A is one extreme (fast but not very close). The other extreme is P = A (too close). Splitting the matrix A gives a new form of Ax = b: Splitting P x = (P − A)x + b . (1) This form suggests an iteration, in which every vector xk leads to the next xk+1: Iteration P xk+1 = (P − A)xk + b . (2) Starting from any x0, the first step finds x1 from P x1 = (P − A)x0 + b. The iteration continues to x2 with the same matrix P , so it often helps to know its triangular factors in P = LU . Sometimes P itself is triangular, or L and U are approximations to the triangular factors of A. Two conditions on P make the iteration successful: 1. The new xk+1 must be quickly computable. Equation (2) must be fast to solve. 2. The errors ek = x − xk should approach zero as rapidly as possible. Subtract equation (2) from (1) to find the error equation. It connects ek to ek+1: Error P ek+1 = (P − A)ek which means ek+1 = (I − P−1A)ek = M ek . (3) The right side b disappears in this error equation. Each step multiplies the error vector ek by M . The speed of convergence of xk to x (and of ek to zero) depends entirely on M . The test for convergence is given by the eigenvalues of M : M = I − P −1A |�(M )| < 1.Convergence test Every eigenvalue of must have The largest eigenvalue (in absolute value) is the spectral radius λ(M) = max |�(M )|. Convergence requires λ(M ) < 1. The convergence rate is set by the largest eigen-value. For a large problem, we are happy with λ(M ) = .9 and even λ(M ) = .99. eWhen the initial error e0 happens to be an eigenvector of M , the next error is 1 = M e0 = �e0. At every step the error is multiplied by �. So we must have |�| < 1. Normally e0 is a combination of all the eigenvectors. When the iteration multiplies by M , each eigenvector is multiplied by its own eigenvalue. After k steps those multipliers are �k , and the largest is (λ(M ))k . If we don’t use a preconditioner then M = I − A. All the eigenvalues of A must be inside a unit circle centered at 1, for convergence. Our second difference matrices A = K would fail this test (I − K is too large). The first job of a preconditioner is to get the matrix decently scaled. Jacobi will now give λ(I − 1 K) < 1, and a really 2 good P will do more.6.2. ITERATIVE METHODS c�2006 Gilbert Strang Jacobi Iterations For preconditioner we first propose a simple choice: P = D of AJacobi iteration diagonal part Typical examples have spectral radius λ(M) = 1 − cN−2 , where N counts meshpoints in the longest direction. This comes closer and closer to 1 (too close) as the mesh is refined and N increases. But Jacobi is important, it does part of the job. For our tridiagonal matrices K, Jacobi’s preconditioner is just P = 2I (the diago-1nal of K). The Jacobi iteration matrix becomes M = I − D−1A = I − 2 K: � ⎡ 0 1 ⎢Iteration matrix 1 1 � 1 0 1 M = I − K = � ⎢ . (4)for a Jacobi