Heath - Scientific Computing (523150), страница 93
Текст из файла (страница 93)
Technically, to preserve symmetry, we should apply the conjugate gradientalgorithm to L−1 AL−T instead of M −1 A, where M = LLT . However, the algorithmcan be suitably rearranged so that only M is used and the corresponding matrix L is notrequired explicitly. The resulting preconditioned conjugate gradient algorithm is given here.Starting with an initial guess x0 and taking r0 = b − Ax0 and s0 = M −1 r0 , the following11.5. ITERATIVE METHODS FOR LINEAR SYSTEMS347steps are repeated for k = 0, 1, . . . until convergence:1.2.3.4.5.αk = rkT M −1 rk /sTk Ask .xk+1 = xk + αk sk .rk+1 = rk − αk Ask .T M −1 rT−1 r .βk+1 = rk+1k+1 /rk Mk−1sk+1 = M rk+1 + βk+1 sk .Note that in addition to the one matrix-vector multiplication, Ask , per iteration, we mustalso apply the preconditioner, M −1 rk , once per iteration.The choice of an appropriate preconditioner depends on the usual trade-off between thegain in the convergence rate and the increased cost per iteration that results from applyingthe preconditioner.
Many different choices of preconditioner have been proposed, and thistopic is an active area of research. Some of the types of preconditioning most commonlyused are:• Diagonal (also called Jacobi): M is taken to be a diagonal matrix with diagonal entriesequal to those of A.• Block diagonal (or block Jacobi): If the indices 1, . . .
, n are partitioned into mutuallydisjoint subsets, then mij = aij if i and j are in the same subset, and mij = 0 otherwise.Natural choices include partitioning along lines or planes in two- or three-dimensionalgrids, respectively, or grouping together physical variables that correspond to a commonnode, as in many finite element problems.• SSOR: Using a matrix splitting of the form A = L + D + LT as in Section 11.5.1, we cantake M = (D + L)D −1 (D + L)T , or, introducing the SSOR relaxation parameter ω,M (ω) =1111( D + L)( D)−1 ( D + L)T .2−ω ωωωWith optimal choice of ω, the SSORp preconditioner is capable of reducing the conditionnumber to cond(M −1 A) = O( cond(A) ), but as usual, obtaining knowledge of thisoptimal value may be impractical.• Incomplete factorization: Ideally, one would like to solve the linear system directly usingthe Cholesky factorization A = LLT , but this may incur unacceptable fill (see Section 11.4.1).
One may instead compute an approximate factorization A ≈ L̂L̂T thatallows little or no fill (e.g., restricting the nonzero entries of L̂ to be in the same positionsas those of the lower triangle of A), then use M = L̂L̂T as a preconditioner.• Polynomial : M −1 is taken to be a polynomial in A that approximates A−1 . One wayto obtain a suitable polynomial is to use a fixed number of steps of a stationary iterative method to solve the preconditioning system M zk = rk at each conjugate gradientiteration.• Approximate inverse: M −1 is determined by using an optimization algorithm to minimizethe residualkI − AM −1 kor kI − M −1 Akin some norm, with M −1 restricted to have a prescribed pattern of nonzero entries.348CHAPTER 11. PARTIAL DIFFERENTIAL EQUATIONSNote that some of these preconditioners require a significant amount of work to form theminitially, and this work must also be included in the cost trade-off mentioned earlier.
Theconjugate gradient method is rarely used without some form of preconditioning. Since diagonal preconditioning requires almost no extra work or storage, at least this much preconditioning is always advisable, and more sophisticated preconditioners are often worthwhile.The conjugate gradient method is generally applicable only to symmetric positive definite systems. If the matrix A is indefinite or nonsymmetric, then the algorithm may breakdown both theoretically (e.g., the corresponding optimization problem may not have a minimum) and practically (e.g., the formula for α may fail).
The method can be generalized tosymmetric indefinite systems, as in the SYMMLQ algorithm of Paige and Saunders [198],for example. The conjugate gradient method cannot be generalized to nonsymmetric systems, however, without sacrificing at least one of the two properties—the short recurrenceproperty and the minimum residual property—that largely account for its effectiveness.Nevertheless, in recent years a number of related algorithms have been formulated for solving nonsymmetric linear systems, including GMRES, QMR, CGS, BiCG, Bi-CGSTAB, andothers. These algorithms tend to be significantly less robust or require considerably morestorage or work than the conjugate gradient algorithm that inspired them, but in manycases they are still the most effective methods available for solving very large nonsymmetricsystems.Example 11.7 Iterative Methods for Linear Systems.
We illustrate various iterativemethods by using them to solve the 4 × 4 linear system for the Laplace equation on the unitsquare in Example 11.4. In each case we take x(0) = 0 as starting guess.The Jacobi method gives the following sequence of iterates for this problem:k0123456789x10.0000.0000.0620.0940.1090.1170.1210.1230.1240.125x20.0000.0000.0620.0940.1090.1170.1210.1230.1240.125x30.0000.2500.3120.3440.3590.3670.3710.3730.3740.375x40.0000.2500.3120.3440.3590.3670.3710.3730.3740.375As expected, the Gauss-Seidel method converges somewhat faster, giving the followingsequence of iterates for this problem:k0123456x10.0000.0000.0620.1090.1210.1240.125x20.0000.0000.0940.1170.1230.1250.125x30.0000.2500.3440.3670.3730.3750.375x40.0000.3120.3590.3710.3740.3750.37511.5. ITERATIVE METHODS FOR LINEAR SYSTEMS349The optimal acceleration parameter in the SOR method turns out to be ω = 1.072 for thisproblem, which is so close to 1 that it converges only slightly faster than Gauss-Seidel,giving the following sequence of iterates:k012345x10.0000.0000.0720.1190.1230.125x20.0000.0000.1080.1210.1240.125x30.0000.2680.3560.3710.3740.375x40.0000.3350.3650.3730.3750.375Finally, the conjugate gradient method converges in only two iterations for this problem,giving the following sequence of iterates:k01211.5.6x10.0000.0000.125x20.0000.0000.125x30.0000.3330.375x40.0000.3330.375Rate of ConvergenceExample 11.7 is too small for the results to be representative of the relative performance ofthe methods for problems of practical size.
Recall from the discussion of convergence ratesin Section 5.1.2 that, asymptotically, a linearly convergent sequence with constant C gains− log10 (C) decimal digits per iteration. Thus, the quantity R = − log10 (ρ(G)) serves as auseful quantitative measure of the speed of convergence of a stationary iterative methodwith iteration matrix G.
In this context, R is sometimes called the rate of convergence, butthis term should not be confused the convergence rate r defined in Section 5.1.2.If we use the same five-point finite difference approximation as in the previous examplefor the Laplace equation on the unit square, but with an arbitrary k × k grid of interiormesh points with mesh size h = 1/(k + 1), then the spectral radius ρ(G) and approximaterate of convergence R for the stationary iterative methods are as shown in Table 11.1.
Wesee that the rates of convergence for Jacobi and Gauss-Seidel are proportional to the squareof the mesh size, or equivalently, that the number of iterations per digit of accuracy gainedis proportional to the number of mesh points. The constants of proportionality also tell usthat Gauss-Seidel is asymptotically twice as fast as Jacobi for this model problem. OptimalSOR, on the other hand, is an order of magnitude faster than either of the other methods, asits rate of convergence is proportional to the mesh size, and hence the number of iterationsper digit gained is proportional to the number of mesh points along one side of the grid.Table 11.1: Spectral radius and rate of convergence for k × k grid problemMethodρ(G)R2Jacobicos(πh)(π / log 10)h2 /22Gauss-Seidelcos (πh)(π 2 / log 10)h2Optimal SOR (1 − sin(πh))/(1 + sin(πh))(2π/ log 10)h350CHAPTER 11.
PARTIAL DIFFERENTIAL EQUATIONSTo make these results more concrete, in Table 11.2 are listed the spectral radius andrate of convergence for each method for a range of values of k. We see that the spectralradius is extremely close to 1 for large values of k, and hence all three methods convergevery slowly. From the rate of convergence R, we see that for k = 10 (a linear system oforder 100), Jacobi requires more than 50 iterations to gain a single decimal digit of accuracy,Gauss-Seidel requires more than 25 iterations, and optimal SOR requires about 4 iterations.For k = 100 (a linear system of order 10,000), to gain a single decimal digit of accuracyJacobi requires about 5000 iterations, Gauss-Seidel about 2500, and optimal SOR about37.
Thus, the Jacobi and Gauss-Seidel methods are impractical for a problem of this size,and optimal SOR, though perhaps reasonable for this problem, also becomes prohibitivelyslow for still larger problems. Moreover, the performance of SOR depends on knowledgeof the optimal value for the relaxation parameter ω, which is known analytically for thissimple model problem to be ω = 2/(1 + sin(πh)), but which is much harder to determinein general.Table 11.2: Spectral radius of iteration matrix and rate of convergence for k × k gridJacobiGauss-SeidelOptimal SORk ρ(G)Rρ(G)Rρ(G)R10 0.95950.0180.92060.0360.56040.25250 0.99810.00080.99620.00160.88400.0535100 0.99950.00020.99900.00040.93970.0270500 0.99998 0.0000085 0.99996 0.000017 0.98754 0.005447The convergence behavior of the nonstationary conjugate gradient method is more complicated, but roughly speaking the error is reduced at each iteration by a factor of√κ−1√κ+1on average, whereκ = cond(A) = kAk2 · kA−1 k2 =λmax (A).λmin (A)When the matrix A is well-conditioned (κ ≈ 1), convergence is rapid; but if A is veryill-conditioned (κ 1), then convergence can be arbitrarily slow.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
