4 Iterative Methods for Solving Linear Systems

1 4 IterativeMethodsforSolvingLinearSystemsI terative Methods formally yield the solution x of a Linear system after aninfinite number of steps. At each step they require the computationof theresidual of the system. In the case of a full matrix, their computationalcost istherefore of the order ofn2operations for each iteration, to be compared withan overall cost of the order of23n3operations needed by direct Methods . Itera-tive Methods can therefore become competitivewith direct Methods providedthe number of iterations that are required toconverge (within a prescribedtolerance) is either independent ofn or scales sublinearly with respect the case of large sparse matrices, as discussed in Section , directmethods may be unconvenient due to the dramatic fill-in, although extremelyefficient direct solverscan be devised on sparse matrices featuringspecialstructures like, for example, those encountered in the approximation ofpartialdifferential equations (see Chapters 12 and 13).

2 Finally, we notice that,when A is ill-conditioned, a combined use of directand Iterative methodsis made possible by preconditioning techniques thatwill be addressed in Section basic idea of Iterative Methods is to construct a sequence of vectors x(k)that enjoy the propertyof convergencex = limk x(k),( )where x is the solution to ( ).In practice, the Iterative process is stopped atthe minimum value ofn such that x(n) x < , where is a fixed toleranceand is any convenient vector norm. However, since the exact solution isobviously not available,it is necessary to introduce suitable stopping criteriato monitor the convergence of the iteration (see Section ).126 4 Iterative Methods for Solving Linear SystemsTo start with, we consider Iterative Methods of the formgiven x(0), x(k+1)= Bx(k)+ f,k 0,( )having denoted by B ann n square matrix called the iteration matrix andby f a vector that is obtained from the right handside An Iterative method of the form ( ) is said to be consistentwith ( ) if f and B are such that x = Bx + f.

3 Equivalently,f = (I B)A 1b. Having denoted bye(k)= x(k) x( )the error at thek-th step of the iteration, the condition for convergence ( )amounts to requiring that limk e(k)= 0 for any choice of theinitial datumx(0)(often called the initial guess).Consistency alone doesnot suffice to ensure theconvergence of the iterativemethod ( ), as shownin the following To solve the Linear system 2Ix =b, consider the iterativemethodx(k+1)= x(k)+b,which is obviously consistent. This scheme isnot convergent for anychoice of theinitial guess. If, for instance,x(0)=0, the method generatesthe sequencex(2k)=0,x(2k+1)=b,k = 0, 1,..On the other hand, ifx(0)=12b the method is convergent. Theorem Let ( ) be a consistent , the sequence of vectors x(k) converges to the solution of ( ) for any choice of x(0)iff (B)< ( ) and the consistency assumption, the recursive relatione(k+1)=Be(k)is obtained.

4 Therefore,e(k)= Bke(0), k = 0, 1,..( )Thus, thanks to Theorem , it follows that limk Bke(0)=0 for anye(0)iff (B)< , suppose that (B)> 1, then there exists at least one eigenvalue (B) with module greater than 1. Lete(0)be an eigenvector associated with ; thenBe(0)= e(0)and, therefore,e(k)= ke(0). As a consequence,e(k)cannot tend to0 ask , since| |> 1. From ( ) and Theorem it follows that asufficient condition forconver-gence to hold is that B < 1, for any consistent matrix norm. It is On the Convergence of Iterative Methods 127to expect that the convergence is faster when (B) is smaller so that an es-timate of (B) might provide a sound indication of the convergence of thealgorithm.

5 Other remarkable quantities in convergence analysis are containedin the following Let B be the iterationmatrix. We call:1. Bm the convergence factor afterm steps of the iteration;2. Bm 1/mthe average convergence factor afterm steps; (B) = 1mlog Bm the average convergence rate afterm steps. These quantities are too expensive to computesince they require evaluatingBm. Therefore, it is usuallypreferred to estimate the asymptotic convergencerate, which is defined asR(B) = limk Rk(B) = log (B),( )where Property has been accounted particular, if B were sym-metric, we would haveRm(B) = 1mlog Bm 2= log (B).In the case of nonsymmetric matrices, (B) sometimes provides an overop-timistic estimate of Bm 1/m(see [Axe94], Section ).

6 Indeed, although (B)< 1, the convergence tozero of the sequence Bm might be non-monotone (see Exercise1). We finally notice that, due to ( ), (B) is theasymptotic convergencefactor. Criteria for estimating the quantities definedso far will be addressedin Section The iterations introduced in ( ) are a special instance ofiterative Methods of the formx(0)= f0(A, b),x(n+1)= fn+1(x(n), x(n 1),.., x(n m), A, b), forn m,where fiand x(m),.., x(1)are given functions andvectors, respectively. Thenumber of steps whichthe current iteration depends on is called the order ofthe method. If the functions fiare independent of the step indexi, the methodis called stationary, otherwise it is nonstationary.

7 Finally, if fidepends linearlyon x(0),.., x(m), the method is called Linear , otherwise it is the light of these definitions, the Methods considered so far are thereforestationary Linear Iterative Methods of first order. In Section , examples ofnonstationary Linear Methods will be provided. 128 4 Iterative Methods for Solving Linear general technique todevise consistent lineariterative Methods is based onan additive splitting of the matrix A of theform A=P N, where P and N aretwo suitable matrices and P is nonsingular. For reasons that will be clear inthe later sections, P iscalled preconditioning matrix or , given x(0), one can compute x(k)fork 1, Solving the systemsPx(k+1)= Nx(k)+ b,k 0.

8 ( )The iteration matrix ofmethod ( ) is B = P 1N, while f = P 1b. Alterna-tively, ( ) can be written in the formx(k+1)= x(k)+ P 1r(k),( )wherer(k)= b Ax(k)( )denotes the residual vector at stepk. Relation ( ) outlines the fact that alinear system, with coefficient matrix P, must besolved to update the solutionat stepk +1. Thus P, besides being nonsingular, oughtto be easily invertible,in order to keep the overall computational cost low. (Notice that, if Pwereequal to A and N=0, method ( ) would converge in one iteration, but atthe same cost of a direct method).Let us mention two results that ensure convergence of the iteration ( ),provided suitable conditions on the splitting ofA are fulfilled (for theirproof,we refer to [Hac94]).

9 Property Let A = P N, with A and P symmetric and the matrix 2P A is positive definite, then the Iterative method defined in( ) is convergent for any choice of the initial datum x(0)and (B) = B A= B P< , the convergence of the iteration ismonotone with respectto thenorms Pand A( , e(k+1) P< e(k) Pand e(k+1) A< e(k) Ak = 0,1,..).Property Let A = P N with A being symmetric and positive the matrix P+PT A is positive definite, then P is invertible, the iterativemethod defined in ( ) is monotonically convergent with respect to norm Aand (B) B A< Jacobi, Gauss-Seidel and Relaxation MethodsIn this section we consider some classical Linear Iterative the diagonal entries of A are nonzero, we cansingle out in each equationthe corresponding unknown, obtaining the equivalent Linear Linear Iterative Methods 129xi=1aii bi n j=1j =iaijxj ,i = 1.

10 ,n.( )In the Jacobi method,once an arbitrarily initial guess x(0)has been chosen,x(k+1)is computed by the formulaex(k+1)i=1aii bi n j=1j =iaijx(k)j ,i = 1,..,n.( )This amounts to performing the following splitting for AP = D, N = D A = E + F,where D is the diagonal matrix of the diagonal entries of A, E is thelowertriangular matrix of entrieseij= aijifi>j,eij= 0 ifi j, and F is theupper triangular matrix of entriesfij= aijifj>i,fij= 0 ifj i. As aconsequence, A = D (E + F).The iteration matrix ofthe Jacobi method is thus given byBJ= D 1(E + F) = I D 1A.( )A generalization of the Jacobi method is the over-relaxation method(or JOR), in which, having introduced a relaxation parameter , ( ) isreplaced byx(k+1)i= aii bi n j=1j =iaijx(k)j + (1 )x(k)i,i = 1.