Krylov Subspace Solvers and Preconditioners

. In these lecture notes an introduction to Krylov subspace solvers and preconditioners is presented. After a discretization of partial diﬀerential equations large, sparse systems of linear equations have to be solved. Fast solution of these systems is very urgent nowadays. The size of the problems can be 10 13 unknowns and 10 13 equations. Iterative solution methods are the methods of choice for these large linear systems. We start with a short introduction of Basic Iterative Methods. Thereafter preconditioned Krylov subspace methods, which are state of the art, are describeed. A distinction is made between various classes of matrices. At the end of the lecture notes many references are given to state of the art Scientiﬁc Computing methods. Here, we will discuss a number of books which are nice to use for an overview of background material. First of all the books of Golub and Van Loan [19] and Horn and Johnson [26] are classical works on all aspects of numerical linear algebra. These books also contain most of the material, which is used for direct solvers. Varga [50] is a good starting point to study the theory of basic iterative methods. Krylov subspace methods and multigrid are discussed in Saad [38] and Trottenberg, Oosterlee and Sch¨uller [42]. Other books on Krylov subspace methods are [1,6,21,34,39].


Introduction and model problem
Problems coming from discretized partial differential equations lead in general to large sparse systems of equations Ax = b, where A ∈ R n×n is a non-singular matrix.Direct solution methods, based on an LU -decomposition can be impractical if A is large and sparse.In this case the factors L ∈ R n×n (a lower triangular matrix with ones on the main diagonal) and U ∈ R n×n (an upper triangular matrix) can be dense.This is especially the case for 3D problems.So a considerable amount of memory is required and even the solution of the triangular system costs many floating point operations.
In contrast to the direct methods are the iterative methods.These methods generate a sequence of approximate solutions x (k) and essentially involve the matrix A only in the context of matrix-vector multiplication.The evaluation of an iterative method invariably focuses on how quickly the iterates x (k) converge.The study of round off errors is in general not very well developed.A reason for this is that the iterates are only approximations of the exact solution, so round off errors in general only influence the speed of convergence but not the quality of the final approximation.
The use of iterative solution methods is very attractive in time dependent and nonlinear problems.For these problems the solution of the linear system is part of an outer loop: time stepping for a time dependent problem and Newton Raphson (or other nonlinear methods) for a nonlinear problem.So good start vectors are in general available: the solution of the preceding outer iteration, whereas the required accuracy is in general low for these problems.Both properties lead to the fact that only a small number of iterations is sufficient to obtain the approximate solution of the linear system.Before starting the description and analysis of iterative methods we describe a typical model problem obtained from a discretized partial differential equation.The properties and definitions of the given methods are illustrated using this problem.
In the sequel v ij is an approximation of w(x i , y j ) where x i = ih and y j = jh, 0 ≤ i ≤ m + 1, 0 ≤ j ≤ m + 1.
The finite difference approximation may be written as: The ordering of the nodal points is given in Figure 1 for m = 3.The k th component u k of the vector u is the unknown corresponding to the grid point labeled k.When all the boundary terms are moved to the right-hand side, the system of difference equations can be written as: If the unknowns on lines parallel to the x-axis are grouped together the given system can be written as: where the unknowns on line k are denoted by U k .The lines and grid points are given in Figure 2.1.The matrices A i,j are given by: The submatrix A k,l gives the coupling of the unknowns from line k to those on line l.

Basic iterative methods
The basic idea behind iterative methods for the solution of a linear system Ax = b is: starting from a given x (k) , obtain a better approximation x (k+1) of x in a cheap way.Note that b − Ax (k) is small if x (k) is close to x.This motivates the iteration process One immediately verifies that if this process converges, x is a possible solution.So if b − Ax (k) 2 is large we get a large change of x (k) to x (k+1) .The choice of M is crucial in order to obtain a fast converging iterative method.Rewriting of (1) leads to: where the matrix N is given by N = M − A. The formula A = M − N is also known as a splitting of A. It can easily be seen that if x (k+1) → x the vector x should satisfy As a first example we describe the point Gauss Jacobi method.First we define D A = diag (A) and L A and U A which are the respectively strictly lower and the strictly upper part of A. So that A = D A + L A + U A .
The choice M = D A and N = −(L A + U A ) leads to the point Gauss Jacobi iteration.This algorithm can be described by the following formula: One immediately sees that only memory is required to store the matrix A, the right-hand side vector b and the approximation vector x (k) which can be overwritten in the next step.For our model problem it is sufficient to store 7 vectors in memory, because the matrix A has only 7 non-zero elements per row.Furthermore, one iteration costs approximately as much work as a matrix vector product.
Whether or not the iterates obtained by formula (2) converge to x = A −1 b depends upon the eigenvalues of M −1 N .The set of eigenvalues of A is denoted as the spectrum of A. The spectral radius of a matrix A is defined by ρ(A) = max {|λ|, where λ ∈ spectrum of A} .
The size of ρ(M −1 N ) is critical to the convergence behavior of (2).
Theorem 1.1.Suppose b ∈ IR n and A = M − N ∈ IR n×n is nonsingular.If M is nonsingular and the spectral radius of M −1 N is less than 1, then the iterates x (k) defined by (2) converge to x = A −1 b for any starting vector x (0) .
As an example we note that point Gauss Jacobi is convergent if the matrix A is strictly diagonal dominant.To show this we note that Since a strictly diagonal dominant matrix has the property that In many problems an increase of the diagonal dominance leads to a decrease of the number of iterations.
Summarizing the results given above we see that a study of basic iterative methods proceeds along the following lines: -a splitting A = M − N is given where systems of the form M z = d (d given and z unknown) are cheaply solvable in order to obtain a practical iteration method by (2), -classes of matrices are identified for which the iteration matrix M −1 N satisfies ρ(M −1 N ) < 1.
-further results about ρ(M −1 N ) are established to gain intuition about how the error e (k) tends to zero.
For Gauss Jacobi we note that M is a diagonal matrix so that systems of the form M z = d are easily solvable.We have specified a class of matrices such that convergence occurs.The final point given above is used to obtain new methods with a faster rate of convergence or a wider class of matrices such that ρ(M −1 N ) < 1.
Below we first give a block variant of the Gauss Jacobi method.Thereafter other methods are specified.In the remainder of this section we suppose that Ax = b is partitioned in the form where A i,j is an n i × n j submatrix and n 1 + n 2 . . .+ n q = n.Here the X i and B i represent subvectors of order n i .Furthermore, we define Gauss Jacobi (block) For the given matrices D A , L A and U A the Gauss Jacobi method is given by M = D A and N = −(L A + U A ).The iterates can be obtained by the following algorithm: end for For the special case n i = 1, i = 1, . . ., n we get the point Gauss Jacobi back.In our example a natural block ordering is obtained if we take n i = 3.In that case the diagonal block matrices A i,i are tridiagonal matrices for which cheap methods exist to solve systems of the form A i,i z = d.
Note that in the Gauss Jacobi iteration one does not use the most recently available information when computing is already known.If we revise the Gauss Jacobi iteration so that we always use the most recent estimate of the exact X i then we obtain: Gauss Seidel (block) The Gauss Seidel method is given by M = D A + L A and N = −U A .The iterates can be obtained by

end for
Note that again the solution of systems of the form M z = d are easily obtained since M = D A + L A and L A is a strictly lower triangular block matrix.As an example for the convergence results that can be proved for the point Gauss Seidel method (A i,i = a i,i ), we give the following theorem: Theorem 1.2.If A ∈ IR n×n is symmetric and positive definite, then the point Gauss Seidel iteration converges for any x (0) .Proof: see [19], p. 512, Theorem 10.1.2.This result is frequently applicable because many of the matrices that arise from discretized elliptic partial differential equations are symmetric positive definite.Gauss Seidel converges faster than Gauss Jacobi.Furthermore, block versions converge faster than point versions.For these results and further detailed convergence proofs we refer to [51], [54], [24], and [3].
The Gauss Seidel iteration is in general a slowly converging process.Inspections of the iterates show that the approximations are monotonously increasing or decreasing.Hence, we expect an improvement of the rate of convergence if we use an extrapolation in the direction of the expected change.This leads to the so called successive over-relaxation (SOR) method:

SOR (block)
The successive over-relaxation method is obtained by choosing a constant ω ∈ IR and compute the iterates by The idea behing the use of ω is to damp (if ω > 1) or increase (if ω > 1) the computed adaptation.The following algorithm can be used:

end for
It can be shown that for 0 < ω < 2 the SOR method converges if A is symmetric and positive definite.For ω < 1 we have underrelaxation and for ω > 1 we have overrelaxation.In most examples overrelaxation is used.For a few structured (but important) problems such as our model problem, the value of the relaxation parameter ω that minimizes ρ(M −1 ω N ω ) is known.Moreover, a significant reduction of ρ(M −1 ωopt N ωopt ) with respect to ρ(M −1 1 N 1 ) can be obtained.Note that for ω = 1 we get the Gauss Seidel method.As an example it appears that for the model problem the number of Gauss Seidel iterations is proportional to 1 h 2 whereas the number of SOR iterations with optimal ω is proportional to 1 h .So for small values of h a considerable gain in work results.However, in more complicated problems it may be necessary to perform a fairly sophisticated eigenvalue analysis in order to determine an appropriate ω.A complete survey of "SOR theory" appeared in [54].Some practical schemes for estimating the optimum ω are discussed in [5], [7], and [53].

Chebyshev
The SOR method is presented as an acceleration of the Gauss Seidel method.Another method to accelerate the convergence of an iterative method is the Chebyshev method.Suppose x (1) , . . ., x (k) have been obtained via (2), and we wish to determine coefficients γ j (k), j = 0, . . ., k such that and consider how to choose the γ j (k) so that the error y (k) − x is minimized.It follows from the proof of Theorem 1.1 that e (k) = (M −1 N ) k e (0) where e (k) = x (k) − x.This implies that Using the 2-norm we look for γ j (k) such that y (k) − x 2 is minimal.To simplify this minimization problem we use the following inequality: where p k (z) = k j=0 γ j (k)z j and p k (1) = 1.We now try to minimize p k (M −1 N ) 2 for all polynomials satisfying p k (1) = 1.Another simplification is the assumption that M −1 N is symmetric with eigenvalues λ i that satisfy 0 < α ≤ λ n . . .≤ λ 1 ≤ β < 1.Using these assumptions we see that So to make the norm of p k (M −1 N ) small we need a polynomial p k (z) that is small on [α, β] subject to the constraint that p k (1) = 1.This is a minimization problem of polynomials on the real axis.The solution of this problem is obtained by Chebyshev polynomials.These polynomials c j (z) can be generated by the following recursion These polynomials satisfy |c j (z)| ≤ 1 on [−1, 1] but grow rapidly in magnitude outside this interval.As a consequence the polynomial , and tends to be small on [α, β].The last property can be explained by the fact that numerator is less than 1 in absolute value, whereas the denominator is large in absolute value since 1+2 1−β β−α > 1.This polynomial combined with (8) leads to Calculation of the approximation y (k) by formula (5) costs much time and memory, since all the vectors x (0) , . . ., x (k) should be kept in memory.Furthermore, to calculate y (k) one needs to add k + 1 vectors, which for the model problem costs for k ≥ 5 more work than one matrix vector product.Using the recursion of the Chebyshev polynomials it is possible to derive a three term recurrence among the y (k) .It can be shown that the vectors y (k) can be calculated as follows: solve z (0) from M z (0) = b − Ay (0) then y (1) is given by ) is given by 1) .
We refer to this scheme as the Chebyshev semi-iterative method associated with M y (k+1) = N y (k) + b.Note that only 4 vectors are needed in memory and the extra work per iteration consists of the addition of 4 vectors.
In order for the acceleration to be effective it is necessary to have good lower and upper bounds of α and β.These parameters may be difficult to obtain.Chebyshev semi-iterative methods are extensively analyzed in [51], [20] and [24].
In deriving the Chebyshev acceleration we assumed that the iteration matrix M −1 N is symmetric.Thus our simple analysis does not apply to the SOR iteration matrix M −1 ω N ω because this matrix is not symmetric.To repair this, Symmetric SOR (SSOR) is proposed.In SSOR one SOR step is followed by a backward SOR step.In this backward step the unknowns are updated in reversed order.For further details see [19], Section 10.1.5.Finally we present some theoretical results for the Chebyshev method 1 .Suppose that the matrix M −1 A is symmetric and positive definite and that the eigenvalues µ i are ordered as follows 0 < µ 1 ≤ µ 2 . . .≤ µ n .It is then possible to prove the following theorem: Theorem 1.3.If the Chebyshev method is applied and M −1 A is symmetric positive definite then we see that the eigenvalues satisfy the following relation: This combined with (9) leads to the inequality: So it remains to estimate the denominator.Note that The Chebyshev polynomial can also be given by This expression can be used to show that The condition number K 2 (M −1 A) is equal to µn µ1 .Together with (10) and (11) this leads to Chebyshev type methods which are applicable to a wider range of matrices are given in the literature.In [28] a Chebyshev method is given for matrices with the property that their eigenvalues are contained in an ellipse in the complex plane, and the origin is no element of this ellipse.For a general theory of semi-iterative methods of Chebyshev type we refer to [8].

Starting vectors and termination criteria
Starting vectors All the given iterative solution methods used to solve Ax = b start with a given vector x (0) .In this subsection we shall give some ideas how to choose a good starting vector x (0) .These choices depend on the problem to be solved.If no further information is available one always starts with x (0) = 0.The solution of a nonlinear problem is in general approximated by the solution of a number of linear systems.In such a problem the final solution of the iterative method at a given outer iteration can be used as a starting solution for the iterative method used to solve the next linear system.Suppose we have a time dependent problem.The solution of the continuous problem is denoted by u (n) .In every time step this solution is approximated by a discrete solution u (n) h satisfying the following linear system These systems are approximately solved by an iterative method where the iterates are denoted by x (n,k) .An obvious choice for the starting vector is x (n+1,0) = x (n,kn) where k n denotes the number of iterations in the n th time step.A better initial estimate can be obtained by the following extrapolation: where t .This leads to the following starting vector Finally starting vectors can sometimes be obtained by solution of related problems, e.g., analytic solution of a simplified problem, a solution computed by a coarser grid, a numerical solution obtained by a small change in one of the parameters etc.

Termination criteria
In Subsection 1.2 we have specified iterative methods to solve Ax = b.However, no criteria to stop the iterative process have been given.The iterative method should be stopped if the approximate solution is accurate enough.
To measure the accuracy we can use the norm of the error: x − x (i) .However, to compute the error the exact solution x should be known, which is impossible to have during the iterations.A good termination criterion is very important for an iterative method, because if the criterion is too weak the approximate solution is useless, whereas if the criterion is too severe the iterative solution method never stops or costs too much work.
We start by giving a termination criterion for a linear convergent process.An iterative method is linear convergent if the iterates satisfy the following equation: and 12) is easily checked during the computation.In general, initially (12) is not satisfied but after some iterations the quantity x (i) −x (i−1) 2 x (i−1) −x (i−2) 2 converges to r.The Gauss Jacobi, Gauss Seidel and SOR methods are linear convergent.
Theorem 1.4.For a linear convergent process we have the following inequality Proof Using (12) we obtain the following inequality for k ≥ i + 1.
Since lim The result of this theorem can be used to give a stopping criterion for linear convergent methods.Sometimes the iterations are stopped if x (i) − x (i−1) 2 is small enough.If r is close to one this may lead to inaccurate results since r 1−r x (i) − x (i−1) 2 and thus x − x (i) 2 may be large.A safe stopping criterion is: stop if If this condition holds then the relative error is less than ε: Furthermore, Theorem 1.4 yields the following result: So assuming that the expression (13) can be replaced by an equality log x − x (i) 2 = i log (r) + log This implies that the curve log x − x (i) 2 is a straight line as function of i.This is the motivation for the term linear convergent process.Given the quantity r, which is also known as the rate of convergence, or reduction factor, the required accuracy and x (1) − x (0) 2 it is possible to estimate the number of iterations to achieve this accuracy.In general r may be close to one and hence a small increase of r may lead to a large increase of the required number of iterations.
For iterative methods, which have another convergence behavior, most stopping criteria are based on the norm of the residual.Below we shall give some of these criteria and comment on their properties.
The main disadvantage of this criterion is that it is not scaling invariant.This implies that if b − Ax (i) 2 < this does not hold for 100(b − Ax (i) ) 2 .Although the accuracy of x (i) remains the same.So a correct choice of depends on properties of the matrix A.
The remaining criteria are all scaling invariant.
The number of iterations is independent of the initial estimate x (0) .This may be a drawback since a better initial estimate does not lead to a decrease of the number of iterations.

≤
This is a good termination criterion.The norm of the residual is small with respect to the norm of the right-hand side.Replacing by /K 2 (A) we can show that the relative error in x is less than .It follows that: In general A 2 and A −1 2 are not known.Some iterative methods gives approximations of these quantities.
This criterion is closely related to Criterion 3. In many cases this criterion also implies that the relative error is less than : Sometimes physical relations lead to other termination criteria.This is the case if the residual has a physical meaning, for instance the residual is equal to some energy or the deviation from a divergence free vector field etc.
It is well known that due to rounding errors the solution x represented on a computer has a residual which may be of the magnitude of u b 2 , where u is the machine precision.So we cannot expect that a computed solution by an iterative solution method has a smaller norm of the residual.In a good implementation of an iterative method, a warning is given if the required accuracy is too high.If for instance the termination criterion is b − Ax (i) 2 ≤ and is chosen less than 1000u b 2 a warning should be given and ε should be replaced by = 1000u b 2 .The arbitrary constant 1000 is used for safety reasons.

Exercises
(1) Suppose E < 1 for some matrix E ∈ IR n×n .Show that (2) Show that if A is strictly diagonal dominant then the Gauss Seidel method converges.
(3) Suppose that A is symmetric and positive definite.
(a) Show that one can write Show that P is symmetric.(c) Show that T g can also be written as (e) Show that P = Q T D A Q and P is symmetric and positive definite.(f) Let λ be an eigenvalue of T g with eigenvector x.Use part (b) to show that x T P x > 0 implies that |λ| < 1. (g) Show that the Gauss Seidel method converges.(4) Extend the method of proof in Exercise 3 to the SOR method with 0 < ω < 2.
(5) Suppose that μ1 is an estimate for µ 1 and μn for µ n .
(a) Show that in general the Chebyshev method converges slower if 0 < μ1 < µ 1 and μn > µ n if μ1 and μn are used in the Chebyshev method.(b) Show that divergence can occur if μn < µ n .(6) (a) Do two iterations with Gauss Jacobi to the system: Note that the second iterate is equal to the exact solution.(b) Is the following claim correct?
The Gauss Jacobi method converges in mostly n iterations if A is a lower triangular matrix 2. A Krylov subspace method for systems with a symmetric positive definite matrix

Introduction
In the basic iterative solution methods we compute the iterates by the following recursion: Writing out the first steps of such a process we obtain: This implies that The subspace K i (A; r 0 ) := span r 0 , Ar 0 , . . ., A i−1 r 0 is called the Krylov-space of dimension i corresponding to matrix A and initial residual r 0 .An x i calculated by a basic iterative method is an element of In the preceding chapter we tried to accelerate convergence by the Chebyshev method.In this method one approximates the solution x by a vector is minimized in a certain way.One of the drawbacks of that method is that information on the eigenvalues of M −1 A should be known.In this chapter we shall describe the Conjugate Gradient method.This method minimizes the error x − x i in an adapted norm, without having to know any information about the eigenvalues.In Section 2.3 we give theoretical results concerning the convergence behavior of the CG method.

The Conjugate Gradient (CG) method
In this section we assume that M = I, and x 0 = 0 so r 0 = b.These assumptions are only needed to facilitate the formulas.They are not necessary for the CG method itself.Furthermore, we assume that A satisfies the following condition.

Condition 3.2.1
The matrix A is symmetric (A = A T ) and positive definite (x T Ax > 0 for x = 0).This condition is crucial for the derivation and success of the CG method.Later on we shall derive extensions to non-symmetric matrices.
The first idea could be to construct a vector x i ∈ K i (A, r 0 ) such that x − x i 2 is minimal.The first iterate x 1 can be written as x 1 = α 0 r 0 where α 0 is a constant which has to be chosen such that x − x 1 2 is minimal.This leads to The norm given in ( 15) is minimized if α 0 = r T 0 x r T 0 r0 .Since x is unknown this choice cannot be determined, so this idea does not lead to a useful method.Note that Ax = b is known so using an adapted inner product implying A could lead to an α 0 which is easy to calculate.To follow this idea we define the following inner product and related norm.

Definition 3.2.2
The A-inner product is defined by (y, z) A = y T Az, and the A-norm by y A = (y, y) A = y T Ay.
It is easy to show that if A satisfies Condition 3.2.1 (., .)A and .A satisfy the rules for inner product and norm respectively.In order to obtain x 1 such that x − x 1 A is minimal we note that . We see that this new inner product leads to a minimization problem, which can be easily solved.In the next iterations we compute x i such that The solution of this minimization problem leads to the conjugate gradient method.First we specify the CG method, thereafter we summarize some of its properties.
Conjugate Gradient method The first description of this algorithm is given in [25].For more recent results see [47].Besides the two vectors x k , r k and matrix A only two extra vectors p k and Ap k should be stored in memory.Note that the vectors from the previous iteration can be overwritten.One iteration of CG costs one matrix vector product and 10 n flops for vector calculations.If the CG algorithm is used in a practical application the termination criterion should be replaced by one of the criteria given in Section 1.3.In this algorithm r k is computed from r k−1 by the equation r k = r k−1 − α k Ap k .This is done in order to save one matrix vector product for the original calculation r k = b − Ax k .In some applications the updated residual obtained from the CG algorithm can deviate much from the exact residual b − Ax k due to rounding errors.So it is strongly recommended to recompute b − Ax k after the termination criterion is satisfied for the updated residual and compare the norm of the exact and updated residual.If the exact residual does not satisfy the termination criterion, the CG method should be restarted with x k as its starting vector.
The vectors defined in the CG method have the following properties: Theorem 2.1.
Remarks on the properties given in Theorem 2.1 -It follows from ( 17) and ( 18) that the vectors r 0 , . . ., r k−1 form an orthogonal basis of K k (A; r 0 ).
-In theory the CG method is a finite method.After n iterations the Krylov subspace is identical to IR n .Since x − y A is minimized over K n (A; r 0 ) = IR n the norm is equal to zero and x n = x.However in practice this property is never utilized for two reasons: firstly in many applications n is very large so that it is not feasible to do n iterations, secondly even if n is small, rounding errors can spoil the results such that the properties given in Theorem 2.1 do not hold for the computed vectors.-The sequence x − x k A is monotone decreasing, so This follows from (21) and the fact that K k (A; r 0 ) ⊂ K k+1 (A; r 0 ).In practice x − x k A is not easy to compute since x is unknown.The norm of the residual is given by r This sequence is not necessarily monotone decreasing.In applications it may occur that r k+1 2 is larger than r k 2 .This does not mean that the CG process becomes divergent.The inequality shows that r k 2 is less than the monotone decreasing sequence A 2 x−x k A , so after some iterations the norm of the residual decreases again.
-The search vector p j is A-orthogonal or A-conjugate to all p i with index i less than j.This is the motivation for the name of the method: the search directions or gradients of the updates are mutually conjugate.-In the algorithm we see two ratios, one to calculate β k and the other one for α k .If the denominator is equal to zero, the CG method breaks down.With respect to β k this implies that r T k−2 r k−2 = 0, which implies r k−2 = 0 and thus x k−2 = x.The linear system is solved.The denominator of α k is zero if p T k Ap k = 0 so p k = 0. Using property (17) this implies that r k−1 = 0 so again the problem is already solved.Conclusion: If the matrix A satisfies Condition 3.2.1 then the CG method is robust.
In a following chapter we shall give CG type methods for general matrices A. But first we shall extend Condition 3.2.1 in such a way that also singular matrices are permitted.If the matrix A is symmetric and positive semi definite (x T Ax ≥ 0) the CG method can be used to solve the linear system Ax = b, provided b is an element of the column space of A (range(A)).This is a natural condition because if it does not hold there is no vector x such that Ax = b.For further details and references see [27].

The convergence behavior of the CG method
An important research topic is the rate of convergence of the CG method.The optimality property enables one to obtain easy to calculate upper bounds of the distance between the k th iterate and the exact solution.
Theorem 2.2.The iterates x k obtained from the CG algorithm satisfy the following inequality: Proof We shall only give a sketch of the proof.It is easily seen that x − x k can be written as a polynomial, say p k (A) with p k (0) = 1, times the initial residual (compare to the Chebyshev method) Due to the minimization property every other polynomial q k (A) with q k (0) = 1 does not decrease the error measured in the A-norm: x The right-hand side can be written as Taking q k (A) equal to the Chebyshev polynomial gives the desired result.
Note that the rate of convergence of CG is comparable to that of the Chebyshev method, however it is not necessary to estimate or calculate eigenvalues of the matrix A. Furthermore, increasing diagonal dominance leads to a better rate of convergence.
Initially the CG method was not very popular.The reason for this is that the convergence can be slow for systems where the condition number K 2 (A) is very large.On the other hand the fact that the solution is found after n iteration is also not useful in practice.Firstly n may be very large, secondly the property does not hold in the presence of rounding errors.To illustrate this we consider the following classical example: Example 1 The linear system Ax = b should be solved where n = 40 and b = (1, 0, . . ., 0) T .The matrix A is given by This can be seen as a finite difference discretization of the bending beam equation: u = f .The eigenvalues of this matrix are given by: The matrix A is symmetric positive definite so the CG method can be used to solve the linear system.The condition number of A is approximately equal to 82 π 4 .The resulting rate of convergence given by is close to one.This explains a slow convergence of the CG method for the first iterations.However after 40 iterations the solution should be found.In Figure 2 the convergence behavior is given where the rounding error is equal to 10 −16 , [18].This example suggests that CG has only a restricted range of applicability.These ideas however changed after the publication of [35].Herein it is shown that the CG method can be very useful for a class of linear systems, not as a direct method, but as an iterative method.These problems originate from discretized partial differential equations.It appears that it is not the size of the matrix that is important for convergence but rather the extreme eigenvalues of A.
One of the results which is based on the extreme eigenvalues is given in Theorem 2.2.This inequality is an upper bound for the error of the CG iterates, and suggests that the CG method is a linearly convergent process (see Figure 3).However, in practice the convergence behavior looks like the one given in Figure 4.This is called a superlinear convergence behavior.So the upper bound is only sharp for the initial iterates.It seems that after some iterations the condition number in Theorem 2.2 is replaced by a smaller "effective" condition number.To illustrate this we give the following example: Using Theorem 2.2 it appears that 280 iteration are necessary to ensure that Computing the solution it appears that CG iterates satisfy the given termination criterion after 120 iterations.So in this example the estimate given in Theorem 2.2 is not sharp.
To obtain a better idea of the convergence behavior we have a closer look at the CG method.We have seen that CG minimizes x − x i A over the Krylov subspace.This can also be seen as the construction of a polynomial q i of degree i with q i (0) = 1 such that Suppose that the orthonormal eigensystem of A is given by: {λ j , y j } j=1,...,n where Ay j = λ j y j , λ j ∈ R, y j 2 = 1, y T j y i = 0, j = i, and 0 < λ 1 ≤ λ 2 . . .≤ λ n .The initial errors can be written as x − x 0 = n j=1 γ j y j , which implies that If for instance λ 1 = λ 2 and γ 1 = 0 and γ 2 = 0 it is always possible to change y 1 and y 2 into ỹ1 and ỹ2 such that γ1 = 0 but γ2 = 0.This combined with equation (22) implies that if q i (λ j ) = 0 for all different λ j then x i = x.
So if there are only m < n different eigenvalues the CG method stops at most after m iterations.Furthermore, the upper bound given in Theorem 2.2 can be sharpened.

Remark
For a given linear system Ax = b and a given x 0 (note that x − x 0 = n j=1 γ j y j ) the quantities α and β are defined by: α = min {λ j |γ j = 0} , β = max {λ j |γ j = 0} .It is easy to show that the following inequality holds: The ratio β α is called the effective condition number of A.
It follows from Theorem 2.1 that r 0 , . . ., r k−1 forms an orthogonal basis for K k (A; r 0 ).So the vectors ri = r i / r i 2 , i = 0, . . ., k − 1 form an orthonormal basis for K k (A; r 0 ).We define the following matrices The Ritz matrix T k can be seen as the projection of A on K k (A; r 0 ).It follows from Theorem 2.1 that T k is a tridiagonal symmetric matrix.The coefficients of T k can be calculated from the α i 's and β i 's of the CG process.
The eigenvalues θ i of the matrix T k are called Ritz values of A with respect to K k (A; r 0 ).If z i is an eigenvector of T k so that T k z i = θ i z i and z i 2 = 1 then R k z i is called a Ritz vector of A. Ritz values and Ritz vectors are approximations of eigenvalues and eigenvectors and play an important role in a better understanding of the convergence behavior of CG.The important properties of the Ritz values are: -the rate of convergence of a Ritz value to its limit eigenvalue depends on the distance of this eigenvalue to the rest of the spectrum -in general the extreme Ritz values converge the fastest and their limits are α and β.
In practical experiments we see that, if Ritz values approximate the extreme eigenvalues of A, then the rate of convergence seems to be based on a smaller effective condition number (the extreme eigenvalues seem to be absent).We first give a heuristic explanation.Thereafter an exact result from the literature is cited.
From Theorem 2.1 it follows that r k = A(x − x k ) is perpendicular to the Krylov subspace K k (A; r 0 ).If a Ritz vector is a good approximation of an eigenvector y j of A this eigenvector is nearly contained in the subspace K k (A; r 0 ).These two combined yields that (A(x − x k )) T y j ∼ = 0.The exact solution and the approximation can be written as From (A(x − x k )) T y j = (x − x k ) T λ j y j ∼ = 0 it follows that x T y j ∼ = x T k y j .So the error x − x k has a negligible component in the eigenvector y j .This suggest that λ j does no longer influence the convergence of the CG process.
For a more precise result we define a comparison process.The iterates of this process are comparable to that of the original process, but its condition number is less than that of the original process.Definition Let x i be the i-th iterate of the CG process for Ax = b.For a given integer i let x j denote the j-th iterate of the comparison CG process for this equation, starting with x 0 such that x − x 0 is the projection of x − x i onto span{y 2 , . . ., y n }.
Note that for the comparison process the initial error has no component in the eigenvector y 1 .
Theorem 2.3.[44] Let x i be the i-th iterate of CG, and x j the j-th iterate of the comparison process.Then for any j there holds: , where θ 1 is the smallest Ritz value in the i-th step of the CG process.
The theorem shows that from any stage i on for which θ (i) 1 does not coincide with an eigenvalue λ k , the error reduction in the next j steps is at most the fixed factor F i worse than the error reduction in the first j steps of the comparison process in which the error vector has no y 1 -component.As an example we consider the case that λ 1 < θ (i) 1 < λ 2 we then have , which is a kind of relative convergence measure for θ 1 −λ1 λ2−λ1 < 0.1 then we have F i < 1.25.Hence, already for this modest degree of convergence of θ 1 the process virtually converges as well as the comparison process (as if the y 1 -component was not present).For more general results and experiments we refer to [44].

Exercises
(1) Show that (y, z) A = y T Az is an inner product if A is symmetric and positive definite.
( (6) Suppose that A is nonsingular, symmetric, and indefinite.Give an example to show that the Conjugate Gradients method can break down.

Preconditioning of Krylov subspace methods
We have seen that the convergence behavior of Krylov subspace methods depends strongly on the eigenvalue distribution of the coefficient matrix.A preconditioner is a matrix that transforms the linear system such that the transformed system has the same solution but the transformed coefficient matrix has a more favorable spectrum.As an example we consider a matrix M which resembles the matrix A. The transformed system is given by and has the same solution as the original system Ax = b.The requirements on the matrix M are the following: -the eigenvalues of M −1 A should be clustered around 1, -it should be possible to obtain M −1 y at a low cost.
Most of this chapter contains preconditioners for symmetric positive definite systems (Section 3.1).For nonsymmetric systems the ideas are analogously, so in Section 3.2 we give some details, which can be used only for non-symmetric systems.

The Preconditioned Conjugate Gradient (PCG) method
In Section 2.3 we observed that the rate of convergence of CG depends on the eigenvalues of A. Initially the condition number λn λ1 determines the decrease of the error.After a number of iterations the λn λ1 is replaced by the effective condition number λn−1 λ1 etc.So the question arises, is it possible to change the linear system Ax = b in such a way that the eigenvalue distribution becomes more favorable with respect to the CG convergence?This is indeed possible and the approach is known as: preconditioning of a linear system.Consider the n × n symmetric positive definite linear system Ax = b.The idea behind Preconditioned Conjugate Gradients is to apply the "original" Conjugate Gradient method to the transformed system Ãx = b , where Ã = P −1 AP −T , x = P −T x and b = P −1 b, and P is a nonsingular matrix.The matrix M defined by M = P P T is called the preconditioner.The resulting algorithm can be rewritten in such a way that only quantities without a ˜sign occurs.
Preconditioned Conjugate Gradient method k = 0 ; Observations and properties for this algorithm are: -it can be shown that the residuals and search directions satisfy: -The denominators r T k−2 z k−2 = z T k−2 M z k−2 never vanish for r k−2 = 0 because M is a positive definite matrix.With respect to the matrix P we have the following requirements: -the multiplication of P −T P −1 by a vector should be cheap.(comparable with a matrix vector product using A).Otherwise one iteration of PCG is much more expensive than one iteration of CG and hence preconditioning leads to a costlier algorithm.-The matrix P −1 AP −T should have a favorable distribution of the eigenvalues.It is easy to show that the eigenvalues of P −1 AP −T are the same as for P −T P −1 A and AP −T P −1 .So we can choose one of these matrices to study the spectrum.In order to give more details on the last requirement we note that the iterate x k obtained by PCG satisfies x − So a small condition number of P −1 AP −T leads to fast convergence.Two extreme choices of P show the possibilities of PCG.Choosing P = I we get the original CG method back, whereas if P T P = A the iterate x 1 is equal to x so PCG converges in one iteration.For a classical paper on the success of PCG we refer to [29].
In the following pages some typical preconditioners are discussed.
Diagonal scaling A simple choice for P is a diagonal matrix with diagonal elements p ii = √ a ii .Note that a ii > 0 for i = 1, . . ., n since A is SPD.In [43] it has been shown that this choice minimizes the condition number of P −1 AP −T if P is restricted to be a diagonal matrix.For this preconditioner it is advantageous to apply CG to Ãx = b.The reason is that P −1 AP −T is easily calculated.Furthermore, diag ( Ã) = 1 which saves n multiplications in the matrix vector product.

Basic iterative method
The basic iterative methods described in Section 1.2 use a splitting of the matrix A = M − N .In the beginning of Section 2.2 we show that the k-th iterate y k from a basic method is an element of x 0 + K k (M −1 A, M −1 r 0 ).Using this matrix M in the PCG method we see that the iterate x k obtained by PCG satisfies the following inequality: This implies that x − x k A ≤ x − y k A , so measured in the .A norm, the error of a PCG iterate is less than the error of a corresponding result of a basic iterative method.The extra costs to compute a PCG iterate with respect to the basic iterate are in general negligible.This leads to the notion that any basic iterative method based on the splitting = M − N can be accelerated by the Conjugate Gradient method as long as M (the preconditioner) is symmetric and positive definite.

Incomplete decomposition
This type of preconditioner is a combination of an iterative method and an approximate direct method.As an illustration we use the model problem defined in Section 1.1.The coefficient matrix of this problem A ∈ IR n×n is a matrix with at most 5 nonzero elements per row.Furthermore, the matrix is symmetric and positive definite.
The nonzero diagonals are numbered as follows: m is the number of grid points in the x-direction.
An optimal choice with respect to convergence is to take a lower triangular matrix L such that A = L T L and P = L (L is the Cholesky factor).However it is well known that the zero elements in the band of A become non zero elements in the band of L. So the amount of work to construct L is large.With respect to memory we note that A can be stored in 3n memory positions, whereas L needs m .n memory positions.For large problems the memory requirements are not easily fulfilled.
If the Cholesky factor L is calculated one observes that the absolute value of the elements in the band of L decreases considerably if the "distance" to the non zero elements of A increases.The non zero elements of L on positions where the elements of A are zero are called "fill-in" elements.The observation of the decrease of fill-in motivates to discard fill in elements entirely, which leads to an incomplete Cholesky decomposition of A. Since the Cholesky decomposition is very stable this is possible without break down for a large class of matrices.To specify this in a more precise way we use the following definition:

Definition
The matrix A = (a ij ) is an M -matrix if a ij ≤ 0 for i = j, the inverse A −1 exists and has positive elements (A −1 ) ij ≥ 0.
The matrix of our model problem is an M -matrix.Furthermore, we give a notation for the elements of L which should be kept to zero.The set of all pairs of indices of off-diagonal matrix entries is denoted by The subset Q of Q n are the places (i, j) where L should be zero.Now the following theorem can be proved: Theorem 3.1.If A is a symmetric M -matrix, there exists for each Q ⊂ Q n having the property that (i, j) ∈ Q implies (j, i) ∈ Q, a uniquely defined lower triangular matrix L and a symmetric nonnegative matrix R with l ij = 0 if (i, j) ∈ Q and r ij = 0 if (i, j) ∈ /Q, such that the splitting A = LL T − R leads to a convergent iterative process LL T x i+1 = Rx i + b for every choice x 0 , where Proof (see [29]; p.151.) After the matrix L is constructed it is used in the PCG algorithm.Note that in this algorithm multiplications by L −1 and L −T are necessary.This is never done by forming L −1 or L −T .It is easy to see that L −1 is a full matrix.If for instance one wants to calculate z = L −1 r we compute z by solving the linear system Lz = r.This is cheap since L is a lower triangular matrix so the forward substitution algorithm can be used.

Example
We consider the model problem and compute a slightly adapted incomplete Cholesky decomposition: A = LD −1 L T − R where the elements of the lower triangular matrix L and diagonal matrix D satisfy the following rules: a) l ij = 0 for all (i, j) where If the elements of L are given as follows: it is easy to see that (using the notation as given in (26)) where elements that are not defined should be replaced by zeros.For this example the amount of work for P −T P −1 times a vector is comparable to the work to compute A times a vector.The combination of this incomplete Cholesky decomposition process with Conjugate Gradients is called the ICCG(0) method ( [29]; p. 156).The 0 means that no extra diagonals are used for fill in.Note that this variant is very cheap with respect to memory: only one extra vector to store D is needed.
Another successfull variant is obtained by a smaller set Q.In this variant the matrix L has three more diagonals than the lower triangular part of the original matrix A. This preconditioning is obtained for the choice For the formula's to obtain the decomposition we refer to ( [29]; p. 156).This preconditioner combined with PCG is known as the ICCG(3) method.A drawback is that all the elements of L are different from the eigenvalues are of order O(h 2 ) and the gaps between them are relatively large.So also for ICCG(0) the condition number is O( 1 h 2 ).Faster convergence can be explained by the fact that the constant before 1 h 2 is less for the ICCG(0) preconditioned system than for A and the distribution of the internal eigenvalues is much better so super linear convergence sets in after a small number of iterations.
The success of the ICCG method has led to many variants.In the following we describe two of them MICCG(0) given in [23] (MICCG means Modified ICCG) and RICCG(0) given in [2] (RICCG means Relaxed ICCG).

MICCG
In the MICCG method the MIC preconditioner is constructed by slightly adapted rules.Again A is splitted as follows A = LD −1 L T − R, where L and D satisfy the following rules: a) l ij = 0 for all (i, j) where c) rowsum (LD −1 L T )=rowsum(A) for all rows and (LD −1 L T ) ij = a ij for all (i, j) where a ij = 0 i > j .Using the same notation of L as given in ( 27) we obtain It can be proved that for this preconditioning there is a cluster of small eigenvalues in the vicinity of 1 and the largest eigenvalues are of order 1 h and have large gap ratio's.So the condition number is O(1/h).
In many problems the initial iterations of MICCG(0) converge faster than ICCG(0).Thereafter for both methods super linear convergence sets in.Using MICCG the largest Ritz values are good approximations of the largest eigenvalues of the preconditioned matrix.A drawback of MICCG is that due to rounding errors nonzero components in eigenvectors related to large eigenvalues appear again after some iterations.This deteriorates the rate of convergence.So if many iterations are needed ICCG can be better than MICCG.
In order to combine the advantages of both methods the RIC preconditioner is proposed in [2], which is an average of the IC and MIC preconditioner.For the details we refer to [2].Only the algorithm is given: choose However the question remains: how to choose α?In Figure 8 which is copied from [45] a typical convergence behavior as function of α is given.This motivates the choice α = 0.95, which leads to a very good rate of convergence on a wide range of problems.The above given preconditioners IC, MIC and RIC can be optimized with respect to work.One way to do this is to look at the explicitly preconditioned system: Applying CG to this system one has to solve lower triangular systems of equations with matrix LD −1/2 .The main diagonal of this matrix is equal to D 1/2 .It saves time if we can change this in such a way that the main diagonal is equal to the identity matrix.One idea could be to replace (33) by Note that Lii = 1 for i = 1, ..., n.PCG now is the application of CG to this preconditioned system.

Eisenstat implementation
In this section we restrict ourselves to the IC(0), MIC(0) and RIC(0) preconditioner.We have already noted that the amount of work of one PCG iteration is approximately twice as large as for a CG iteration.In [9] it is shown that much of the extra work can be avoided.If CG is applied to (34) products of the following form are calculated: For the preconditioners used, the off diagonal part of L is equal to the off-diagonal part of Ã.Using this we obtain: So v j+1 can be calculated by a multiplication by L−T and L−1 and some vector operations.The saving in CPU time is the time to calculate the product of A times a vector.Now one iteration of PCG costs approximately the same as one iteration of CG.

General stencils
In practical problems the stencils in finite element methods may be larger or more irregularly distributed than for finite difference methods.The same type of preconditioners can be used.However there are some differences.We restrict ourselves to the IC(0) preconditioner.For the five point stencil we see that the off diagonal part of L is equal to the strictly lower triangular part of A. For general stencils this property does not hold.Drawbacks are: All the elements of L should be stored, so the memory requirements of PCG are twice as large as for CG.Furthermore, the Eisenstat implementation can no longer be used.This motivates another preconditioner constructed by the following rules: ICD (Incomplete Cholesky restricted to the Diagonal).
A is again splitted as A = LD −1 L T − R and L and D satisfy the following rules: a) l ij = 0 for all (i, j) where This enables us to save memory (only the diagonal D should be stored) and CPU time (since now Eisenstat implementation can be used) per iteration.For large problems the rate of convergence of ICD is worse than for IC.Also MICD and RICD preconditioners can be given.

Preconditioning for general matrices
The preconditioning for non-symmetric matrices goes along the same lines as for symmetric matrices.There is a large amount of literature for generalization of the incomplete Cholesky decompositions.In general it is much more difficult to prove that the decomposition does not break down or that the resulting preconditioned system has a spectrum which leads to fast convergence.Since symmetry is no longer important the number of possible preconditioners is much larger.Furthermore, if we have an incomplete LU decomposition of A, we can apply the iterative methods from 4.3.3 to the following three equivalent systems of equations: or The rate of convergence is approximately the same for all variants.When the Eisenstat implementation is applied one should use (37).Otherwise we prefer (38) because then the stopping criterion is based on r 2 = b − Ax k 2 whereas for (36) it is based on U −1 L −1 r k 2 , and for (37) it is based on L −1 r k 2 .

Exercises
(1) Derive the preconditioned CG method using the CG method applied to Ãx = b.
(2) (a) Show that the formula's given in (28)  (c) Prove that the LD −1 L T decomposition (28) exists if a i = a, b i = b, c i = c and A is diagonally dominant.
(4) A practical exercise Use as test matrices: [a, f ] = poisson(30, 30, 0, 0, central ) (a) Adapt the matlab cg algorithm such that preconditioning is included.Use a diagonal preconditioner and compare the number of iterations with cg without preconditioner.(b) Use the formula's given in (28) to obtain an incomplete LD −1 L T decomposition of A. Make a plot of the diagonal elements of D. Can you understand this plot?(c) Use the LD −1 L T preconditioner in the method given in (a) and compare the convergence behavior with that of the diagonal preconditioner.
4. Krylov subspace methods for general matrices

Introduction
In a preceding chapter we discuss the Conjugate Gradient method.This Krylov subspace method can only be used if the coefficient matrix is symmetric and positive definite.In this chapter we discuss Krylov subspace methods for an increasing class of matrices.For these classes we give different iterative methods.At this moment there is no method which is the best for all cases.This is in contrast with the symmetric positive definite case.In Subsection 4.3.4we give some guidelines for choosing an appropriate method for a given system of equations.In Section 4.2 we consider symmetric indefinite systems.General real valued matrices are the subject of Section 4.3.We end this chapter with a section containing iterative methods for complex valued linear systems.

Indefinite symmetric matrices
In this section we relax the condition that A should be positive definite (Chapter 2), and only assume that A is symmetric.This means that x T Ax > 0 for some x and possibly y T Ay < 0 for some y.For the real eigenvalues this implies that A has positive and negative eigenvalues.For this type of matrices .A defines no longer a norm.Furthermore, CG may have a serious break down since p T k Ap k may be zero whereas p k 2 is not zero.In this section we give two different (but related) methods to overcome these difficulties.These methods are defined in [31].

SYMMLQ
In the CG method we have defined the orthogonal matrix R k ∈ IR n×k where the j th column of R k is equal to the normalized residual r j / r j 2 .It appears that the following relation holds where Tk ∈ IR k+1×k is a tridiagonal matrix.This decomposition is also known as the Lanczos algorithm for the tridiagonalisation of A ( [19]; p.477).This decomposition is always possible for symmetric matrices, also for indefinite ones.The CG iterates are computed as follows: Note that T k consists of the first k rows of Tk .If A is positive definite then T k is positive definite.It appears further that in the CG process an LDL T factorization of T k is used to solve (41).This can lead to break down because T k may be indefinite in this section (compare [19]; Section 9.3.1,9.3.2,10.2.6).In the SYMMLQ method [31] problem (41) is solved in a stable way by using an LQ decomposition.So T k is factorized in the following way: with Lk lower triangular.For more details to obtain an efficient code we refer to [31] section 5.

MINRES
In SYMMLQ we have solved (41) in a stable way and obtain the "CG" iteration.However since .A is no longer a correct norm, the optimality properties of CG are lost.To get these back we can try to use the decomposition and minimize the error in the .A T A norm.This is a norm if A is nonsingular.This leads to the following approximation: where Note that AR k y − b 2 = R k+1 Tk y − b 2 using (39).
Starting with x 0 = 0 implies that r 0 = b.Since R k+1 is an orthogonal matrix and b = R k+1 r 0 2 e 1 , where ∈ IR k+1 , we have to solve the following least squares problem Again for more details see [31]; Section 6.7.A comparison of both methods is given in [31].In general the rate of convergence of SYMMLQ or MINRES for indefinite symmetric systems of equations is much worse than of CG for definite systems.Preconditioning techniques for these methods are specified in [36].

Iterative methods for general matrices
In this section we consider iterative methods to solve Ax = b where the only requirement is that A ∈ IR n×n is nonsingular.In the symmetric case we have seen that CG has the following three nice properties: -the approximation x i is an element of K i (A; r 0 ), -optimality, the error is minimal measured in a certain norm, -short recurrences, only the results of one foregoing step is necessary so work and memory do not increase for an increasing number of iterations.
It is shown in [13] that it is impossible to obtain a Krylov method based on K i (A; r 0 ), which has these properties for general matrices.So either the method has an optimality property but long recurrences, or no optimality and short recurrences, or it is not based on K i (A; r 0 ).Some surveys on general iteration methods have been published: [6], [19] Section 10.4, [15], [39], [21], [3].
It appears that there are essentially three different ways to solve non-symmetric linear systems, while maintaining some kind of orthogonality between the residuals: (1) Solve the normal equations A T Ax = A T b with Conjugate Gradients.
(2) Construct a basis for the Krylov subspace by a 3-term bi-orthogonality relation.
(3) Make all the residuals explicitly orthogonal in order to have an orthogonal basis for the Krylov subspace.
These classes form the subject of the following subsections.An introduction and comparison of these classes is given in [30].

CG applied to the normal equations
The first idea to apply CG to the normal equations or AA T y = b with x = A T y (47) is obvious.The first method is known as the CGNR method, whereas the second method is known as the CGNE method.When A is nonsingular A T A is symmetric and positive definite.So all the properties and theoretical results for CG can be used.There are however some drawbacks first the rate of convergence now depends on K 2 (A T A) = K 2 (A) 2 .In many applications K 2 (A) 2 is very large so the convergence of CG applied to ( 46) is very slow.Another difference is that CG applied to Ax = b depends on the eigenvalues of A whereas CG applied to (46) depends on the eigenvalues of A T A, which are equal to the square of the singular values of A.
Per iteration a multiplication with A and A T is necessary, so the amount of work is approximately twice as much as for the CG method.Furthermore, in several (FEM, parallel) applications Av is easily obtained but A T v not, due to the unstructured grid and the data structure used.Finally, not only the convergence depends on K 2 (A) 2 but also the error due to rounding errors.To improve the numerical stability it is suggested in [4] to replace inner products like p T A T Ap by (Ap) T Ap.Another improvement is the method LSQR proposed by [32].This method is based on the application of the Lanczos method to the auxiliary system This is a very reliable algorithm.It uses reliable stopping criteria and estimates of standard errors for x and the condition number of A.

Bi-CG type methods
In this type of methods we have short recurrences but no optimality property.We have seen that CG is based on the Lanczos algorithm.The Lanczos algorithm for non-symmetric matrices is called the bi-Lanczos algorithm.Bi-CG type methods are based on bi-Lanczos.In the Lanczos method we try to find a matrix Q such that Q T Q = I and Q T AQ = T tridiagonal.In the Bi-Lanczos algorithm we construct a similarity transformation X such that To obtain this matrix we construct a basis r 0 , ..., r i−1 , which are the residuals, for K i (A; r 0 ) such that r j ⊥K j (A T ; s 0 ) and s 0 , ..., s i−1 form a basis for K i (A T ; s 0 ) such that s j ⊥K j (A; r 0 ), so the sequences {r i } and {s i } are biorthogonal.Using these properties and the definitions R k = [r 0 ...r k−1 ], S k = [s 0 ...s k−1 ] the following relation can be proved [46]: Since S T k R k is a diagonal matrix with diagonal elements r T j s j , we find, that if all these diagonal elements are nonzero, We see that this algorithm fails when a diagonal element of S T k R k becomes (nearly) zero, because these elements are used to normalize the vectors s j (compare [19] §9.3.6).This is called a serious (near) break down.The way to get around this difficulty is the so-called look-ahead strategy.For details on look-ahead we refer to [33], and [16].Another way to avoid break down is to restart as soon as a diagonal element gets small.This strategy is very simple, but one should realize that at a restart the Krylov subspace that has been built up so far, is thrown away, which destroys possibilities for faster (super linear) convergence.(The description of the methods given below is based on those given in [46].)Bi-CG As has been shown for Conjugate Gradients, the LU decomposition of the tridiagonal system T k can be updated from iteration to iteration and this leads to a recursive update of the solution vector.This avoids having to save all intermediate r and s vectors.This variant of Bi-Lanczos is usually called Bi-Conjugate Gradients, or shortly Bi-CG [14].Of course one can in general not be sure that an LU decomposition (without pivoting) of the tridiagonal matrix T k exists, and if it does not exist then a serious break down of the Bi-CG algorithm occurs.This break down can be avoided in the Bi-Lanczos formulation of this iterative solution scheme.The algorithm is given as follows: Bi-CG x 0 is given; r 0 = b − Ax 0 ; r0 is an arbitrary vector (r 0 , r 0 ) = 0 possible choice r0 = r 0 ; ρ 0 = 1 p0 = p 0 = 0 for i = 1, 2, ... ρ i = (r i−1 , r i−1 ) ; Note that for symmetric matrices, Bi-Lanczos generates the same solution as Lanczos, provided that s 0 = r 0 , and under the same condition, Bi-CG delivers the same iterates as CG, for symmetric positive definite matrices.However, the Bi-orthogonal variants do so at the cost of two matrix vector operations per iteration.

QMR
The QMR method [17] relates to Bi-CG in a similar way as MINRES relates to CG.For stability reasons the basis vectors r j and s j are normalized (as is usual in the underlying Bi-Lanczos algorithm).
If we group the residual vectors r j , for j = 0, ..., i − 1 in a matrix R i , then we can write the recurrence relations as Similar as for MINRES we would like to construct the x i , with for which is minimal.However, in this case that would be quite an amount of work since the columns of R i+1 are not necessarily orthogonal.In [17] it is suggested to solve the minimum norm least squares problem min y∈I R i Ti y − r 0 2 e 1 2 . ( This leads to the simplest form of the QMR method.A more general form arises if the least squares problem ( 50) is replaced by a weighted least squares problem.No strategies are yet known for optimal weights, however.
In [17] the QMR method is carried out on top of a look-ahead variant of the bi-orthogonal Lanczos method, which makes the method more robust.Experiments suggest that QMR has a much smoother convergence behavior than Bi-CG, but it is not essentially faster than Bi-CG.For the algorithm we refer to [3] page 24.

CGS
For the bi-conjugate gradient residual vectors it is well-known that they can be written as r j = P j (A)r 0 and rj = P j (A T )r 0 , where P j is a polynomial of degree j such that P j (0) = 1.From the bi-orthogonality relation we have that (r j , ri ) = (P j (A)r 0 , P i (A T )r 0 ) = (P i (A)P j (A)r 0 , r0 ) = 0 , for i < j.
The iteration parameters for bi-conjugate gradients are computed from innerproducts like above.Sonneveld observed in [41] that we can also construct the vectors r j = P 2 j (A)r 0 , using only the latter form of the innerproduct for recovering the bi-conjugate gradients parameters (which implicitly define the polynomial P j ).By doing so, it can be avoided that the vectors rj have to be formed, nor is there any multiplication with the matrix A T .The resulting CGS [41] method works in general very well for many non-symmetric linear problems.It converges often faster than Bi-CG (about twice as fast in some cases).However, CGS usually shows a very irregular convergence behavior.This behavior can even lead to cancellation and a spoiled solution [45].
The following scheme carries out the CGS process for the solution of Ax = b, with a given preconditioner K:

Conjugate Gradient Squared method
x 0 is an initial guess; r 0 = b − Ax 0 ; r0 is an arbitrary vector, such that (r 0 , r0 ) = 0 , e.g., r0 = r 0 ; ρ 0 = (r 0 , r0 ) ; In exact arithmetic, the α j and β j are the same constants as those generated by Bi-CG.Therefore, they can be used to compute the Petrov-Galerkin approximations for eigenvalues of A.
Bi-CGSTAB Bi-CGSTAB [46] is based on the following observation.Instead of squaring the Bi-CG polynomial, we can construct other iteration methods, by which x i are generated so that r i = Pi (A)P i (A)r 0 with other i th degree polynomials Pi .An obvious possibility is to take for Pi a polynomial of the form and to select suitable constants ω j ∈ IR.This expression leads to an almost trivial recurrence relation for the Q i .In Bi-CGSTAB ω j in the j th iteration step is chosen so as to minimize r j , with respect to ω j , for residuals that can be written as r j = Q j (A)P j (A)r 0 .The preconditioned Bi-CGSTAB algorithm for solving the linear system Ax = b, with preconditioning K reads as follows: Bi-CGSTAB method x 0 is an initial guess; r 0 = b − Ax 0 ; r0 is an arbitrary vector, such that (r 0 , r 0 ) = 0, e.g., r0 = r 0 ; ρ The matrix K in this scheme represents the preconditioning matrix and the way of preconditioning [46].The above scheme in fact carries out the Bi-CGSTAB procedure for the explicitly postconditioned linear system AK −1 y = b , but the vector y i has been transformed back to the vector x i corresponding to the original system Ax = b.Compared to CGS two extra innerproducts need to be calculated per iteration.
In exact arithmetic, the α j and β j have the same values as those generated by Bi-CG and CGS.Hence, they can be used to extract eigenvalue approximations for the eigenvalues of A (see Bi-CG).
An advantage of these methods is that they use short recurrences.A disadvantage is that there is only a semioptimality property.As a result of this, more matrix vector products are needed and no convergence properties have been proved.In experiments we see that the convergence behavior looks like CG for a large class of problems.However, the influence of rounding errors is much more important than for CG.Small changes in the algorithm can lead to instabilities.Finally, it is always necessary to compare the norm of the updated residual to the exact residual b − Ax k 2 .If "near" break down had occurred, these quantities may be different by several orders of magnitude.In such a case the method should be restarted.

GMRES-type methods
These methods are based on long recurrences, and have certain optimality properties.The long recurrences imply that the amount of work per iteration and required memory grow for increasing number of iterations.Consequently in practice one cannot afford to run the full algorithm, and it becomes necessary to use restarts or to truncate vector recursions.In this section we describe GMRES, GCR and a combination of both GMRESR.
In GMRES (General Minimal RESidual method) the approximate solution As a consequence of (51) it appears that r k is orthogonal to AK k (A; r 0 ), so r k ⊥K k (A; Ar 0 ).If A is symmetric the GMRES method is equivalent to the MINRES method as described in [31].Using the matrix Hk it follows that AV k = V k+1 Hk where the n × k matrix V k is defined by With this equation it is shown in [40] that x k = x 0 + V k y k where y k is the solution of the following least squares problem: with β = r 0 2 and e 1 is the first unit vector in IR k+1 .GMRES is a stable method and no break down occurs, if h j+1,j = 0 than x j = x so this is a "lucky" break down (see [40]; Section 3.4).
Due to the optimality (see inequality (51)) convergence proofs are possible [40].If the eigenvalues of A are real the same bounds on the norm of the residual can be proved as for the CG method.For a more general eigenvalue distribution we give a result in the following theorem.Let P m be the space of all polynomials of degree less than m and let σ represent the spectrum of A.
Theorem 4.1.Suppose that A is diagonalizable so that A = XDX −1 , σ = {λ 1 , ..., λ n }, and let Then the residual norm of the m-th iterate satisfies: where K(X) = X 2 X −1 2 .If furthermore all eigenvalues are enclosed in a circle centered at C ∈ IR with C > 0 and having radius R with C > R, then Proof: see [40]; p. 866.
For GMRES we see in many cases a super linear convergence behavior comparable to CG.The same type of results are proved for GMRES [48].As we have already noted in the beginning, work per iteration and memory requirements increase for an increasing number of iterations.In this algorithm the Arnoldi process requires k vectors in memory in the k-th iteration.Furthermore, 2k 2 •n flops are needed for the total Gram Schmidt process.
To restrict work and memory requirements one stops GMRES after m iterations, forms the approximate solution and uses this as a starting vector for a following application of GMRES.This is denoted by the GMRES(m) procedure.Not restarted GMRES is denoted by full GMRES.However restarting destroys many of the nice properties of full GMRES, for instance the optimality property is only valid inside a GMRES(m) step and the super linear convergence behavior is lost.This is a severe drawback of the GMRES(m) method [11,22].

GCR
Slightly earlier than GMRES, [10] proposed the GCR method (Generalized Conjugate Residual method).The algorithm is given as follows: GCR algorithm choose x 0 , compute r 0 = b − Ax 0 for i = 1, 2, ... do s i = r i−1 , v i = As i , for j = 1, ..., i − 1 do α = (v j , v i ) , s i := s i − αs j , v i := v i − αv j , end for s i := s i / v i 2 , v i := v i / v i 2 x i := x i−1 + (v i , r i−1 )s i ; r i := r i−1 − (v i , r i−1 )v i ; end for The storage of s i and v i costs two times as much memory as for GMRES.The rate of convergence of GCR and GMRES are comparable.However, there are examples where GCR breaks down.So comparing full GMRES and full GCR the first one is preferred in many applications.When the required memory is not available GCR can be restarted.Furthermore, another strategy is possible which is known as truncation.An example of this is to replace the j-loop by for j = i − m, ..., i − 1 do Now 2m vectors are needed in memory.Other truncation variants to discard search direction are possible.In general we see that truncated methods have a better convergence behavior especially if super linear convergence plays an important role.So if restarting or truncation is necessary truncated GCR is in general better than restarted GMRES.For convergence results and other properties we refer to [10].

GMRESR
A number of methods are proposed to diminish the disadvantages of restarting and or truncation.One of these methods is GMRESR proposed in [49] and further investigated in [52].This method consists of an outer-and an inner loop.In the inner loop we approximate the solution of a linear system with GMRES to find a good search direction.Thereafter in the outer loop the minimal residual approximation using these search directions is calculated by a GCR approach.GMRESR algorithm choose x 0 and m, compute r 0 = b − Ax 0 for i = 1, 2, ... do s i = P m,i−1 (A)r i−1 , v i = As i , for j = 1, ..., i − 1 do α = (v j , v i ) , s i := s i − αs j , v i := v i − αv j , end for s i := s i / v i 2 , v i := v i / v i 2 x i := x i−1 + (v i , r i−1 )s i ; r i := r i−1 − (v i , r i−1 )v i ; end for The notation s i = P m,i−1 (A)r i−1 denotes that one applies one iteration of GMRES(m) to the system As = r i−1 .The result of this operation is s i .For m = 0 we get GCR, whereas for m → ∞ one outer iteration is sufficient and we get GMRES.For the amount of work we refer to [49], where also optimal choices of m are given.In many problems the rate of convergence of GMRESR is comparable to full GMRES, whereas the amount of work and memory is much less.In the following picture we visualize the strong point of GMRESR in comparison with GMRES(m) and truncated GCR.A search direction is indicated by the symbol v i .We see for GMRES(3) that after 3 iterations all information is thrown away.For GCR(3) a window of the last 3 vectors moves from left to right.For GMRESR the information after 3 inner iterations is condensed into one search direction so "no" information gets lost.Also for GMRESR restarts and truncation are possible [52].In the inner loop other iterative methods can be used.Several of these choices lead to a good iterative method.In theory we can call the same loop again, which motivates the name GMRES Recursive.A comparable approach is the FGMRES method give in [37].Herein the outer loop consists of a slightly adapted GMRES algorithm.Since FGMRES and GMRESR are comparable in work and memory, but FGMRES can not be truncated, we prefer the GMRESR method.

Choice of an iterative method
For non-symmetric matrices it is very difficult to decide which iterative method should be used.All the methods treated here have their own type of problems for which they are winners.Furthermore, the choice depends on the computer used and the availability of memory.In general CGS and Bi-CGSTAB are easy to implement and reasonably fast for a large class of problems.If break down or bad convergence appear, GMRES like methods are better.Finally, LSQR always converges but can take a large number of iterations.
In [52] we have tried to specify some easy to obtain parameters to facilitate a choice.First one should have a crude idea of the total number of iterations (mg) using full GMRES, secondly one should measure the ratio f which is defined as f = the CPU time used for one preconditioned matrix vector product the CPU time used for a vector update .
Note that f depends on the used computer.Under certain assumptions given in [52] we obtain Figure 10.This figure gives only qualitative information.It illustrates the dependence of the choice on f and mg.If mg is large and f is small, Bi-CGSTAB is the best method.For large values of f and small values of mg the GMRES method is optimal and for intermediate values GMRESR is the best method.In [3] a flowchart is given with suggestions for the selection of a suitable iterative method.

Iterative methods for complex matrices
There are in practice, important applications that lead to linear systems where the coefficient matrix has complex valued entries.Examples are: complex Helmholtz equations, Schrödinger's equation, under water acoustics etc [12].If the resulting system is Hermitian the methods of Chapter 2 can be used.In these algorithms the inner product x T y should be replaced by the complex inner product xT y.For non Hermitian matrices iterative methods as given in Sections 4.3.1 to 4.3.3can be used.Again they should be adapted to use the correct inner product.In many applications the resulting complex linear systems have additional structure that can be exploited.For instance matrices of the following form arise: Another special case that arises frequently in applications are complex symmetric matrices A = A T .
For example, the complex Helmholtz equations leads to complex symmetric systems.For methods to solve these systems we refer to [16]; section 2.2, 2.3, 6. (3) In the GCR algorithm the vector r i is obtained from vectorupdates.Show that the relation r i = b − Ax i is valid.(4) Prove the following properties for the GMRES method:

2 Figure 1 .
Figure 1.The natural (lexicographic) ordering of the nodal points

Figure 2 .
Figure 2. The convergence behavior of CG applied to Example 1.


are going to solve Ax = b.(a)Show that Conjugate Gradients applied to this system should convergence in 1 or 2 iterations (using the convergence theory).(b)Choose x (0) and do 2 iterations with the Conjugate Gradients method.
means that if Ax = b and x and/or b are slowly varying vectors this incomplete Cholesky decomposition is a very good approximation for the inverse of A with respect to x and/or b.

2 Figure 7 .
Figure 7. Effect of number of equations on the rate of convergence the average parameter α ∈ [0, 1] then di , bi and ci are given by:

Figure 8 .
Figure 8. Convergence in relation with α are correct.(b) Show that the formula's given in (31) are correct.(3) (a) Suppose that a i = 4 and b i = −1.Show that lim i→∞ di = 2 + √ 3, where di is as defined in (28).(b) Do the same for a i = 4, b i = −1 and c i = −1 with m = 10, and show that lim i→∞ di = 2 + √ 2.
and S T k (Ax k − b) = 0 .Using (48), r 0 = b and x k = R k y we obtain

( 1 ) 1 −1 2 .
Show that the solution y x satisfies A T Ax = A T b. (2) Take the following matrix 2 −(a) Suppose that GCR is applied to the system Ax = b.Show that GCR converges in 1 iteration if x − x 0 = cr 0 , where c = 0 is a scalar and r 0 = b − Ax 0 .(b) Apply GCR for the choices b = Do the same for x 0 = 1 0 .