# Milestones in the Development of Iterative Solution Methods

Authors: Owe Axelsson
Journal: Journal of Electrical and Computer Engineering
 (2010)

---

## Abstract

Iterative solution methods to solve linear systems of equations were originally formulated as basic iteration methods of defect-correction type, commonly referred to as Richardson's iteration method. These methods developed further into various versions of splitting methods, including the successive overrelaxation (SOR) method. Later, immensely important developments included convergence acceleration methods, such as the Chebyshev and conjugate gradient iteration methods and preconditioning methods of various forms. A major strive has been to find methods with a total computational complexity of optimal order, that is, proportional to the degrees of freedom involved in the equation. Methods that have turned out to have been particularly important for the further developments of linear equation solvers are surveyed. Some of them are presented in greater detail.


---

## Body

## 1. Introduction

In many applications of quite different types appearing in various sciences, engineering, and finance, large-scale linear algebraic systems of equations arise. A particular type of problems appear in signal processing. This also includes nonlinear systems of equation, which are normally solved by linearization at each outer nonlinear iteration step, but they will not be further discussed in this paper.Due to their high demand of computer memory and computer time, which can grow rapidly with increasing problem size, direct solution methods, such as Gaussian elimination, are in general not feasible unless the size of the problem is relatively small. In the early computer age, when available size of computer central memories was very small and the speed of arithmetic operations slow, this was found to be the case even for quite modest-sized problems.Even for modern computers with exceedingly large memories and very fast arithmetics it is still an issue because nowadays one wants to solve much more involved problems of much larger sizes, for instance to enable a sufficient resolution of partial differential equation problems with highly varying (material) coefficients, such as is found in heterogeneous media. Presently problems with up to billions of degrees of freedom (d.o.f.) are solved. For instance, if an elliptic equation of elasticity type is discretized and solved on a 3D mesh with 512 meshpoints in each coordinate direction, then an equation of that size arises.A basic iteration method to solve a linear system(1)Ax=b,
where A is nonsingular, has the following form.Given an initial approximationx0, for k=0,1,… until convergence, let rk=Axk-b,ek=-τrk,  xk+1=xk+ek.  Here, τ>0 is a parameter to be chosen.This method can be described either as a defect(rk)—correction (ek)method or, alternatively, as a method to compute the stationary solution of the evolution equation(2)dx(t)dt=Ax(t)-b,t>0,x(0)=x0,
by timestepping with time-step τ, that is,(3)x(t+τ)=x(t)-τ(Ax(t)-b),t=0,τ,….Such methods are commonly referred to as Richardson iteration methods (e.g., see [1–4]). However, already in 1823 Gauss [5] wrote, “Fast jeden Abend mache ich eine neue Auflage des Tableau, wo immer leicht nachzuhelfen ist. Bei der Einförmigkeit des Messungsgeschäfts gibt dies immer eine angenehme Unterhaltung; man sieht daran auch immer gleich, ob etwas Zweifelhaftes eingeschlichen ist, was noch wünschenswert bleibt usw. Ich empfehle Ihnen diesen Modus zur Nachahmung. Schwerlich werden Sie je wieder direct eliminieren, wenigstens nicht, wenn Sie mehr als zwei Unbekannte haben. Das indirecte Verfahren läßt sich halb im Schlafe ausführen oder man kann während desselben an andere Dingen denken.” (Freely translated, “I recommend this modus operandi. You will hardly eliminate directly anymore, at least not when you have more than two unknowns. The indirect method can be pursued while half asleep or while thinking about other things.”)It holds that(4)xk+1=xk-τ(Axk-b),
or(5)ek+1=(I-τA)ek,
where ek=x-xkis the iteration error and x is the solution of (1).Hence,(6)ek=(I-τA)ke0,k=0,1,….For convergence of the method, that isek→0, the parameter τ must in general be chosen such that ρ:=∥I-τA∥<1, where ∥·∥is a matrix norm, subordinate to the chosen vector norm. (We remark here that this is not possible if Ais indefinite.)Letρ(·)denote the spectral radius of a matrix, that is, the maximal absolute value of the eigenvalues of the matrix.IfAis self–adjoint, then it can be shown that ρ(A)=∥A∥2=ρ(A*A), where ∥·∥2denotes the matrix norm subordinate to the Euclidian vector norm.  For general, nonsymmetric matrices it has been shown (e.g., see [6, page 162]) that there exist matrix norms that are arbitrarily close to the spectral radius. These can, however, correspond to an unnatural scaling of the matrix.The rate of convergence is determined by the convergence factorρ. For symmetric positive definite matrices, the optimal value of τto minimize ρ, is τ=2/(λ1+λn), where λ1,  λnare the extreme eigenvalues of A. Normally, however, the eigenvalues are not available.As an example, for second order elliptic diffusion type of problems inΩd(d=2,3)using a standard central difference or a finite element method, the spectral condition number λn/λ1=O(h-2), where his the (constant) meshsize parameter. Hence, the number of iterations to reach a relative accuracy  ɛ is of order O(h-2)|logɛ|),h→0.Since each iteration usesO(h-d)elementary arithmetic operations, this shows that the total number of operations needed to reduce the error to a given tolerance is of order O(h-d-2). This is in general smaller than for a direct solution method when d≥2, but still far from the optimal order, O(h-d), that we aim at.To improve on this, often a splitting of the matrixA is used. It is readily shown that for any initial vector, the number of iterations required to get a relative residual, ∥rk∥/∥r0∥<ɛ, for some ɛ,0<ɛ<1, is at most kit=⌈ln(1/ɛ)/ln(1/ρ)+1⌉, where ⌈⌉ denotes the integer part. Frequently, ρ=1-cδr, where c is a constant, ris a positive integer, often r=2and δis a small number, typically δ=1/n, which decreases with increasing problems size n. This implies, that the number of iterations is propotional to (1/δ)r, which number increases rapidly when δ→0.Forτ=1, the splitting A=C-R of A in two terms is used, where C is nonsingular. The iterative method (4) then takes the form(7)Cxk+1=Rxk+b,k=0,1,….
Method (7) is convergent if ρ(C-1R)<1. Splitting methods will be discussed in Section 2.LetB=C-1R. If ∥B∥ is known and ∥B∥<1, we can use the following theorem to get a test when the iteration error is small enough, that is, when to stop the iterations.Theorem 1.
Let∥B∥<1,  B=C-1R, and xm be defined by (7). Then,
(8)‖x-xm‖≤‖B‖1-‖B‖‖xm-xm-1‖,m=1,2,….Proof.
From (7) follows xm+1-xm=B(xm-xm-1) and, by recursion,(9)xm+k+1-xm+k=Bk+1(xm-xm-1),k=0,1,….
Note now that xm+p-xm=∑k=0p-1(xm+k+1-xm+k). Hence, by the triangle inequality and (9)
(10)‖xm+p-xm‖≤∑k=0p-1‖Bk+1‖‖xm-xm-1‖≤‖B‖-‖B‖p+11-‖B‖‖xm-xm-1‖.
Letting p→∞ and noting that xm+p→x̂, (8) follows.The basic iteration method (4) or the splitting methods in Section 2, can be improved in various ways. This will be the major topic of this paper.Note first that application of the splitting in (7) requires in general that the matrix R is given in explicit form, which can make the method less viable.The most natural way to improve (4) is to introduce an approximation C of A, to be used when the correction ek in (4) is computed. The relation ek=-τrk is then replaced by Cek=-τrk.Such a matrix is mostly called preconditioner since, by a proper choice, it can significantly improve the condition number𝒦 of A, that is,(11)K(C-1A)≪K(A),
where 𝒦(B)=∥B∥∥B-1∥.Clearly, in practice, the matrixC must be chosen such that the linear systems with C can be solved with relatively little expense compared to a solution method for A. In particular, the expense for C is much smaller than that for A using a direct solution method.For badly scaled matricesA a simple, but often practically useful, choice of C is the (block) diagonal part D of A. Much more efficient choices will be discussed later in the paper.Early suggestions to use such a matrixC can be found in papers by D'Yakonov [7] and Gunn [8]. For an automatic scaling procedure, see [9] and references therein.In the present paper, we will survey various choices ofC which have proven to be useful in practice. The paper attempts to give a more personal account of the development of iterative solution methods. It is also not our ambition to present the present state- of- the- art but rather to describe the unfolding of the field.In the remainder of the paper, we discuss, in order, methods based on splitting of the given matrix, the accelerated iterative methods of the Chebyshev and (generalized) conjugate gradient types, pointwise and block incomplete factorization preconditioning methods, symmetrized preconditioners of SSOR and ADI type, approximate inverses, and augmented subspace preconditioning methods. If space had allowed it, it would have been followed by presentation of geometric and algebraic multigrid methods, two-level and multilevel methods, elementwise preconditioning methods, and domain decomposition methods. Also, iteration error estimates and influence of rounding errors, and preconditioners for matrices of saddle point type would have been included. The paper ends with some concluding remarks.The following relations will be used; ifA=[aij], then AT=[aji] denotes the transpose of A, while A*=[a¯ji], denotes the Hermitian transpose.


## 2. Splitting Methods

A comprehensive, early presentation of splitting methods, and much more on iterative solution methods, is found in Varga [10]. Often there is a natural splitting of the given matrix as(12)A=C-R,whereC is nonsingular. This can be used to formulate an iterative solution method in the form(13)Cxk+1=Rxk+b,k=0,1,….This method converges ifρ(C-1R)<1.Definition 1.
(a)A matrix C is said to be monotone if C is nonsingular and C-1≥0 (componentwise).
(b)A=C-R is called a regular splitting [10], if C is monotone and R≥0.
(c) aweak regular splitting [11], if C is monotone and C-1R≥0.
(d) anonnegative splitting [12], if C is nonsingular and C-1R≥0.
The following holds, see, for example, [6].Theorem 2.
LetA=C-R be a nonnegative splitting of A. Then, the following properties are equaivalent: (a)
ρ(B)<1, that is, A=C-R is a convergent splitting,(b)
I-B is monotone,(c)
A is nonsingular and G=A-1R≥0.(d)
A is nonsingular and ρ(B)=ρ(G)/[1+ρ(G)], where G=A-1R.Corollary 1.
IfA=C-R is a weak regular splitting, then the splitting is convergent if and only if A is monotone.Proof.
(see, e.g., [6]).A splitting method that became popular in the fifties is the SOR method. Here,A=D-L-U, where D is the (block) diagonal and L,  U are the (block) lower and upper triangular parts of A, respectively. The successive relaxation method takes the form(14)(1ωD-L)xk+1=[(1ω-1)D+U]xk+b,k=0,1,…,
where ω≠0 is a parameter, called the relaxation parameter. For ω=1 one gets the familiar Gauss-Seidel method (Gauss 1823 [5] and Seidel 1814 [13]) and for ω>1 the successive overrelaxation (SOR) method (Frankel 1950 [14] and Young 1950 [15]).For the iteration matrix in (14),(15)Lω=(1ωD-L)-1((1ω-1)D+U),it holds thatρ(ℒω)≤|ω-1|, where the upper bound is sharp. Therefore, the relaxation method is divergent for ω≤0 and ω≥2 (see, e.g., [6, 16]).An optimal value ofω can be determined as follows. Assume then that A has property (Aπ), that is, there exists a permutation matrix P such that PAPT is a block tridiagonal matrix. The following Lemma holds.Lemma 1 (see [15]).
Assume thatA has property (Aπ) and let ω≠0. Let B:=D-1(L+U). Then,(a)
ifλ≠0 is an eigenvalue of ℒω and μ satisfies
(16)μ2=(λ+ω-1)2(ω2λ),
then, μ is an eigenvalue of B,(b)
ifμ is an eigenvalue of B and λ satisfies
(17)λ+ω-1=ωμλ1/2,
then, λ is an eigenvalue of ℒω.Proof.
For a short proof, see [6].Theorem 3.
Assume that(a)
A has property (Aπ), and(b)
the block matrixB=I-D-1A has only real eigenvalues. Then, the SOR method converges for any initial vector if and only if ρ(B)<1 and 0<ω<2. Further, we have
(18)ωopt=21+1-ρ(B)2,
for which the asymptotic convergence factor is given as
(19)min⁡ωρ(Lω)=ρ(Lωopt)=ωopt-1=(1-1-ρ(B)2)(1+1-ρ(B)2).Proof.
Fort a short proof, see [6].The eigenvalues ofC-1A are in general complex, and for ω=ωopt it can be shown that they are distributed around a circle in the complex plane. This implies that the method can not be polynomially accelerated. (See Section 3 for a presentation of polynomial acceleration methods.) Further, the efficiency of the SOR method turns out to be critically dependent on the choice of ω.A similar result as in Theorem3 has been shown in [6], see also [17], that holds even if A does not have property (Aπ), but is Hermitian.Theorem 4.
LetA be Hermitian and positive definite and let
(20)L̃ω=D1/2LωD-1/2=(1ωI-L̃)-1((1ω-1)I-L̃*),
where L̃=D-1/2LD1/2, and let 0<ω<2. Then,
(21)ρ(Lω)2=ρ(L̃ω)2≤1-2/ω-1(1/ω-1/2)2δ-1+γ+1/ω,
where
(22)γ=sup⁡x≠0{[(|x,L̃x|2/(x,x))-1/4(x,x)](Ãx,x)},δ=λmin⁡(Ã)=min⁡x≠0(Ãx,x)(x,x).
Further, if |(x,L̃x)|≤1/2(x,x), then
(23)ω*=21+2δ,
minimizes the upper bound of ρ(Lω) and we have
(24)ρ(Lω*)2=1-δ/21+δ/2.For a proof, see [6].In Section4, we present a symmetric version of the SOR method where acceleration is possible.


## 3. Accelerated Iterative Methods

In this section, the important Chebyshev and conjugate gradient iteration methods are presented.Consider first the iterative method (4) with variable time-steps τk,(25)xk+1=xk-τkC-1rk,rk=Axk-b,k=0,1,….
Here, {τk} is a sequence of iteration (acceleration) parameters. If τk=τ,  k=0,1,…,  we talk about a stationary iterative method, otherwise about a nonstationary or semiiterative method.Letek=x-xk,   the iteration error. Then, it follows from (25) that ek+1=(I-τkC-1A)ek,  k=0,1,…, so   em=Pm(C-1A)e0 (and  rm=APm(C-1A)A-1r0=Pm(AC-1)r0).  Here,  Pm(λ)=Πk=0m(1-τkλ) a polynomial of degree m having zeros at 1/τk and satisfying Pm(0)=1.We want to choose the parameters{τk} such that ∥em∥ is minimized. However, this would mean that in general the parameters would depend on e0,  which is not known. Also the eigenvalues of C-1A are not known. We then take the approach of minimizing ∥em∥/∥e0∥ for all e0; that is, we want to minimize ∥Pm(C-1A)r0∥.


### 3.1. The Chebyshev Iterative Method

In case the eigenvalues ofC-1A are real and positive and if a positive lower (a) and (b) an upper bound are known of the spectrum, then we see that {τk} should be chosen such that maxa≤λ≤b|Pm(λ)| is minimized over all Pm∈Πm0, that is, over the set of polynomials of degree m satisfying Pm(0)=1.The solution to this min-max problem is well known,(26)Pm(λ)=Tm((b+a-2λ)/(b-a))Tm((b+a)/(b-a)),where Tm(z)=(1/2)[(z+(z2-1)1/2)m+(z-(z2-1)1/2)m]=cos(marccosz) are the Chebyshev polynomials of the first kind. The corresponding values of τk satisfy(27)1τk=b-a2cos⁡Θk+b+a2,Θk=2k-12mπ,k=1,2,…,m,which are the zeros of the polynomial. The corresponding method is referred to as the Chebyshev (one-step) acceleration method, see, for example, [10, 18]. It is an easy matter to show that(28)1Tm(b+ab-a)≤2ϱm,whereϱ=(b1/2-a1/2)(b1/2+a1/2).This implies that if the number of iterations satisfiesm≥lnϱ-1ln(2/ɛ), that is, in particular if(29)m≥12(ba)1/2ln⁡⁡(2/ɛ),ɛ>0,
then ∥em∥/∥e0∥≤ɛ.The disadvantage with this method is that to make it effective one needs accurate estimates ofa and b, and we need to determine m beforehand (which, however, can be done by (29)). The method cannot utilize any special distribution of the eigenvalues in the spectrum (as opposed to the conjugate gradient method, see below). More important, however, is that this method is actually numerically unstable (similarly to an explicit time stepping method for initial value problems when the time steps are too large). This is due to the fact that ∥I-τkC-1A∥ is much larger than unity for several of the values τk. However, one may prove that with some particular permutation of the parameters, their instability effect can be avoided.There is an alternative to the choice (25). Namely, there exists a three term form of the Chebyshev acceleration method(30)xk+1=αkxk+(1-αk)xk-1-βkC-1rk,k=1,2,…,wherex1=x0-(1/2)β0C-1r0.Here, the parameters are chosen asβ0=4/(a+b),(31)αk=a+b2βk,βk-1=a+b2-(b-a4)2βk-1k=1,2,….Hence, we do not have to determine the number of steps beforehand. More importantly, it has been shown in [18] that this method is numerically stable. (For some related remarks, see [6]). A similar form of the method was proposed a long time ago, see Golub and Varga [19] and the references cited therein.It is interesting to note that the parameters approach stationary values. IfC-1A=I-B and B has eigenvalues in [-ϱ,ϱ], ϱ=ϱ(B)<1 (the spectral radius of B), then(32)a=1-ϱ,b=1+ϱ,αk=a+b2βk⟶2[1+(1-ϱ2)1/2],which is recognized as the parameter ωopt of the optimal SOR method (see Section 2). Young [20] has proven that the asymptotic rate of convergence is retained even if one uses the stationary values throughout the iterations.For the case of complex eigenvalues ofC-1A with positive real parts and contained in an ellipse one may choose parameters similarly. See [6, 21, 22] for details. For comments on the optimaly of the method, see [23]. For application of the method for nonsymmetric problems, see [6, 24].Perhaps the main thrust during the 1970 has been in using the conjugate gradient method as an acceleration method. Already much has been written on the subject; we refer to [25–27] for a historical account, to [18, 28–33] for an exposition of the preconditioned conjugate gradient PCG method and to [18, 34] for a survey of generalized and truncated gradient methods for nonsymmetric and indefinite matrix problems.The advantage with conjugate gradient methods is that they are self adaptive; the optimal parameters are calculated by the algorithm so that the error in energy norm∥el∥A1/2={(el)TAel}1/2 is minimized. This applies to a problem where C and A are symmetric and positive definite (SPD) or, more generally, if C-1A is similarly equivalent to an SPD matrix. Hence, there is no need to know any bounds for the spectrum. Since the method converges at least as fast as the Chebyshev method it follows that ∥x-xm∥A1/2≤ɛ∥x-x0∥A1/2, if(33)m=int⁡{12K1/2ln⁡(2ɛ)+1}.We describe now the conjugate gradient method. Thereby we follow the presentations in [29, 35].


### 3.2. The Preconditioned Conjugate Gradient Method

During the past 40 years or so, the preconditioned conjugate gradient method has become the major iterative solution method for linear systems of algebraic equations, in particular those arising in science and engineering. The author of these notes became interested in the method by the beginning of 1970 (cf. [18]).The conjugate gradient algorithm to solve a system of linear equations theAx=b, where A(n×n) is symmetric and positive definite, was originally introduced by Hestenes and Stiefel [25] in 1950. Before we discuss the properties of the CG method, we describe its implementation. Namely, the algorithm consists of the steps in Algorithm 1.Algorithm 1:Standard conjugate gradient algorithm.
Givenx(0),ɛ              initial guess and absolute orrelative stopping toleranceSetx(0),g=Ax-b,δ0=gTgd=-g      initial search directionRepeat       until convergenceh=Adτ=δ0/(dTh)x=x+τd     new approximationg=g+τh     new (iterative) residualδ1=gTgifδ1≤ɛ then stop,  otherwiseβ=δ1/δ0,δ0=δ1d=-g+βd      new search directionWhat one sees from a first glance is that the CG algorithm is quite simple. Each iteration consists of one matrix-vector multiplication, two vector updates and two scalar products. Apart from the initial guessx(0) (which can be taken to be the zero vector) and stopping tolerance, there are no other method parameters to be determined or tuned by the user. Thus, the method is easily programmable, cheap in terms of arithmetic operations and performs as a black box.For some problems the standard (unpreconditioned) CG method performs impressively well and this can be explained by some particular properties of this powerful algorithm.The CG method is best described as a method to minimize the quadratic functional(34)f(x)=12xTAx-bTx+c,
over a set of vectors. If A is nonsingular, then we can rewrite f in the form(35)f(x)=12xT(Ax-b)TA-1(Ax-b)-12bTA-1b+c,
so, minimizing the quadratic functional is equivalent to solving the system Ax=b. If A is singular and A-1 in (35) is replaced by a generalized inverse of A, then the above equivalence still holds if the minimization takes place on a subspace in the orthogonal complement to the null-space of A.Given an initial approximationx(0) and the corresponding residual r(0)=Ax(0)-b, the minimization in the conjugate gradient method takes place successively on a subspace(36)Kk={r(0),Ar(0),A2r(0),…,Ak-1r(0)},of growing dimension. This subspace is referred to as the Krylov set.In the derivation of the algorithm, the next approximate solution is constructed as(37)x(k+1)=x(k)+τkd(k),where τk is chosen(38)τk=-d(k)Tg(k)d(k)TAd(k)=-d(k)T(Ax(k)-b)d(k)TAd(k),which minimizes the functionf(x(k)+τd(k)),  -∞<τ<∞. Also, the gradient of f at x(k+1) is made orthogonal to the search direction d(k). This is seen from the following relations:(39)x(k+1)=x(k)+τkd(k)⟹Ax(k+1)-b=Ax(k)-b+τkAd(k)⟹g(k+1)=g(k)+τkAd(k)⟹d(k)Tg(k+1)=d(k)Tg(k)+τkd(k)TAd(k)=0.As in Fourier type minimization methods, it turns out to be efficient to work with orthogonal (A-orthogonal) search directions d(k) which, since A is symmetric, can be determined from a three-term recursion(40)d(0)=r(0),d(k+1)=-Ad(k)+β̃kd(k),k=1,2,…,or equivalently, from(41)d(k+1)=-r(k+1)+βkd(k).This recursive choice of search directions is done so that at each step the solution has smallest error in theA-norm, ∥x-x(k)∥A={e(k)TAe(k)}1/2, where e(k)=x-x(k) is the iteration error. As mentioned, the minimization takes place over the set of (Krylov) vectors 𝒦k and, as is readily seen(42)Kk={x(1)-x(0),x(2)-x(0),…,x(k)-x(0)}={g(0),g(1),…,g(k-1)}={d(0),d(1),…,d(k-1)}.To summarize, the CG method possesses the following remarkable properties.Theorem 5.
Let the CG Algorithm1 be applied to a symmetric positive definite matrix A. Then, in exact arithmetic the following properties hold: (1)
the iteratively constructed residualsg are mutually orthogonal, that is, g(k)Tg(j)=0, j<k;(2)
the search directionsd are A-orthogonal (or conjugate), that is, d(k)TAd(j)=0, j<k:(3)
as long as the method has not converged, that is,g(k)≠0, the algorithm proceeds with no breakdown and (42) holds;(4)
as long as the method has not converged, the newly constructed approximationx(k) is the unique point in x(0)⊕𝒦k that minimizes ∥e(k)∥A=∥x-x(k)∥A,(5)
the convergence is monotone inA-norm, that is, ∥e(k)∥A<∥e(k-1)∥A and e(m)=0 will be achieved for some m≤n.For a proof of the above theorem consult, for instance, [29] or [6].Since the method is optimal, that is it gives the smallest error on a subspace of growing dimension, itterminates with the exact solution (ignoring round-off errors) in at most n steps (the dimension of the whole vector space x∈Rn). In fact, it can be readily seen that the CG algorithm terminates after m steps, where m is the degree of the minimal polynomial Qm to A with respect to the initial residual vector, in other words, Qm has the smallest degree of all polynomials Q for which Q(A)r(0)=0. Therefore, the CG method can be viewed also as a direct solution method. However, in practice we want convergence to occur to an acceptable accuracy in much fewer steps than n or m. Thus, we use CG as an iterative method.For further discussions of the CG methods, see [6, 34, 35].When one experiments with CG to solve systems with various matrices one observes some phenomena which need special attention. This can be illustrated by a simple example.Consider the solution ofAx=b by the standard conjugate gradient, where(43)A=[2-1-12-1⋱⋱-12-1-11],b=[10000].The exact solution isx̂=[1,1,…,1]T. Starting with x(0)=[0,0,…,0]T one finds that after k iterations(44)x(k)=[kk+1,k-1k+1,…,1k+1,0,…,0]T,for 1≤k≤n-1 and x(n)=x̂. Hence, the information travels one step at a time from left to right and it takes n steps before the last component has changed at all. The algorithm converges exactly in n steps and terminates due to the final recurrence property of the method.Another detail one observes is that the norm of the error,∥x-x(k)∥, can be much larger than the norm of the iteratively computed residual.These examples illustrate the fact that although the method has an optimal order of convergence rate in theenergy norm, its actual convergence rate in spectral norm can be different and depends both on the distribution of the eigenvalues of the (preconditioned) system matrix and on the initial approximation (or residual). (For comparison, we note that the rate of convergence for steepest descent depends only on the ratio of the extremal eigenvalues of A.) Faster convergence for the CG method is expected when the eigenvalues are clustered.One way to get a better eigenvalue distribution is to preconditionA by a proper preconditioner B. Hence, in order to achieve a better eigenvalue distribution it is crucial in practice to use some form of preconditioning, that is, a matrix B which approximates A in some sense, which is relatively cheap to solve systems with and for which the spectrum of B-1A (equivalently B-1/2AB-1/2 if B is s.p.d.) is more favorable for the convergence of the CG method. As it turns out, if B is symmetric and positive definite, the corresponding preconditioned version, the PCG method, is best derived by replacing the inner product with (u,v)=uTBv. It takes the following form, see Algorithm 2.Algorithm 2:Preconditioned conjugate gradient algorithm.
Givenx(0),ɛ         initial guess andstopping toleranceSetx(0),g=Ax-b,h=[B]-1gδ0=gThd=-h        initial search directionRepeat    until convergenceh=Adτ=δ0/(dTh)x=x+τd       new approximationg=g+τh       new (iterative) residualδ1=gTgh=[B]-1g       new pseudoresidualδ1=gThifδ1≤ɛ then stopβ=δ1/δ0,δ0=δ1d=-h+βd        new search directionHere,[B]-1 denotes the action of B-1, that is, one does not multiply with the inverse matrix B-1, but normally solves a linear system with matrix B.In order to understand what is wanted of agood preconditioning matrix, we discuss first some issues of major importance related to the rate of convergence of the CG method. Thereby it becomes clear that the standard spectral condition number is often too simple to explain the detailed convergence behaviour. In particular we discuss the sub- and superlinear convergence phases frequently observed in the convergence history of the conjugate gradient method.A preconditioner can be applied in two different manners, namely, asB-1A or BA. The first form implies the necessity to solve a system with B at each iteration step while the second form implies a matrix-vector multiplication with B (a multiplicative preconditioner). In the latter, case B can be seen as an approximate inverse of A. One can also use a hybrid form αB1-1+βB2.The presentation here is limited to symmetric positive semidefinite matrices. It is based mainly on the articles [29, 31].


### 3.3. On the Rate of Convergence Estimates of the Conjugate Gradient Method

LetA be symmetric, positive semidefinite and consider the solution of Ax=b by a preconditioned conjugate gradient method. In order to understand how an efficient preconditioner to A should be chosen we must first understand some general properties of the rate of convergence of conjugate gradient methods.


#### 3.3.1. Rate of Convergence Estimates Based on  Minimax Approximation

As is well known (see e.g., [6, 30]), the conjugate gradient method is a norm minimizing method. For the preconditioned standard CG method, we have(45)‖ek‖A=min⁡Pk∈πk‖Pk(B)e0‖A,
where ∥u∥A={uTAu}1/2, ek=x-xk is the iteration error and πk denotes the set of polynomials of degree k which are normalized at the origin, that is, Pk(0)=1. This is a norm on the subspace orthogonal to the mullspace of A, that is, on the whole space, if A is nonsingular.Consider theC-innerproduct (u,v)=uTCv, and note that B=C-1A is symmetric with respect to this innerproduct, let v1,v2,…,vn be orthonormal eigenvectors and let λi,  i=1,…,n be the corresponding eigenvalues of B. Let(46)e0=∑j=1nαjvj,be the eigenvector expansion of the initial vector where αj=(e0,vi),  i=1,…,n. Note further that the eigenvectors are both A- and C-orthogonal. Then, by the construction of the CG method, it follows(47)ek=∑j=1nαjPk(λj)vj,and, using the nonnegativity of the eigenvalues, we find(48)‖ek‖A=‖∑j=1nαjPk(λj)vj‖A={∑j=1nαj2λjPk2(λj)}1/2≤{∑1,λj>0nαj2λj}1/2max⁡1≤i≤nλi>0|Pk(λi)|=max⁡1≤i≤nλi>0|Pk(λi)|‖e0‖A.
Here we have used the A-orthogonality of the eigenvectors. Similarly, using the C-orthogonality, we find(49)‖ek‖C={∑j=1nαj2Pk2(λj)}1/2≤max⁡1≤i≤n|Pk(λj)|‖e0‖C.Due to the minimization property (45) there follows from (48) the familiar bound(50)‖ek‖A≤min⁡Pk∈πk1max⁡1≤i≤nλi>0|Pk(λi)|‖e0‖A.
Estimate (50) is sharp in the respect that for every k there exists an initial vector for which equality is attained. In fact, for such a vector we necessarily have that αj≠0 if and only if αj belongs to a set of k+1 points ( the so-called Haar condition) where maxi|Pk(λi)| is taken. For such an initial vector (49) shows that, if the eigenvalues are positive, we have also(51)‖ek‖C=min⁡Pk∈πk1max⁡1≤i≤n|Pk(λi)|‖e0‖C.The rate of convergence of the iteration error∥ek∥A is measured by the average convergence factor(52){‖ek‖A‖e0‖A}1/k.Inequality (50) shows that this can be majorized with an estimate of the rate of convergence of a best polynomial approximation problem (namely the best approximation of the function ≡0, of polynomials in πk1) in maximum norm on the discrete set formed by the spectrum of B. Clearly, multiple eigenvalues are treated as single so the actual approximation problem is(53)min⁡Pk∈πk1max⁡1≤i≤m|Pk(λ̃i)|,
where the disjoint positive eigenvalues λ̃j have been ordered in increasing value, 0<λ̃i<⋯<λ̃m, and m is the number of such eigenvalues. However, the solution of this problem requires knowledge of the spectrum, which is not available in general. Even if it is known, the estimate (53) can be troublesome in practice, since it involves approximation on a general discrete set of points.Besides being costly to apply, such estimates do not give any qualitative insight in the behaviour of the conjugate gradient method for various typical eigenvalue distributions.That is why we make some further assumptions on the spectrum in order to simplify the approximation problem and at the same time, present estimates which can be used both to estimate the number of iterations and to give some insight in the qualitative behaviour of the iteration method.


#### 3.3.2. Standard Condition Number: Linear Convergence

If the eigenvalues are (densely) located in an interval[a,b] where a>0, we can majorize the best approximation problem on the discrete set with the best approximation problem on the interval and frequently still get a good estimate. We have(54)min⁡Pk∈πk1max⁡1≤i≤m|Pk(λ̃i)|≤min⁡Pk∈πk1max⁡a≤x≤b|Pk(x)|.The solution to thisminmax problem is well known and uses Chebyshev polynomials. One finds that(55)min⁡Pk∈πk1max⁡a≤x≤b|Pk(x)|=1Tk((b+a)/(b-a))=2σk1+σ2k,
where σ=(1-a/b)/(1+a/b), and Tk(x)=(1/2)[(x+x2-1)+(x-x2-1)k],  a=λ̃1,  b=λ̃m. Hence, the average rate of convergence of the upper bound approaches σ as k→∞. Also, it is readily found (see [35]) that the relative iteration error ∥ek∥A/∥e0∥A≤ɛ if(56)k=k*(a,b,ɛ)=⌈ln⁡⁡((1/ɛ)+(1/ɛ2)-1)ln⁡⁡σ-1⌉.
Here ⌈ξ⌉ denotes the smallest integer not less than ξ.It turns out that the above holds more generally ifA is nonsymmetric but the eigenvalues are contained in an ellipse with foci a, b, where b≥a>0, if one replaces σ with σ̂=σ(1+δ)/(1-δ), where δ is the eccentricity of the ellipse, (i.e., the ratio of the semiaxes) and δ<2a/b/(1+a/b). Also, in a similar way, the case of eigenvalues contained in two separate intervals or ellipses can be analysed, see, for example, [35] for further details.Whenb/a→∞,  δ=0, and ɛ→0 the following upper bound becomes an increasingly accurate replacement of (56),(57)k*≤⌈12baln⁡⁡2ɛ⌉.
The above estimate of the rate of convergence and of the number of iterations show that they depend only on the condition number b/a and on the eccentricity of the ellipse, containing the eigenvalues. Therefore, except in special cases, this estimate is not very accurate. When we use a more detailed information of the spectrum and the initial error vector, sometimes substantially better estimates can be derived. This holds for instance when there are well separated small and/or large eigenvalues. Before we consider this important case, we mention briefly another similar minimax result which holds when we use different norms for the iteration error vector and for the initial vector.By (48), we have(58)‖ek‖A={∑αj2λjPk2(λj)}1/2={∑αj2λj1-2sλj2sPk2(λj)}1/2≤min⁡Pk∈πk1max⁡1≤λj≤m|λjsPk(λj)|{∑αj2λj(1-2s)}1/2=min⁡Pk∈πk1max⁡1≤λj≤m|λjsPk(λj)|‖e0‖A1-2s.If the initial vector is such that Fourier coefficients for the highest eigenvalue modes are dominating, then∥e0∥A1-2s may exist and take not too large values even for some s≥1/2. We consider the interesting case where s≥1/2, for which the following theorem holds (see [6, 36]).Theorem 6.
Letπk1 denote the set of polynomials of degree k such that Pk(0)=1. Then for k=1,2… and for any s≥1/2 such that 2s is an integer, it holds
(59)‖ek‖A‖e0‖A1-2s≤min⁡Pk∈πk1max⁡0≤x≤1|xsPk(x)|≤(sk+s)2s.Remark 1.
Fors=1/2 it holds
(60)max⁡0≤x≤1|x1/2Pk(x)|=12k+1,
for Pk(x)=U2k(1-x) and for s=1, it holds
(61)max⁡0≤x≤1|xPk(x)|=1k+1tanπ4k+4<1(k+1)2,
for Pk(x)=(x-1(-1)k)/(k+1)tan(π/(4k+4))Tk+1((1+cos(π/(2k+2)))x-cos(π/(2k+2))) where Tk(x) and Uk(x) are the Chebyshev polynomials of kth degree of the first and second kind, respectively.For other values (59) is an upper bound only, that is, not sharp. At any rate, it shows that the error ∥ek∥A converges (initially) at least as fast as (s/(k+s))2s, that is, as 1/(2k+1) for s=1/2 and as (1/(k+1))2 for s=1.Note that this convergence rate does not depend on the eigenvalues, in particular not on the spectral condition number.Conclusion.
By computing the initial approximation vector from a coarse mesh, the components fore0 for the first Fourier modes will be small and ∥e0∥A1-2s may take on values that are not very large even when s=1 or bigger. Therefore, there is an initial decay of the residual as O(k2s), independent of the condition number. Note, however, that the actual errors may not have decayed sufficiently even if the residual has.We consider now upper bound estimates which show how the convergence history may enter a superlinear rate.An Estimate to Show a Superlinear Convergence Rate Based on the K-Condition Number
A somewhat rough, but simple and illustrative superlinear convergence estimate can be obtained in terms of the so calledK-condition number, (see [37, 38])
(62)K=K(B)=((1/n)tr⁡(B))ndet⁡(B)=(1n∑i=1nλi)n(Πi=1nλi)-1,
where we assume that B is s.p.d.
Note thatK1/n equals the quotient between the arithmetic and geometric averages of the eigenvalues. This quantity is similar to the spectral condition number κ(B) in that it is never smaller than 1, and is equal to 1 if and only if B=αI,  α>0 (recall that B is symmetrizable).
Based on theK-condition number, a superlinear convergence result can be obtained as follows.Theorem 7.
Letk<n be even and k≥3lnK. Then,
(63)∥ek∥A∥e0∥A≤(3ln⁡Kk)k/2.Proof.
Letk=2m and the polynomial Pk be of a simplest possible form, that is, let it vanish at the m smallest and m largest eigenvalues of B. As follows from (48), we have then
(64)‖ek‖A‖e0‖A≤max⁡λ1≤λ≤λn|Πi=1m(1-λλi)(1-λλn+1-i)|=Πi=1mmax⁡λi≤λ≤λn+1-i(λλi-1)(1-λλn+1-i)=Πi=1m((λi+λn+1-i)24λiλn+1-i-1)≤((Πi=1m(λi+λn+1-i)24λiλn+1-i)1/m-1)m.
The latter follows from (Πi=1m(1-Θi))1/m+(Πi=1mΘi)1/m≤1 with Θi=4λiλn+1-i(λi+λn+1-i)2. Using now twice the inequality between the arithmetic and geometric mean values, one has
(65)Πi=1m(λi+λn+1-i)24λiλn+1-i≤((1/n)-2m∑i=m+1n-mλi)n-2mΠi=m+1n-mλi×Πi=1m(λi+λn+1-i/2)2Πi=1mλiλn+1-i≤(1/n(∑i=m+1n-mλi+∑i=1m(λi+λn+1-i)))nΠi=m+1n-mλiΠi=1mλiλn+1-i=K.
Using exp (2x)-1≤2x/(1-x),  x<1, we get the required estimate,
(66)‖ek‖A∥e0∥A≤(K2/k-1)k/2≤(2ln⁡Kk-ln⁡K)k/2≤(2ln⁡K2k/3(k/3-ln⁡K))k/2≤(3ln⁡Kk)k/2.A somewhat better result of the same type exists. The estimate is of similar type, that is,(67)∥rk∥C-1∥r0∥C-1≤(K1/k-1)k/2,
where rk=Aek was obtained using more complicated techniques, see [6, 37], and the references quoted therein. Note that here as k/lnK→∞, we have(68)∥ek∥C-1∥e0∥C-1≤(e1/kln⁡K-1)k/2≈(1+ln⁡K2k)k/2(ln⁡Kk)k/2≈K1/4(ln⁡Kk)k/2,that is, the simpler upper bound in (63) is asymptotically worse than this bound (albeit in a different norm) by the factor 3k/2/lnK1/4.The upper bounds in the above estimates involve a convergence factor which decreases with increasing iteration number and show hence a superlinear rate of convergence. Note, however, thatK1/n≤κ(B)≲4K (see [6]) where κ(B)=λn/λ1 is the spectral condition number, so the K-condition number may take very large values when κ(B) is large.The estimates based on theK-condition number involve only “integral” characteristics of the preconditioned matrix (the trace and the determinant). Sometimes, it is possible to obtain a practical estimate of K(B) which can be useful for the a priori construction of good preconditioners and for the a posteriori assessment of their quality, see Section 5 for further details.The estimate(69)‖rk‖C-1‖r0‖C-1≤(K1/k-1)k/2≤ɛshows that(70)K1/2(1-K-1/k)k/2≤ɛor(71)k2log⁡211-K-1/k≥log⁡2K1/2ɛ,which holds if(72)k>log⁡2K+2log⁡21ɛ.
Hence, when K≫ɛ-2 the estimated number of iterations depends essentially only on log2K, that is, depends little on the relative accuracy ɛ, which indicates a fast superlinear convergence, when k>log2K.When actually estimating the number of iterations, Theorem6 shows a useful result only when k>O(lnK)=n(lnK1/n), that is, the quotient between the arithmetic and geometric averages of the eigenvalues, which equals K1/n, must be close to unity and the eigenvalues must be very well clustured so K1/n=1+O(n-ɛ) for some ɛ>0; otherwise the estimated number of iterations will be O(n), which is normally a useless result. The next example illustrates this further.Example 1.
Consider a geometric distribution of eigenvalues ofA,   λj=js,   j=1,2,…,n for some positive s. Here, asymptotically
(73)tr⁡(A)=∑1njs~1s+1ns+1,n⟶∞.
Using Stirling's formula, we find
(74)det⁡(A)=∏1nλj=(∏1nj)s~(2πn)s/2(ne)ns,n⟶∞,
so,
(75)κ(A)1/n~ess+1,n⟶∞.
Hence, s must be sufficiently small for the estimate in Theorem 6. to be useful. On the other hand, the spectral condition number (i) κ(A)=ns, and the simple estimate based on κ(A) leads to k~O(ns/2) and gives hence, asymptotically, a smaller upper bound when s<2. For further discussions on superlinear rate of convergence, see [39].


### 3.4. Generalized Conjugate Gradient Methods

The rate of convergence estimates as given above, holds for a restricted class of matrices, symmetric or, more generally, for normal matrices.To handle more general classes of problems for which such optimal rate of convergence results as in (45) holds, one needs more involved methods. Much work has been devoted to this problem. This includes methods like generalized minimum residual (GMRES), see Saad and Schultz [40], generalized conjugate residual (GCR), and generalized conjugate gradient (GCG), see [6] and for further details [41]. As opposed to the standard conjugate and gradient method, they require a long version of updates for the search directions, as the newest search direction at each stage is not in general, automatically (in exact precision) orthogonal to the previous search directions, but must be orthogonalized at each step. This makes the computational expense per step grow linearly and the total expense grows quadratically with the iteration index. In addition, due to finite precision, there is a tendency of loss of orthogonality, even for symmetric problems when many iterations are required. One remedy which has been suggested is to use the method only for a few steps, say 10, and restart the method with the current approximation as initial approximation.Clearly, however, in this way, the optimal convergence property of the whole Krylov set of vector is lost. For this and other possibilities, see, for example, [42].Another important version of the generalized conjugate gradient methods occurs when one uses variable preconditioners. Variable preconditioners, that is, a preconditioner changed from one iteration to the next iteration step, are used in many contexts.For instance, one can use variable drop tolerance, computed adaptively, in an incomplete factorization method (see Section4). When the given matrix is partitioned in two by two blocks, it can be efficient to use inner iterations when solving arising systems for one, or both, of the diagonal block matrices, see, for example, [43], and the flexible conjugate gradient method in Saad, [44, 45].Due to space limitations, the above topics will not be further discussed in this paper.


## 3.1. The Chebyshev Iterative Method

In case the eigenvalues ofC-1A are real and positive and if a positive lower (a) and (b) an upper bound are known of the spectrum, then we see that {τk} should be chosen such that maxa≤λ≤b|Pm(λ)| is minimized over all Pm∈Πm0, that is, over the set of polynomials of degree m satisfying Pm(0)=1.The solution to this min-max problem is well known,(26)Pm(λ)=Tm((b+a-2λ)/(b-a))Tm((b+a)/(b-a)),where Tm(z)=(1/2)[(z+(z2-1)1/2)m+(z-(z2-1)1/2)m]=cos(marccosz) are the Chebyshev polynomials of the first kind. The corresponding values of τk satisfy(27)1τk=b-a2cos⁡Θk+b+a2,Θk=2k-12mπ,k=1,2,…,m,which are the zeros of the polynomial. The corresponding method is referred to as the Chebyshev (one-step) acceleration method, see, for example, [10, 18]. It is an easy matter to show that(28)1Tm(b+ab-a)≤2ϱm,whereϱ=(b1/2-a1/2)(b1/2+a1/2).This implies that if the number of iterations satisfiesm≥lnϱ-1ln(2/ɛ), that is, in particular if(29)m≥12(ba)1/2ln⁡⁡(2/ɛ),ɛ>0,
then ∥em∥/∥e0∥≤ɛ.The disadvantage with this method is that to make it effective one needs accurate estimates ofa and b, and we need to determine m beforehand (which, however, can be done by (29)). The method cannot utilize any special distribution of the eigenvalues in the spectrum (as opposed to the conjugate gradient method, see below). More important, however, is that this method is actually numerically unstable (similarly to an explicit time stepping method for initial value problems when the time steps are too large). This is due to the fact that ∥I-τkC-1A∥ is much larger than unity for several of the values τk. However, one may prove that with some particular permutation of the parameters, their instability effect can be avoided.There is an alternative to the choice (25). Namely, there exists a three term form of the Chebyshev acceleration method(30)xk+1=αkxk+(1-αk)xk-1-βkC-1rk,k=1,2,…,wherex1=x0-(1/2)β0C-1r0.Here, the parameters are chosen asβ0=4/(a+b),(31)αk=a+b2βk,βk-1=a+b2-(b-a4)2βk-1k=1,2,….Hence, we do not have to determine the number of steps beforehand. More importantly, it has been shown in [18] that this method is numerically stable. (For some related remarks, see [6]). A similar form of the method was proposed a long time ago, see Golub and Varga [19] and the references cited therein.It is interesting to note that the parameters approach stationary values. IfC-1A=I-B and B has eigenvalues in [-ϱ,ϱ], ϱ=ϱ(B)<1 (the spectral radius of B), then(32)a=1-ϱ,b=1+ϱ,αk=a+b2βk⟶2[1+(1-ϱ2)1/2],which is recognized as the parameter ωopt of the optimal SOR method (see Section 2). Young [20] has proven that the asymptotic rate of convergence is retained even if one uses the stationary values throughout the iterations.For the case of complex eigenvalues ofC-1A with positive real parts and contained in an ellipse one may choose parameters similarly. See [6, 21, 22] for details. For comments on the optimaly of the method, see [23]. For application of the method for nonsymmetric problems, see [6, 24].Perhaps the main thrust during the 1970 has been in using the conjugate gradient method as an acceleration method. Already much has been written on the subject; we refer to [25–27] for a historical account, to [18, 28–33] for an exposition of the preconditioned conjugate gradient PCG method and to [18, 34] for a survey of generalized and truncated gradient methods for nonsymmetric and indefinite matrix problems.The advantage with conjugate gradient methods is that they are self adaptive; the optimal parameters are calculated by the algorithm so that the error in energy norm∥el∥A1/2={(el)TAel}1/2 is minimized. This applies to a problem where C and A are symmetric and positive definite (SPD) or, more generally, if C-1A is similarly equivalent to an SPD matrix. Hence, there is no need to know any bounds for the spectrum. Since the method converges at least as fast as the Chebyshev method it follows that ∥x-xm∥A1/2≤ɛ∥x-x0∥A1/2, if(33)m=int⁡{12K1/2ln⁡(2ɛ)+1}.We describe now the conjugate gradient method. Thereby we follow the presentations in [29, 35].


## 3.2. The Preconditioned Conjugate Gradient Method

During the past 40 years or so, the preconditioned conjugate gradient method has become the major iterative solution method for linear systems of algebraic equations, in particular those arising in science and engineering. The author of these notes became interested in the method by the beginning of 1970 (cf. [18]).The conjugate gradient algorithm to solve a system of linear equations theAx=b, where A(n×n) is symmetric and positive definite, was originally introduced by Hestenes and Stiefel [25] in 1950. Before we discuss the properties of the CG method, we describe its implementation. Namely, the algorithm consists of the steps in Algorithm 1.Algorithm 1:Standard conjugate gradient algorithm.
Givenx(0),ɛ              initial guess and absolute orrelative stopping toleranceSetx(0),g=Ax-b,δ0=gTgd=-g      initial search directionRepeat       until convergenceh=Adτ=δ0/(dTh)x=x+τd     new approximationg=g+τh     new (iterative) residualδ1=gTgifδ1≤ɛ then stop,  otherwiseβ=δ1/δ0,δ0=δ1d=-g+βd      new search directionWhat one sees from a first glance is that the CG algorithm is quite simple. Each iteration consists of one matrix-vector multiplication, two vector updates and two scalar products. Apart from the initial guessx(0) (which can be taken to be the zero vector) and stopping tolerance, there are no other method parameters to be determined or tuned by the user. Thus, the method is easily programmable, cheap in terms of arithmetic operations and performs as a black box.For some problems the standard (unpreconditioned) CG method performs impressively well and this can be explained by some particular properties of this powerful algorithm.The CG method is best described as a method to minimize the quadratic functional(34)f(x)=12xTAx-bTx+c,
over a set of vectors. If A is nonsingular, then we can rewrite f in the form(35)f(x)=12xT(Ax-b)TA-1(Ax-b)-12bTA-1b+c,
so, minimizing the quadratic functional is equivalent to solving the system Ax=b. If A is singular and A-1 in (35) is replaced by a generalized inverse of A, then the above equivalence still holds if the minimization takes place on a subspace in the orthogonal complement to the null-space of A.Given an initial approximationx(0) and the corresponding residual r(0)=Ax(0)-b, the minimization in the conjugate gradient method takes place successively on a subspace(36)Kk={r(0),Ar(0),A2r(0),…,Ak-1r(0)},of growing dimension. This subspace is referred to as the Krylov set.In the derivation of the algorithm, the next approximate solution is constructed as(37)x(k+1)=x(k)+τkd(k),where τk is chosen(38)τk=-d(k)Tg(k)d(k)TAd(k)=-d(k)T(Ax(k)-b)d(k)TAd(k),which minimizes the functionf(x(k)+τd(k)),  -∞<τ<∞. Also, the gradient of f at x(k+1) is made orthogonal to the search direction d(k). This is seen from the following relations:(39)x(k+1)=x(k)+τkd(k)⟹Ax(k+1)-b=Ax(k)-b+τkAd(k)⟹g(k+1)=g(k)+τkAd(k)⟹d(k)Tg(k+1)=d(k)Tg(k)+τkd(k)TAd(k)=0.As in Fourier type minimization methods, it turns out to be efficient to work with orthogonal (A-orthogonal) search directions d(k) which, since A is symmetric, can be determined from a three-term recursion(40)d(0)=r(0),d(k+1)=-Ad(k)+β̃kd(k),k=1,2,…,or equivalently, from(41)d(k+1)=-r(k+1)+βkd(k).This recursive choice of search directions is done so that at each step the solution has smallest error in theA-norm, ∥x-x(k)∥A={e(k)TAe(k)}1/2, where e(k)=x-x(k) is the iteration error. As mentioned, the minimization takes place over the set of (Krylov) vectors 𝒦k and, as is readily seen(42)Kk={x(1)-x(0),x(2)-x(0),…,x(k)-x(0)}={g(0),g(1),…,g(k-1)}={d(0),d(1),…,d(k-1)}.To summarize, the CG method possesses the following remarkable properties.Theorem 5.
Let the CG Algorithm1 be applied to a symmetric positive definite matrix A. Then, in exact arithmetic the following properties hold: (1)
the iteratively constructed residualsg are mutually orthogonal, that is, g(k)Tg(j)=0, j<k;(2)
the search directionsd are A-orthogonal (or conjugate), that is, d(k)TAd(j)=0, j<k:(3)
as long as the method has not converged, that is,g(k)≠0, the algorithm proceeds with no breakdown and (42) holds;(4)
as long as the method has not converged, the newly constructed approximationx(k) is the unique point in x(0)⊕𝒦k that minimizes ∥e(k)∥A=∥x-x(k)∥A,(5)
the convergence is monotone inA-norm, that is, ∥e(k)∥A<∥e(k-1)∥A and e(m)=0 will be achieved for some m≤n.For a proof of the above theorem consult, for instance, [29] or [6].Since the method is optimal, that is it gives the smallest error on a subspace of growing dimension, itterminates with the exact solution (ignoring round-off errors) in at most n steps (the dimension of the whole vector space x∈Rn). In fact, it can be readily seen that the CG algorithm terminates after m steps, where m is the degree of the minimal polynomial Qm to A with respect to the initial residual vector, in other words, Qm has the smallest degree of all polynomials Q for which Q(A)r(0)=0. Therefore, the CG method can be viewed also as a direct solution method. However, in practice we want convergence to occur to an acceptable accuracy in much fewer steps than n or m. Thus, we use CG as an iterative method.For further discussions of the CG methods, see [6, 34, 35].When one experiments with CG to solve systems with various matrices one observes some phenomena which need special attention. This can be illustrated by a simple example.Consider the solution ofAx=b by the standard conjugate gradient, where(43)A=[2-1-12-1⋱⋱-12-1-11],b=[10000].The exact solution isx̂=[1,1,…,1]T. Starting with x(0)=[0,0,…,0]T one finds that after k iterations(44)x(k)=[kk+1,k-1k+1,…,1k+1,0,…,0]T,for 1≤k≤n-1 and x(n)=x̂. Hence, the information travels one step at a time from left to right and it takes n steps before the last component has changed at all. The algorithm converges exactly in n steps and terminates due to the final recurrence property of the method.Another detail one observes is that the norm of the error,∥x-x(k)∥, can be much larger than the norm of the iteratively computed residual.These examples illustrate the fact that although the method has an optimal order of convergence rate in theenergy norm, its actual convergence rate in spectral norm can be different and depends both on the distribution of the eigenvalues of the (preconditioned) system matrix and on the initial approximation (or residual). (For comparison, we note that the rate of convergence for steepest descent depends only on the ratio of the extremal eigenvalues of A.) Faster convergence for the CG method is expected when the eigenvalues are clustered.One way to get a better eigenvalue distribution is to preconditionA by a proper preconditioner B. Hence, in order to achieve a better eigenvalue distribution it is crucial in practice to use some form of preconditioning, that is, a matrix B which approximates A in some sense, which is relatively cheap to solve systems with and for which the spectrum of B-1A (equivalently B-1/2AB-1/2 if B is s.p.d.) is more favorable for the convergence of the CG method. As it turns out, if B is symmetric and positive definite, the corresponding preconditioned version, the PCG method, is best derived by replacing the inner product with (u,v)=uTBv. It takes the following form, see Algorithm 2.Algorithm 2:Preconditioned conjugate gradient algorithm.
Givenx(0),ɛ         initial guess andstopping toleranceSetx(0),g=Ax-b,h=[B]-1gδ0=gThd=-h        initial search directionRepeat    until convergenceh=Adτ=δ0/(dTh)x=x+τd       new approximationg=g+τh       new (iterative) residualδ1=gTgh=[B]-1g       new pseudoresidualδ1=gThifδ1≤ɛ then stopβ=δ1/δ0,δ0=δ1d=-h+βd        new search directionHere,[B]-1 denotes the action of B-1, that is, one does not multiply with the inverse matrix B-1, but normally solves a linear system with matrix B.In order to understand what is wanted of agood preconditioning matrix, we discuss first some issues of major importance related to the rate of convergence of the CG method. Thereby it becomes clear that the standard spectral condition number is often too simple to explain the detailed convergence behaviour. In particular we discuss the sub- and superlinear convergence phases frequently observed in the convergence history of the conjugate gradient method.A preconditioner can be applied in two different manners, namely, asB-1A or BA. The first form implies the necessity to solve a system with B at each iteration step while the second form implies a matrix-vector multiplication with B (a multiplicative preconditioner). In the latter, case B can be seen as an approximate inverse of A. One can also use a hybrid form αB1-1+βB2.The presentation here is limited to symmetric positive semidefinite matrices. It is based mainly on the articles [29, 31].


## 3.3. On the Rate of Convergence Estimates of the Conjugate Gradient Method

LetA be symmetric, positive semidefinite and consider the solution of Ax=b by a preconditioned conjugate gradient method. In order to understand how an efficient preconditioner to A should be chosen we must first understand some general properties of the rate of convergence of conjugate gradient methods.


### 3.3.1. Rate of Convergence Estimates Based on  Minimax Approximation

As is well known (see e.g., [6, 30]), the conjugate gradient method is a norm minimizing method. For the preconditioned standard CG method, we have(45)‖ek‖A=min⁡Pk∈πk‖Pk(B)e0‖A,
where ∥u∥A={uTAu}1/2, ek=x-xk is the iteration error and πk denotes the set of polynomials of degree k which are normalized at the origin, that is, Pk(0)=1. This is a norm on the subspace orthogonal to the mullspace of A, that is, on the whole space, if A is nonsingular.Consider theC-innerproduct (u,v)=uTCv, and note that B=C-1A is symmetric with respect to this innerproduct, let v1,v2,…,vn be orthonormal eigenvectors and let λi,  i=1,…,n be the corresponding eigenvalues of B. Let(46)e0=∑j=1nαjvj,be the eigenvector expansion of the initial vector where αj=(e0,vi),  i=1,…,n. Note further that the eigenvectors are both A- and C-orthogonal. Then, by the construction of the CG method, it follows(47)ek=∑j=1nαjPk(λj)vj,and, using the nonnegativity of the eigenvalues, we find(48)‖ek‖A=‖∑j=1nαjPk(λj)vj‖A={∑j=1nαj2λjPk2(λj)}1/2≤{∑1,λj>0nαj2λj}1/2max⁡1≤i≤nλi>0|Pk(λi)|=max⁡1≤i≤nλi>0|Pk(λi)|‖e0‖A.
Here we have used the A-orthogonality of the eigenvectors. Similarly, using the C-orthogonality, we find(49)‖ek‖C={∑j=1nαj2Pk2(λj)}1/2≤max⁡1≤i≤n|Pk(λj)|‖e0‖C.Due to the minimization property (45) there follows from (48) the familiar bound(50)‖ek‖A≤min⁡Pk∈πk1max⁡1≤i≤nλi>0|Pk(λi)|‖e0‖A.
Estimate (50) is sharp in the respect that for every k there exists an initial vector for which equality is attained. In fact, for such a vector we necessarily have that αj≠0 if and only if αj belongs to a set of k+1 points ( the so-called Haar condition) where maxi|Pk(λi)| is taken. For such an initial vector (49) shows that, if the eigenvalues are positive, we have also(51)‖ek‖C=min⁡Pk∈πk1max⁡1≤i≤n|Pk(λi)|‖e0‖C.The rate of convergence of the iteration error∥ek∥A is measured by the average convergence factor(52){‖ek‖A‖e0‖A}1/k.Inequality (50) shows that this can be majorized with an estimate of the rate of convergence of a best polynomial approximation problem (namely the best approximation of the function ≡0, of polynomials in πk1) in maximum norm on the discrete set formed by the spectrum of B. Clearly, multiple eigenvalues are treated as single so the actual approximation problem is(53)min⁡Pk∈πk1max⁡1≤i≤m|Pk(λ̃i)|,
where the disjoint positive eigenvalues λ̃j have been ordered in increasing value, 0<λ̃i<⋯<λ̃m, and m is the number of such eigenvalues. However, the solution of this problem requires knowledge of the spectrum, which is not available in general. Even if it is known, the estimate (53) can be troublesome in practice, since it involves approximation on a general discrete set of points.Besides being costly to apply, such estimates do not give any qualitative insight in the behaviour of the conjugate gradient method for various typical eigenvalue distributions.That is why we make some further assumptions on the spectrum in order to simplify the approximation problem and at the same time, present estimates which can be used both to estimate the number of iterations and to give some insight in the qualitative behaviour of the iteration method.


### 3.3.2. Standard Condition Number: Linear Convergence

If the eigenvalues are (densely) located in an interval[a,b] where a>0, we can majorize the best approximation problem on the discrete set with the best approximation problem on the interval and frequently still get a good estimate. We have(54)min⁡Pk∈πk1max⁡1≤i≤m|Pk(λ̃i)|≤min⁡Pk∈πk1max⁡a≤x≤b|Pk(x)|.The solution to thisminmax problem is well known and uses Chebyshev polynomials. One finds that(55)min⁡Pk∈πk1max⁡a≤x≤b|Pk(x)|=1Tk((b+a)/(b-a))=2σk1+σ2k,
where σ=(1-a/b)/(1+a/b), and Tk(x)=(1/2)[(x+x2-1)+(x-x2-1)k],  a=λ̃1,  b=λ̃m. Hence, the average rate of convergence of the upper bound approaches σ as k→∞. Also, it is readily found (see [35]) that the relative iteration error ∥ek∥A/∥e0∥A≤ɛ if(56)k=k*(a,b,ɛ)=⌈ln⁡⁡((1/ɛ)+(1/ɛ2)-1)ln⁡⁡σ-1⌉.
Here ⌈ξ⌉ denotes the smallest integer not less than ξ.It turns out that the above holds more generally ifA is nonsymmetric but the eigenvalues are contained in an ellipse with foci a, b, where b≥a>0, if one replaces σ with σ̂=σ(1+δ)/(1-δ), where δ is the eccentricity of the ellipse, (i.e., the ratio of the semiaxes) and δ<2a/b/(1+a/b). Also, in a similar way, the case of eigenvalues contained in two separate intervals or ellipses can be analysed, see, for example, [35] for further details.Whenb/a→∞,  δ=0, and ɛ→0 the following upper bound becomes an increasingly accurate replacement of (56),(57)k*≤⌈12baln⁡⁡2ɛ⌉.
The above estimate of the rate of convergence and of the number of iterations show that they depend only on the condition number b/a and on the eccentricity of the ellipse, containing the eigenvalues. Therefore, except in special cases, this estimate is not very accurate. When we use a more detailed information of the spectrum and the initial error vector, sometimes substantially better estimates can be derived. This holds for instance when there are well separated small and/or large eigenvalues. Before we consider this important case, we mention briefly another similar minimax result which holds when we use different norms for the iteration error vector and for the initial vector.By (48), we have(58)‖ek‖A={∑αj2λjPk2(λj)}1/2={∑αj2λj1-2sλj2sPk2(λj)}1/2≤min⁡Pk∈πk1max⁡1≤λj≤m|λjsPk(λj)|{∑αj2λj(1-2s)}1/2=min⁡Pk∈πk1max⁡1≤λj≤m|λjsPk(λj)|‖e0‖A1-2s.If the initial vector is such that Fourier coefficients for the highest eigenvalue modes are dominating, then∥e0∥A1-2s may exist and take not too large values even for some s≥1/2. We consider the interesting case where s≥1/2, for which the following theorem holds (see [6, 36]).Theorem 6.
Letπk1 denote the set of polynomials of degree k such that Pk(0)=1. Then for k=1,2… and for any s≥1/2 such that 2s is an integer, it holds
(59)‖ek‖A‖e0‖A1-2s≤min⁡Pk∈πk1max⁡0≤x≤1|xsPk(x)|≤(sk+s)2s.Remark 1.
Fors=1/2 it holds
(60)max⁡0≤x≤1|x1/2Pk(x)|=12k+1,
for Pk(x)=U2k(1-x) and for s=1, it holds
(61)max⁡0≤x≤1|xPk(x)|=1k+1tanπ4k+4<1(k+1)2,
for Pk(x)=(x-1(-1)k)/(k+1)tan(π/(4k+4))Tk+1((1+cos(π/(2k+2)))x-cos(π/(2k+2))) where Tk(x) and Uk(x) are the Chebyshev polynomials of kth degree of the first and second kind, respectively.For other values (59) is an upper bound only, that is, not sharp. At any rate, it shows that the error ∥ek∥A converges (initially) at least as fast as (s/(k+s))2s, that is, as 1/(2k+1) for s=1/2 and as (1/(k+1))2 for s=1.Note that this convergence rate does not depend on the eigenvalues, in particular not on the spectral condition number.Conclusion.
By computing the initial approximation vector from a coarse mesh, the components fore0 for the first Fourier modes will be small and ∥e0∥A1-2s may take on values that are not very large even when s=1 or bigger. Therefore, there is an initial decay of the residual as O(k2s), independent of the condition number. Note, however, that the actual errors may not have decayed sufficiently even if the residual has.We consider now upper bound estimates which show how the convergence history may enter a superlinear rate.An Estimate to Show a Superlinear Convergence Rate Based on the K-Condition Number
A somewhat rough, but simple and illustrative superlinear convergence estimate can be obtained in terms of the so calledK-condition number, (see [37, 38])
(62)K=K(B)=((1/n)tr⁡(B))ndet⁡(B)=(1n∑i=1nλi)n(Πi=1nλi)-1,
where we assume that B is s.p.d.
Note thatK1/n equals the quotient between the arithmetic and geometric averages of the eigenvalues. This quantity is similar to the spectral condition number κ(B) in that it is never smaller than 1, and is equal to 1 if and only if B=αI,  α>0 (recall that B is symmetrizable).
Based on theK-condition number, a superlinear convergence result can be obtained as follows.Theorem 7.
Letk<n be even and k≥3lnK. Then,
(63)∥ek∥A∥e0∥A≤(3ln⁡Kk)k/2.Proof.
Letk=2m and the polynomial Pk be of a simplest possible form, that is, let it vanish at the m smallest and m largest eigenvalues of B. As follows from (48), we have then
(64)‖ek‖A‖e0‖A≤max⁡λ1≤λ≤λn|Πi=1m(1-λλi)(1-λλn+1-i)|=Πi=1mmax⁡λi≤λ≤λn+1-i(λλi-1)(1-λλn+1-i)=Πi=1m((λi+λn+1-i)24λiλn+1-i-1)≤((Πi=1m(λi+λn+1-i)24λiλn+1-i)1/m-1)m.
The latter follows from (Πi=1m(1-Θi))1/m+(Πi=1mΘi)1/m≤1 with Θi=4λiλn+1-i(λi+λn+1-i)2. Using now twice the inequality between the arithmetic and geometric mean values, one has
(65)Πi=1m(λi+λn+1-i)24λiλn+1-i≤((1/n)-2m∑i=m+1n-mλi)n-2mΠi=m+1n-mλi×Πi=1m(λi+λn+1-i/2)2Πi=1mλiλn+1-i≤(1/n(∑i=m+1n-mλi+∑i=1m(λi+λn+1-i)))nΠi=m+1n-mλiΠi=1mλiλn+1-i=K.
Using exp (2x)-1≤2x/(1-x),  x<1, we get the required estimate,
(66)‖ek‖A∥e0∥A≤(K2/k-1)k/2≤(2ln⁡Kk-ln⁡K)k/2≤(2ln⁡K2k/3(k/3-ln⁡K))k/2≤(3ln⁡Kk)k/2.A somewhat better result of the same type exists. The estimate is of similar type, that is,(67)∥rk∥C-1∥r0∥C-1≤(K1/k-1)k/2,
where rk=Aek was obtained using more complicated techniques, see [6, 37], and the references quoted therein. Note that here as k/lnK→∞, we have(68)∥ek∥C-1∥e0∥C-1≤(e1/kln⁡K-1)k/2≈(1+ln⁡K2k)k/2(ln⁡Kk)k/2≈K1/4(ln⁡Kk)k/2,that is, the simpler upper bound in (63) is asymptotically worse than this bound (albeit in a different norm) by the factor 3k/2/lnK1/4.The upper bounds in the above estimates involve a convergence factor which decreases with increasing iteration number and show hence a superlinear rate of convergence. Note, however, thatK1/n≤κ(B)≲4K (see [6]) where κ(B)=λn/λ1 is the spectral condition number, so the K-condition number may take very large values when κ(B) is large.The estimates based on theK-condition number involve only “integral” characteristics of the preconditioned matrix (the trace and the determinant). Sometimes, it is possible to obtain a practical estimate of K(B) which can be useful for the a priori construction of good preconditioners and for the a posteriori assessment of their quality, see Section 5 for further details.The estimate(69)‖rk‖C-1‖r0‖C-1≤(K1/k-1)k/2≤ɛshows that(70)K1/2(1-K-1/k)k/2≤ɛor(71)k2log⁡211-K-1/k≥log⁡2K1/2ɛ,which holds if(72)k>log⁡2K+2log⁡21ɛ.
Hence, when K≫ɛ-2 the estimated number of iterations depends essentially only on log2K, that is, depends little on the relative accuracy ɛ, which indicates a fast superlinear convergence, when k>log2K.When actually estimating the number of iterations, Theorem6 shows a useful result only when k>O(lnK)=n(lnK1/n), that is, the quotient between the arithmetic and geometric averages of the eigenvalues, which equals K1/n, must be close to unity and the eigenvalues must be very well clustured so K1/n=1+O(n-ɛ) for some ɛ>0; otherwise the estimated number of iterations will be O(n), which is normally a useless result. The next example illustrates this further.Example 1.
Consider a geometric distribution of eigenvalues ofA,   λj=js,   j=1,2,…,n for some positive s. Here, asymptotically
(73)tr⁡(A)=∑1njs~1s+1ns+1,n⟶∞.
Using Stirling's formula, we find
(74)det⁡(A)=∏1nλj=(∏1nj)s~(2πn)s/2(ne)ns,n⟶∞,
so,
(75)κ(A)1/n~ess+1,n⟶∞.
Hence, s must be sufficiently small for the estimate in Theorem 6. to be useful. On the other hand, the spectral condition number (i) κ(A)=ns, and the simple estimate based on κ(A) leads to k~O(ns/2) and gives hence, asymptotically, a smaller upper bound when s<2. For further discussions on superlinear rate of convergence, see [39].


## 3.3.1. Rate of Convergence Estimates Based on  Minimax Approximation

As is well known (see e.g., [6, 30]), the conjugate gradient method is a norm minimizing method. For the preconditioned standard CG method, we have(45)‖ek‖A=min⁡Pk∈πk‖Pk(B)e0‖A,
where ∥u∥A={uTAu}1/2, ek=x-xk is the iteration error and πk denotes the set of polynomials of degree k which are normalized at the origin, that is, Pk(0)=1. This is a norm on the subspace orthogonal to the mullspace of A, that is, on the whole space, if A is nonsingular.Consider theC-innerproduct (u,v)=uTCv, and note that B=C-1A is symmetric with respect to this innerproduct, let v1,v2,…,vn be orthonormal eigenvectors and let λi,  i=1,…,n be the corresponding eigenvalues of B. Let(46)e0=∑j=1nαjvj,be the eigenvector expansion of the initial vector where αj=(e0,vi),  i=1,…,n. Note further that the eigenvectors are both A- and C-orthogonal. Then, by the construction of the CG method, it follows(47)ek=∑j=1nαjPk(λj)vj,and, using the nonnegativity of the eigenvalues, we find(48)‖ek‖A=‖∑j=1nαjPk(λj)vj‖A={∑j=1nαj2λjPk2(λj)}1/2≤{∑1,λj>0nαj2λj}1/2max⁡1≤i≤nλi>0|Pk(λi)|=max⁡1≤i≤nλi>0|Pk(λi)|‖e0‖A.
Here we have used the A-orthogonality of the eigenvectors. Similarly, using the C-orthogonality, we find(49)‖ek‖C={∑j=1nαj2Pk2(λj)}1/2≤max⁡1≤i≤n|Pk(λj)|‖e0‖C.Due to the minimization property (45) there follows from (48) the familiar bound(50)‖ek‖A≤min⁡Pk∈πk1max⁡1≤i≤nλi>0|Pk(λi)|‖e0‖A.
Estimate (50) is sharp in the respect that for every k there exists an initial vector for which equality is attained. In fact, for such a vector we necessarily have that αj≠0 if and only if αj belongs to a set of k+1 points ( the so-called Haar condition) where maxi|Pk(λi)| is taken. For such an initial vector (49) shows that, if the eigenvalues are positive, we have also(51)‖ek‖C=min⁡Pk∈πk1max⁡1≤i≤n|Pk(λi)|‖e0‖C.The rate of convergence of the iteration error∥ek∥A is measured by the average convergence factor(52){‖ek‖A‖e0‖A}1/k.Inequality (50) shows that this can be majorized with an estimate of the rate of convergence of a best polynomial approximation problem (namely the best approximation of the function ≡0, of polynomials in πk1) in maximum norm on the discrete set formed by the spectrum of B. Clearly, multiple eigenvalues are treated as single so the actual approximation problem is(53)min⁡Pk∈πk1max⁡1≤i≤m|Pk(λ̃i)|,
where the disjoint positive eigenvalues λ̃j have been ordered in increasing value, 0<λ̃i<⋯<λ̃m, and m is the number of such eigenvalues. However, the solution of this problem requires knowledge of the spectrum, which is not available in general. Even if it is known, the estimate (53) can be troublesome in practice, since it involves approximation on a general discrete set of points.Besides being costly to apply, such estimates do not give any qualitative insight in the behaviour of the conjugate gradient method for various typical eigenvalue distributions.That is why we make some further assumptions on the spectrum in order to simplify the approximation problem and at the same time, present estimates which can be used both to estimate the number of iterations and to give some insight in the qualitative behaviour of the iteration method.


## 3.3.2. Standard Condition Number: Linear Convergence

If the eigenvalues are (densely) located in an interval[a,b] where a>0, we can majorize the best approximation problem on the discrete set with the best approximation problem on the interval and frequently still get a good estimate. We have(54)min⁡Pk∈πk1max⁡1≤i≤m|Pk(λ̃i)|≤min⁡Pk∈πk1max⁡a≤x≤b|Pk(x)|.The solution to thisminmax problem is well known and uses Chebyshev polynomials. One finds that(55)min⁡Pk∈πk1max⁡a≤x≤b|Pk(x)|=1Tk((b+a)/(b-a))=2σk1+σ2k,
where σ=(1-a/b)/(1+a/b), and Tk(x)=(1/2)[(x+x2-1)+(x-x2-1)k],  a=λ̃1,  b=λ̃m. Hence, the average rate of convergence of the upper bound approaches σ as k→∞. Also, it is readily found (see [35]) that the relative iteration error ∥ek∥A/∥e0∥A≤ɛ if(56)k=k*(a,b,ɛ)=⌈ln⁡⁡((1/ɛ)+(1/ɛ2)-1)ln⁡⁡σ-1⌉.
Here ⌈ξ⌉ denotes the smallest integer not less than ξ.It turns out that the above holds more generally ifA is nonsymmetric but the eigenvalues are contained in an ellipse with foci a, b, where b≥a>0, if one replaces σ with σ̂=σ(1+δ)/(1-δ), where δ is the eccentricity of the ellipse, (i.e., the ratio of the semiaxes) and δ<2a/b/(1+a/b). Also, in a similar way, the case of eigenvalues contained in two separate intervals or ellipses can be analysed, see, for example, [35] for further details.Whenb/a→∞,  δ=0, and ɛ→0 the following upper bound becomes an increasingly accurate replacement of (56),(57)k*≤⌈12baln⁡⁡2ɛ⌉.
The above estimate of the rate of convergence and of the number of iterations show that they depend only on the condition number b/a and on the eccentricity of the ellipse, containing the eigenvalues. Therefore, except in special cases, this estimate is not very accurate. When we use a more detailed information of the spectrum and the initial error vector, sometimes substantially better estimates can be derived. This holds for instance when there are well separated small and/or large eigenvalues. Before we consider this important case, we mention briefly another similar minimax result which holds when we use different norms for the iteration error vector and for the initial vector.By (48), we have(58)‖ek‖A={∑αj2λjPk2(λj)}1/2={∑αj2λj1-2sλj2sPk2(λj)}1/2≤min⁡Pk∈πk1max⁡1≤λj≤m|λjsPk(λj)|{∑αj2λj(1-2s)}1/2=min⁡Pk∈πk1max⁡1≤λj≤m|λjsPk(λj)|‖e0‖A1-2s.If the initial vector is such that Fourier coefficients for the highest eigenvalue modes are dominating, then∥e0∥A1-2s may exist and take not too large values even for some s≥1/2. We consider the interesting case where s≥1/2, for which the following theorem holds (see [6, 36]).Theorem 6.
Letπk1 denote the set of polynomials of degree k such that Pk(0)=1. Then for k=1,2… and for any s≥1/2 such that 2s is an integer, it holds
(59)‖ek‖A‖e0‖A1-2s≤min⁡Pk∈πk1max⁡0≤x≤1|xsPk(x)|≤(sk+s)2s.Remark 1.
Fors=1/2 it holds
(60)max⁡0≤x≤1|x1/2Pk(x)|=12k+1,
for Pk(x)=U2k(1-x) and for s=1, it holds
(61)max⁡0≤x≤1|xPk(x)|=1k+1tanπ4k+4<1(k+1)2,
for Pk(x)=(x-1(-1)k)/(k+1)tan(π/(4k+4))Tk+1((1+cos(π/(2k+2)))x-cos(π/(2k+2))) where Tk(x) and Uk(x) are the Chebyshev polynomials of kth degree of the first and second kind, respectively.For other values (59) is an upper bound only, that is, not sharp. At any rate, it shows that the error ∥ek∥A converges (initially) at least as fast as (s/(k+s))2s, that is, as 1/(2k+1) for s=1/2 and as (1/(k+1))2 for s=1.Note that this convergence rate does not depend on the eigenvalues, in particular not on the spectral condition number.Conclusion.
By computing the initial approximation vector from a coarse mesh, the components fore0 for the first Fourier modes will be small and ∥e0∥A1-2s may take on values that are not very large even when s=1 or bigger. Therefore, there is an initial decay of the residual as O(k2s), independent of the condition number. Note, however, that the actual errors may not have decayed sufficiently even if the residual has.We consider now upper bound estimates which show how the convergence history may enter a superlinear rate.An Estimate to Show a Superlinear Convergence Rate Based on the K-Condition Number
A somewhat rough, but simple and illustrative superlinear convergence estimate can be obtained in terms of the so calledK-condition number, (see [37, 38])
(62)K=K(B)=((1/n)tr⁡(B))ndet⁡(B)=(1n∑i=1nλi)n(Πi=1nλi)-1,
where we assume that B is s.p.d.
Note thatK1/n equals the quotient between the arithmetic and geometric averages of the eigenvalues. This quantity is similar to the spectral condition number κ(B) in that it is never smaller than 1, and is equal to 1 if and only if B=αI,  α>0 (recall that B is symmetrizable).
Based on theK-condition number, a superlinear convergence result can be obtained as follows.Theorem 7.
Letk<n be even and k≥3lnK. Then,
(63)∥ek∥A∥e0∥A≤(3ln⁡Kk)k/2.Proof.
Letk=2m and the polynomial Pk be of a simplest possible form, that is, let it vanish at the m smallest and m largest eigenvalues of B. As follows from (48), we have then
(64)‖ek‖A‖e0‖A≤max⁡λ1≤λ≤λn|Πi=1m(1-λλi)(1-λλn+1-i)|=Πi=1mmax⁡λi≤λ≤λn+1-i(λλi-1)(1-λλn+1-i)=Πi=1m((λi+λn+1-i)24λiλn+1-i-1)≤((Πi=1m(λi+λn+1-i)24λiλn+1-i)1/m-1)m.
The latter follows from (Πi=1m(1-Θi))1/m+(Πi=1mΘi)1/m≤1 with Θi=4λiλn+1-i(λi+λn+1-i)2. Using now twice the inequality between the arithmetic and geometric mean values, one has
(65)Πi=1m(λi+λn+1-i)24λiλn+1-i≤((1/n)-2m∑i=m+1n-mλi)n-2mΠi=m+1n-mλi×Πi=1m(λi+λn+1-i/2)2Πi=1mλiλn+1-i≤(1/n(∑i=m+1n-mλi+∑i=1m(λi+λn+1-i)))nΠi=m+1n-mλiΠi=1mλiλn+1-i=K.
Using exp (2x)-1≤2x/(1-x),  x<1, we get the required estimate,
(66)‖ek‖A∥e0∥A≤(K2/k-1)k/2≤(2ln⁡Kk-ln⁡K)k/2≤(2ln⁡K2k/3(k/3-ln⁡K))k/2≤(3ln⁡Kk)k/2.A somewhat better result of the same type exists. The estimate is of similar type, that is,(67)∥rk∥C-1∥r0∥C-1≤(K1/k-1)k/2,
where rk=Aek was obtained using more complicated techniques, see [6, 37], and the references quoted therein. Note that here as k/lnK→∞, we have(68)∥ek∥C-1∥e0∥C-1≤(e1/kln⁡K-1)k/2≈(1+ln⁡K2k)k/2(ln⁡Kk)k/2≈K1/4(ln⁡Kk)k/2,that is, the simpler upper bound in (63) is asymptotically worse than this bound (albeit in a different norm) by the factor 3k/2/lnK1/4.The upper bounds in the above estimates involve a convergence factor which decreases with increasing iteration number and show hence a superlinear rate of convergence. Note, however, thatK1/n≤κ(B)≲4K (see [6]) where κ(B)=λn/λ1 is the spectral condition number, so the K-condition number may take very large values when κ(B) is large.The estimates based on theK-condition number involve only “integral” characteristics of the preconditioned matrix (the trace and the determinant). Sometimes, it is possible to obtain a practical estimate of K(B) which can be useful for the a priori construction of good preconditioners and for the a posteriori assessment of their quality, see Section 5 for further details.The estimate(69)‖rk‖C-1‖r0‖C-1≤(K1/k-1)k/2≤ɛshows that(70)K1/2(1-K-1/k)k/2≤ɛor(71)k2log⁡211-K-1/k≥log⁡2K1/2ɛ,which holds if(72)k>log⁡2K+2log⁡21ɛ.
Hence, when K≫ɛ-2 the estimated number of iterations depends essentially only on log2K, that is, depends little on the relative accuracy ɛ, which indicates a fast superlinear convergence, when k>log2K.When actually estimating the number of iterations, Theorem6 shows a useful result only when k>O(lnK)=n(lnK1/n), that is, the quotient between the arithmetic and geometric averages of the eigenvalues, which equals K1/n, must be close to unity and the eigenvalues must be very well clustured so K1/n=1+O(n-ɛ) for some ɛ>0; otherwise the estimated number of iterations will be O(n), which is normally a useless result. The next example illustrates this further.Example 1.
Consider a geometric distribution of eigenvalues ofA,   λj=js,   j=1,2,…,n for some positive s. Here, asymptotically
(73)tr⁡(A)=∑1njs~1s+1ns+1,n⟶∞.
Using Stirling's formula, we find
(74)det⁡(A)=∏1nλj=(∏1nj)s~(2πn)s/2(ne)ns,n⟶∞,
so,
(75)κ(A)1/n~ess+1,n⟶∞.
Hence, s must be sufficiently small for the estimate in Theorem 6. to be useful. On the other hand, the spectral condition number (i) κ(A)=ns, and the simple estimate based on κ(A) leads to k~O(ns/2) and gives hence, asymptotically, a smaller upper bound when s<2. For further discussions on superlinear rate of convergence, see [39].


## 3.4. Generalized Conjugate Gradient Methods

The rate of convergence estimates as given above, holds for a restricted class of matrices, symmetric or, more generally, for normal matrices.To handle more general classes of problems for which such optimal rate of convergence results as in (45) holds, one needs more involved methods. Much work has been devoted to this problem. This includes methods like generalized minimum residual (GMRES), see Saad and Schultz [40], generalized conjugate residual (GCR), and generalized conjugate gradient (GCG), see [6] and for further details [41]. As opposed to the standard conjugate and gradient method, they require a long version of updates for the search directions, as the newest search direction at each stage is not in general, automatically (in exact precision) orthogonal to the previous search directions, but must be orthogonalized at each step. This makes the computational expense per step grow linearly and the total expense grows quadratically with the iteration index. In addition, due to finite precision, there is a tendency of loss of orthogonality, even for symmetric problems when many iterations are required. One remedy which has been suggested is to use the method only for a few steps, say 10, and restart the method with the current approximation as initial approximation.Clearly, however, in this way, the optimal convergence property of the whole Krylov set of vector is lost. For this and other possibilities, see, for example, [42].Another important version of the generalized conjugate gradient methods occurs when one uses variable preconditioners. Variable preconditioners, that is, a preconditioner changed from one iteration to the next iteration step, are used in many contexts.For instance, one can use variable drop tolerance, computed adaptively, in an incomplete factorization method (see Section4). When the given matrix is partitioned in two by two blocks, it can be efficient to use inner iterations when solving arising systems for one, or both, of the diagonal block matrices, see, for example, [43], and the flexible conjugate gradient method in Saad, [44, 45].Due to space limitations, the above topics will not be further discussed in this paper.


## 4. Incomplete Factorization Methods

There exist two classes of preconditioning methods that are closely related to direct solution methods. In this paper, we make a survey only of their main ingredients, but delete many of the particular aspects.The first method is based on incomplete factorization were some entries arising during a triangular factorization are neglected to save in memory. The deletion can be based on some drop tolerance criterion or on a normally a priori, chosen sparsity pattern. The factorization based on a drop tolerance takes the following form. During the elimination (or equivalently, triangular factorization), the off-diagonal entries are accepted only if they are not too small. For instance,(76)aij:={aij-airarr-1arjif|aij|≥εaiiajj,0,otherwise.
Here, ε,  0<ε≪1 is the drop-tolerance parameter. Such methods may lead to too much fill-in (i.e., aij≠0 in positions where the original entry was occupied by a zero), because to be robust, they may require near machine-precision drop tolerances. Furthermore, as direct solution methods, they are difficult to parallelize efficiently.The incomplete factorization method can readily be extended to matrices partitioned in block form. Often, instead of a drop tolerance, one prescribes the sparsity pattern of the triangular factors in the computed preconditioner, that is, entries arising outside the chosen pattern are ignored. An early presentation of such incomplete factorization methods was given by Meijerink and van der Vorst [46]. One can make a diagonal compensation of the neglected entries, that is add them to the diagonal entries in the same row, possibly first multiplied by some scalar Θ,0<Θ≤1. For discussions of such approaches, see [29, 30, 47, 48]. This frequently moves small eigenvalues, corresponding to the smoother harmonics, to cluster near the origin, in this way sometimes improving the spectral condition number by an order of magnitude (see [6, 47]).The other class of methods are based on approximate inversesG, for instance such that minimizes a Frobenius norm of the error matrix I-GA, see Section 5 for further details. To be sufficiently accurate these methods lead frequently to nearly full matrices. This can be understood as the matrices we want to approximate are often sparse discretizations of diffusion problems. The inverse of such an operator is a discrete Green's function which, as wellknown, often has a significantly sized support on nearly the whole domain of definition.However, we can use an additive approximation of the inverse involving two, or more, terms which is approximate on different vector subspaces. By defining in this way the preconditioner recursively on a sequence of lower dimensional subspaces, it may preserve the accurate approximation property of the full, inverse method while still needing only actions of sparse operators.Frequently, the given matrices are partitioned in a natural way in a two by two block form. For such matrices, it can be seen that the two approaches are similar. Consider namely(77)A=[A1A12A21A2],
where we assume that A1 and the Schur complement matrix S=A2-A21A1-1A12 are nonsingular. (This holds, in particular, if A is symmetric and positive definite.) We can construct either a block approximate factorization of A or approximate the inverse of A on additive form. As the following shows, the approaches are related. First, a block matrix factorization of A is(78)A=[A10A21S][I1A1-1A120I2],
where I1,  I2 denote the unit matrices of proper order. For its inverse, it holds(79)A-1=[I1-A1-1A120I2][A1-10-S-1A21A1-1S-1]=[A1-1+A1-1A12S-1A21A1-1-A1-1A12S-1-S-1A21A1-1S-1],
or(80)A-1=[A1-1000]+[-A1-1A12I2]S-1[-A21A1-1,I2].
A straightforward computation reveals that AṼ≡ṼTAṼ=S and, hence,(81)A-1=[A1-1000]+Ṽ(ṼTAṼ)-1ṼT=[A1-1000]+ṼAṼ-1ṼT,
where(82)Ṽ=[-A1-1A12I2].LetM1≃A1 be an approximation of A1 ( for which linear systems are simpler to solve than for A1) and let G1≃A1-1 be a sparse approximate inverse. Possibly G1=M1-1. Then,(83)M=[M10A21B2][I1M1-1A120I2]=[M1A12A21A2]+[000B2+A21M1-1A12-A2]
is a preconditioner to Aand(84)B=[G1000]+VB2-1VT
is an approximate inverse, where V=[-G1A12I2] and B2 is an approximation of S. If B2=VTÃV, where Ã=[G1-1A12A21A2], then(85)B=[G1000]+V(VTÃV)-1VT=[G1000]+VS(Ã)VT,
where S(Ã)=A2-A21G1A12. If M1=G1-1, then in this case(86)B=M-1.
Hence, a convergence estimate for one method can be directly applied for the other method as well. For further discussions of block matrix preconditioners, see, for example, [49–52]. As can be seen from the above, Schur complement matrices play a major role in both matrix factorizations. For sparse approximations of  Schur complement matrix, in particular element, e.g., element type approximations, see, for example [53–55].We consider now multilevel extensions of the additive approximate inverse subspace correction method. It is illustrative to consider first the exact inverse and its relation to Gaussian (block matrix) elimination.


### 4.1. The Exact Inverse on Additive Form

Let thenA(0)=A and consider a matrix(87)A(k)=[A1(k)A12(k)A21(k)A2(k)],
in a sequence defined by(88)A(k+1)≡S2(k)=A2(k)-A21(k)A1(k)-1A21(k),k=0,1,…,k0,
where each A1(k) in nonsingular, being a block diagonal of a symmetric positive definite matrix. Hence, the following recursion holds(89)A(k)-1=[A1(k)-1000]+[-A1(k)-1A12(k)I2(k+1)]A(k+1)-1[A21(k)A1(k)-1,I2(k+1)],k=0,  1,…,k0.Here, I2(k+1) is the identity matrix on level k+1. Note that in this example the dimensions decrease with increasing level number and the final matrix (i.e., Schur complement) in the sequence is A(k0)=S2(k0+1). The above recursion can be rewritten in compact form(90)A-1=[A1(0)-1000]+UD-1L,
where the kth column of the block upper triangular matrix U equals [-A1(k)-1A12(k)I2(k+1)] and L=UT. Further,(91)D=[A1(k)0A(2)⋱0A(k0)].Hence, this is the (block matrix) Gaussian elimination method applied directly to form the inverse matrix. In this way, there is no need to first form the factorizationA=L̃DŨ and then A-1=Ũ-1D-1L̃-1. As is wellknown and readily seen, the columns of Ũ-1 and L̃-1 are formed directly with no additional computation, from those of Ũ and L̃, respectively. Note that Ũ-1 is upper (block) triangular and L̃-1 is lower (block) triangular).The matrixD contains the (block) pivot matrices which arise during the factorization. Permutations to increase the stability by finding dominating pivots can be done by replacing A(k+1) with P(k)TÃ(k+1)P(k) where Ã(k+1)=P(k)A(k+1)P(k)T is the permuted matrix on which the next elimination step takes place.An incomplete factorization method for approximate inverses can be defined by approximating each arising Schur complement matrix with some sparse matrixB̃(k+1) and possibly also approximating A1(k)-1 with some matrix B1(k), to yield the approximate inverse(92)B(k)=[B1(k)000]+V(k)B̃(k+1)V(k)T,
where V(k)=[-B1(k)A12(k)I2(k+1)].In forming the approximate Schur complement one can use a simpler matrixD1(k)than B1(k), often a diagonal matrix suffices. The intermediate Schur complement matrix S̃2(k)=A2(k)-A21(k)D1(k)A12(k) can be possibly further approximated by deleting certain off-diagonal entries to preserve sparsity. These entries can be compensated for by modifying the diagonal of S̃2(k) to form the final approximation B̃(k+1). Thereby, it can be important to make the approximate Schur complement exact on some particular vector or vector space. (We do not go into these aspects further here, see [43, 56] for details.)The eigenvectors for the smallest eigenvalues ofA provide efficient column vectors for the matrix V to reduce significantly the condition number of BA as compared to that of A. However, in general the eigenvectors are not known, and even if they are known it would be costly to apply the corresponding preconditioner as V would be a full matrix. A more viable choice is to let V be defined by the basis functions {φi(H)} of a coarse mesh (or coarsened matrix graph) so that(93)Im⁡V=span⁡{φi(H)}.V,   VT acts then, respectively, as prolongation and restriction operators and(94)V=[J12I2],
where J12 is the interpolation matrix from the coarse mesh (ΩH) to the fine mesh (Ωh), and we assume   ΩH⊂Ωh. Further, letting the matrices be variationally defined, as in a finite element method, we have(95)AH=VTAV,
where A is the finite element matrix on the fine mesh.Now, the eigenvectors for the smaller eigenvalues ofAH are normally accurate approximations of the corresponding eigenvectors for A. Furthermore, the eigenvectors of AH are members of ImV. Therefore, the matrix V in (94) acts nearly as well as the eigenvector matrix but, in addition, is sparse. Hence the approximate inverse takes the form(96)B=G+σVAH-1VT,
where σ=λmax(GA),  AH=VTAV.Here, the projection matrix(97)P=VAH-1VTA,
projects vectors on the subspace ImV, containing normally good approximations of the eigenvectors for the smallest eigenvalues of A, that is, those who may cause severe ill-conditioning. Clearly, the approximation is more accurate as closer ΩH is to Ωh. However, the cost of the action of P (mainly the coarse mesh solver for the action of AH-1) increases when ΩH expands. One can balance ΩH to Ωh in order to let the action of P involve the same order of computational complexity as an action of A, that is, O(h-2) for a sparse matrix A. Assuming that the cost of an action of AH-1 is O(H-2.5) in a 2D diffusion problem (e.g., using a modified incomplete factorization method as preconditioner for the conjugate gradient method), we find H=h4/5. The number of outer iterations with preconditioner B depends also on the choice of G. We refer the discussion of how G can be chosen properly to [56].As an example, for a model diffusion problem with constant coefficients on a regular mesh, say for the Laplacian operator on unit square, the eigenvectors for the Laplacian(-Δ) on Ωh with Dirichlet boundary conditions are(98)vk,l(h)=sin⁡kπsin⁡lπy,x,y∈Ωh,
where k,  l=1,2,…,h-1-1, for the eigenvalues(99)λk,l(h)=(2sin⁡kπh2)2+(2sin⁡lπh2)2.
For the coarse mesh, it holds(100)vk,l(H)=sin⁡kπxsin⁡lπy,x,y∈ΩH,
where k,  l,=1,2,…,H-1, and here Vvk,l(H) are good approximations (interpolants) of the eigenvectors vk,l(h) on Ωh for the smallest eigenvalues.An alternative choice of matrixV is to take eigenvectors from a nearby problem, normally defined by taking limit values of some problem parameter, see [56].Multigrid, algebraic multilevel and algebraic multigrid methods have been presented thoroughly in, for  example [29, 43, 57, 58]. Because of space limitations, they can not be presented here.


### 4.2. Symmetrization of Preconditioners; the SSOR and ADI Methods

As we have seen, the incomplete factorization methods require first a factorization step. There exists simpler preconditioning methods that require no factorization but have a form similar to the incomplete factorization methods. We will present two methods of this type. As an introduction, consider first an iterative method of the form(101)M(xl+1-xl)=b-Axl,l=0,1,…,
to solve Ax=b, where A and M are nonsingular. As we saw in Section 2, the asymptotic rate of convergence is determined by the spectral radius of the iteration matrix(102)B=I-M-1A.
For a method such as the SOR method (which also requires no factorization), with optimal overrelaxation parameter ω (assuming that A has property Aπ or A is s.p.d., see Section 2), the eigenvalues of the corresponding iteration matrix B are situated on a circle. No further acceleration is then possible.There is, however, a simple remedy to this, based on taking a step in the forward direction of the chosen ordering, followed by a backward step, that is, a step in the opposite order to the vector components. This method is said to have its origin in the early days of computers when programs were stored on tapes that had to be rewound before a new forward SOR step could begin. It was found that this otherwise useless computer time for the rewinding could be better used for a backward SOR sweep!As we will see, for symmetric and positive definite matrices the combined forward and backward sweeps correspond to a s.p.d. matrix which, contrary to the SOR method, has the advantage that it can be used as a preconditioning matrix in an iterative acceleration method. This method, called the SSOR method, will be defined later.For an early discussion of the SSOR method, used as a preconditioner, see [59]. For discussions about symmetrization of preconditioners, see [6, 60, 61]. More generally, if A is s.p.d, we consider the symmetrization of an iterative method in the form(103)xl+1=xl+M-1(b-Axl).
For the analysis only, we consider the transformed form of (103),(104)yl+1=(I-A1/2M-1A1/2)yl+b̃,
where(105)yl=A1/2xl,b̃=A1/2M-1b̃.
If M is unsymmetric, the iteration matrix I-A1/2M-1A1/2 is also unsymmetric. We will now consider a method using M and another preconditioner chosen so that the iteration matrix for the combined method becomes symmetric. We call this the symmetrization of the method.LetM1,  M2 be two such preconditioning matrices. Let(106)Bi=I-M̃i-1,M̃i=A-1/2MiA-1/2,
and consider the combined iteration matrix B2B1. As we will now see, it arises as an iteration matrix for the combined method(107)M1(xl+1/2-xl)=b-Axl,M2(xl+1-xl+1/2)=b-Axl+1/2,l=0,1,….
For the analysis only, we transform this to the form(108)yl+1/2-yl=b̃(1)-M̃1-1yl,yl+1-yl+1/2=b̃(2)-M̃2-1yl+1/2,
where(109)b̃(i)=A(1/2)Mi-1b.
This iteration takes the form(110)yl+1=b̃(2)+(I-M̃2-1)b̃(1)+(I-M̃1-1yl),
that is,(111)yl+1=b̃(2)+(I-M̃2-1)b̃(1)+(I-M̃2-1)(I-M̃1-1)yl,
or(112)yl+1=b̂+B2B1yl,l=0,1,…,
where(113)b̂=b̃(2)+(I-M̃2-1)b(1).
For the following we need a lemma.Lemma 2.
IfA,   B, and C are Hermitian positive definite and each pair of them commute, then ABC is Hermitian positive definite.Proof.
We have(ABC)*=CBA and use commutativity to find
(114)CBA=BCA+BAC=ABC.
Hence, ABC is Hermitian. Next, we show that the product of two s.p.d matrices that commute is positive definite. We have
(115)A-1/2ABA1/2ABA1/2=A1/2BA1/2,
which is Hermitian positive definite. Hence, by similarity, the eigenvalues of AB are positive and, since
(116)(AB)*=AB,AB is Hermitian positive definite. In the same way, (AB)C is Hermitian positive definite.Lemma 3.
LetA be s.p.d. and assume either of the following additional conditions:(a)
M2*=M1.(b)
M1,  M2 are s.p.d. ρ(A1/2Mi-1A1/2)<1,   i=1,2, and the pair of matrices M1, M2, commutes.Then, the combined iteration method (107) converges if and only if M1+M2-A is s.p.d.Proof.
It is readily seen thaty=b̂+B2B1y (i.e., the iteration method is consistent with Ax=b), where y=A1/2x and x=A-1b. Hence,
(117)y-yl+1=B2B1(y-yl),
and the iteration method (112), and hence (107) converges for any initial vector if and only if ρ(B2B1)<1, where ρ(·) denotes the spectral radius. But
(118)B2B1=I-M̃1-1-M̃2-1+M̃1-1M̃2-1=I-A1/2M1-1(M1+M2-A)M2-1A1/2.
It is readily seen that under either of the given conditions (a) or (b),
(119)M1-1(M1+M2-A)M2-1=M1-1+M2-1-M1-1AM2-1
is symmetric. Further, it is positive definite if and only if M1+M2-A is positive definite. Hence, I-B2B1 is s.p.d. Further, B2B1=(I-M̃2-1)(I-M̃1-1) is symmetric, and a similarity transformation shows that B2B1 is similar to (I-M̃2-1)1/2(I-M̃1-1)(I-M̃2-1)1/2, which is a congruence transformation of I-M̃1-1, whose eigenvalues are positive. Hence, B2B1 has positive eigenvalues, so the eigenvalues of B2B1 are contained in the interval (0,1) and, in particular, ρ(B2B1)<1.The proof of Lemma3 shows that B2B1 is symmetric, so the combined iteration method is a symmetrized version of either of the simple methods.Let us now consider a special class of symmetrized methods. We letA be split as A=D+L+U, where we assume that D is s.p.d., and let(120)V=(1-1ω)D+L,H=(1-1ω)D+U,D̂=(2/ω-1)D, where ω is a parameter, 0<ω<2. (Here, L and U are not necessarily the lower and upper triangular parts of A.) Note that(121)D̂+V+H=A,
so this is also a splitting of A. As an example of a combined, or symmetrized, iteration method, we consider the preconditioning matrix(122)C=(D̂+V)D̂-1(D̂+H),
and show that this leads to a convergent iteration method(123)C(xl+1-xl)=b-Axl,l=0,1,….
This corresponds to choosing M1=D̂-1/2(D̂+H) and M2=(D̂+V)D̂-1/2, and it can be seen that the conditions of Lemma 3 hold if the conditions in the next theorem hold.Theorem 8.
LetA=D+L+U, where D is s.p.d. Let V,H,D̂ be defined by (120), and assume that either (a) or (b) holds, where(a)
U=L*(b)
L,U are s.p.d. and each pair of matrices L, U, D commute. Then the eigenvalues λ of the matrix C-1A, where C is defined in (122), are contained in the interval 0<λ≤1.Proof.
This can be shown either by verifying the conditions in Lemma3 or more directly as follows. As in the proof of Lemma 3, it follows that C is s.p.d. Hence, the eigenvalues of C-1A are positive. Further,
(124)C=D̂+V+H+VD̂-1H,
so, by (121)
(125)C=A+VD̂-1H.
Under either condition (a) or (b), C=VD̂-1H is symmetric and positive semidefinite.
This shows thatx*Cx≥x*Ax for all  x, so the eigenvalues of C-1A are bounded above by 1.We will now show that the matrixC can also efficiently be used as a preconditioning matrix, which for a proper value of the parameter ω, and under an additional condition, can even reduce the order of magnitude of the condition number. In this respect, note that when C is used as a preconditioning matrix for the Chebyshev iterative method, it is not necessary to have C scaled so that λ(C-1A)≤1, because it is suffices then that 0<m≤λ(C-1A)≤M, for some numbers m, M. Hence, the factor 2/ω-1 in D̂-1 can be neglected.Theorem 9.
LetA=D+L+U be a splitting of A, where A and D are s.p.d. and either (a) U=L* or (b) L,  U are s.p.d. and each pair of D,  L,  U commute. Then, the eigenvalues of matrix C-1A, where
(126)C=(1ωD+L)D̂-1(1ωD+U)
and 0<ω<2,  D̂=(2/ω-1)D, are contained in the interval
(127)[(2-ω){1+ω(1/ω-1/2)2δ-1+ωγ},1],
where
(128)δ=min⁡x≠0xTAxxTDxγ=max⁡x≠0xT(LD-1U-1/4D)xxTAx.
Further, if there exists a vector for whichxT(L+U)xT(L+U)x≤0, then γ≥-1/4, and if(129)ρ(L̃Ũ)≤14,
then γ≤0, and if
(130)ρ(L̃Ũ)≤14+O(δ),thenγ≤O(1),δ⟶0.
Here, L̃=D-1/2LD-1/2.Proof.
It is readily seen that(131)C=12-ω(1ωD+L)(1ωD)-1(1ωD+U)=12-ω[A+(1ω-1)D+ωLD-1U]=12-ω[A+ω(1ω-12)2D+ω(LD-1U-14D)].
This shows the lower bound in (127); the upper bound follows by Theorem 8. By choosing a vector for which xT(L+U)x≤0, it follows that
(132)xT(LD-1U-(1/4)D)xxTAx≥xT(LD-1x-(1/4)Dx)xTDx≥-14,
which shows γ≥-1/4. The remainder of the theorem is immediate.


#### 4.2.1. The Condition Number

Theorem9 shows that the optimal value of ω to minimize the upper bound of the condition number of C-1A is the value that minimizes the real-valued  function (133)f(ω)=1+ω(1/ω-1/2)2δ-1+ωγ2-ω.
It is readily seen (see Axelsson and Barker, 1984 [30]), that f(ω) is minimized for(134)ω*=21+2(1/2+γ)δ,min⁡ωf(ω)=f(ω*)=(12+γ)δ-1+12.
In general, δ is not known, but we may know that δ=O(h2), for some problem parameter, h→0 (such as for the step length in second-order elliptic problems). Then, if γ=O(1),   h→0, we let ω=2/(1+ξh) for some ξ>0, in which case(135)f(ω)=O(h-1)=O(δ-1),h⟶0.This means thatC-1A has an order of magnitude smaller condition number than A itself, which latter is O(δ-1).We consider now two applications of Theorem9.


#### 4.2.2. The SSOR Method

In the first case,L is the lower triangular part of A or the lower block triangular part, if A is partitioned in block matrix form and U=L*. Then,(136)C=12-ω(1ωD+L)(1ωD)-1(1ωD+L*),
is a symmetrized version of the SOR method and is called the SSOR (symmetric successive overrelaxation) method.As an example, for an elliptic differential equation of second order it can be seen that the conditionρ(L̃L̃T)≤1/4 holds for problems with Dirichlet boundary conditions and constant coefficients. For extensions of this, see Axelsson and Barker [30]. For the model difference equation on a square domain with side π, we have(137)δ=2(sin⁡h2)2,γ≤0,
and we find(138)ω*=21+2sin⁡h/2~21+h,f(ω*)=12δ+12~h-1+12,h⟶0.


### 4.3. The ADI Method

In the second case of methods of (101), we let L denote the off-diagonal part of the difference operator working in the x-direction and U off-diagonal part of the difference operator in the y-direction. D is its diagonal part. Then, the matrix(139)Ĉ=(1ωD+L)(1ωD)-1(1ωD+U),
is called an alternating direction preconditioning matrix and the corresponding iteration method is called the ADI (alternating direction iteration) method. In this method, we solve alternately one-dimensional difference equations in x-and y-directions. Much has been written on the ADI-method which was originally presented in Peaceman and Rachford [62]; see Varga [10], for an early influential presentation and Birkhoff et al. [63] and Wachspress [64], for instance.As we will see, for the model difference equations we get the same optimal value ofω as in (138). The condition γ=O(1) may be less restrictive, for the ADI-method, but the condition of commutativity is much more restrictive, as the following lemma shows.Lemma 4.
LetA,  B be two Hermitian matrices of order n. Then AB=BA if and only if A and B have a common set of orthonormal eigenvectors.Proof.
If such a common set of eigenvectors{vi} exists, then Avi=σivi,  Bvi=τivi and
(140)ABvi=σiτivii=1,2,…,n.
Since the eigenvector space of an Hermitian matrix is complete, we therefore have
(141)ABx=BAx,∀x∈Cn,
which shows that AB=BA. Conversely, suppose that AB=BA. As A is Hermitian, take U to be a unitary matrix that diagonalizes A, that is
(142)Ã=UAU*=[γ1I10γ2I2⋱0γrIr],
where γ1<γ2<⋯<γr are the distinct eigenvalues of A and Ij is the identity matrix of order nj, the multiplicity of γj. (Here A is posssibly permuted accordingly.) Let B̃=UBU* and partition B̃ corresponding to the partitioning of A, that is,
(143)B̃=[B11B12⋯B1r⋮Br1Br2⋯Brr].
Since AB=BA, we have
(144)ÃB̃=UABU*=UBAU*=B̃Ã.
Carying out the block multiplication ÃB̃=B̃Ã, we find that this, in turn, implies Bij=0,  i≠j, since γi≠γj,  i≠j. Simply stated, a (block) matrix commutes with a (block) diagonal matrix if and only if it is itself (block) diagonal. Hence, B̃ is block diagonal and each Hermitian submatrix Bi,i has ni orthonormal eigenvectors that are also eigenvectors of the submatrix γiIi of Ã. Since ∑i=1rni=n and all eigenvectors are orthonormal, A and B must have the same set of eigenvectors.For the second-order elliptic difference equation in two space dimensions, it turns out that the commutativity ofL and U essentially corresponds to the property that the original problem is separable, that is, that solutions of ℒu=f can be written in the form u=φ(x)ψ(y). This means that the coefficients a(x,y) and b(x,y) in the differential operator ∂/∂x[a(x,y)∂u/∂x]+∂/∂y[b(x,y)∂u/∂x]+c(x,y)u must satisfy a(x,y)=a(x),b(x,y)=b(y), and c(x,y)=c, a constant. Hence, if a(x,y)=b(x,y), then a(x,y)=b(x,y)=a, a constant. Furthermore, the convex closure of the meshpoints must be a rectangle with sides parallel to the coordinate axes (Varga, [10] 1962). If A=A1+A2, the ADI-method can be written in the form(145)(I+τ1A1)xl+1/2=(I-τ1A2)xl+τ1b,(I+τ2A2)xl+1=(I-τ2A1)xl+1/2+τ2b,l=0,1,….This is thePeaceman-Rachford[62]iteration method. The iteration matrix M is similar to  (146)(I-τ2A1)(I+τ1A1)-1(I-τ1A2)(I+τ2A2)-1.
When A1,  A2 are Hermitian positive definite, their eigenvalues λi(1),  λi(2) are positive, and(147)‖(I-τ2A1)(I+τ1A1)-1‖2=ρ((I-τ2A1)(I+τ1A1)-1)=max⁡i|1-τ2λi(1)1+τ1λi(1)|.
Thus,(148)ρ(M)=ρ((I-τ2A1)(I+τ1A1)-1(I-τ1A2)-1(I+τ2A2)-1)≤‖(I-τ2A1)(I+τ1A1)-1(I-τ1A2)(I+τ2A2)-1‖2≤‖(I-τ2A1)(I+τ1A1)-1‖2×‖(I-τ1A2)(I+τ2A2)-1‖2=μ(τ1,τ2),ρ(M)≤μ(τ1,τ2)=max⁡i|1-τ2λi(1)1+τ1λi(1)|max⁡i|1-τ1λi(2)1+τ2λi(2)|.Note that forτ1=τ2=τ>0,  μ(τ1,τ2)=μ(τ,τ)<1, so we have ρ(M)<1, that is convergence for any τ>0. This holds even if A1,A2 do not commute. Note also that when A1 and A2 commute, we have(149)ρ(M)=μ(τ1,τ2).Let us continue the analyses for the general case whereA1,  A2 do not necessarily commute. We want to compute the optimal values of τ1 and τ2 such that μ(τ1,τ2) is minimized. For simplicity, we assume that α, β are the same lower and upper bounds of the eigenvalues of A1 and A2, that is, 0<α≤λi(j)≤β, j=1,2. We have(150)μ(τ1,τ2)≤max⁡{|1-τ2α1+τ1α|,|1-τ2β1+τ1β|}×max⁡{|1-τ1α1+τ2α|,|1-τ1β1+τ2β|}.
We want to choose τ1,  τ2>0 such that this bound is as small as possible. Note, then, that for such values of τ1,  τ2 we must have 1-τiα>0 and 1-τiβ<0. Next note that each factor in the bound (150) is minimized when(151)1-τ2α1+τ1α=τ2β-1τ1β+1,1-τ1α1+τ2α=τ1β-1τ2β+1,
respectively, that is, when(152)τ1τ2-α+β2αβ(τ1-τ2)-1αβ=0,τ1τ2+α+β2αβ(τ1-τ2)-1αβ=0,
respectively. Thus, both factors are simultaneously minimized when τ1=τ2, and then τ1=τ2=1/αβ.Theorem 10.
LetA=A1+A2, where A1,  A2, are s.p.d., and consider the Peaceman-Rachford ADI method (145) to solve Ax=b with τ1=τ2=1/αβ. The spectral radius of the corresponding iteration matrix M satisfies
(153)ρ(M)≤min⁡τ1,τ2μ(τ1,τ2)≤(1-α/β1+α/β)2~1-4αβ,
if α/β→0.Proof.
Forτ1=τ2=1/αβ, we have
(154)μ(τ1,τ2)=(1-α/β1+α/β)2.Remark 2.
For a model difference equation for a second-order elliptic differential equation problem on a square with sideπ, we have with stepsize h,
(155)α=(sin⁡(h/2)h/2)2~1,β=(cos⁡(h/2)h/2)2~4h2,h⟶0.
Then,
(156)μ(τ1,τ2)=(1-tan(h/2)1+tan(h/2))2=1-sin⁡(h)1+sin⁡(h)~1-2h,h⟶0.
Note that this is just the convergence factor we get for the SOR method with an optimal overrelaxation parameter.
Sinceρ(M)≤μ(τ1,τ2), this means that the ADI method with parameters (chosen as above) converges at least as fast as the SOR method. Note, however, that in the ADI-method we must solve two systems of equations with tridiagonal coefficient matrices (I-τAi) on each step, while the pointwise SOR method requires no solution of such systems.


### 4.4. The Commutative Case

Assume now thatA1,  A2 commute. Then, as we have seen, M is symmetric and has real eigenvalues, and we can apply the Chebyshev acceleration method. The eigenvalues of the corresponding preconditioned matrix C̃ are related to the eigenvalues of M by(157)λ(C̃)=1-λ(M).
Since -ρ(M)≤λ(M)≤ρ(M), where(158)ρ(M)=(1-α/β1+α/β)2~1-4αβ,
we have(159)4αβ~1-ρ(M)≤λ(C̃)=1+ρ(M)~2-4αβ~2,αβ⟶0.
The asymptotic rate of convergence of the Chebyshev accelerated method therefore is(160)21-ρ(M)1+ρ(M)~22(αβ)1/4,α/β⟶0.
For the model difference equation, we have the asymptotic rate of convergence(161)~2h1/2,h⟶0.


### 4.5. The Cyclically Repeated ADI Method

The real power of the ADI method is brought forth when we use a sequence of parametrsτl. Assume that A1,  A2 commute, then we choose the parameters τl cyclically. With a cycle of q parameters and with the assumption(162)0<α≤λi(j)≤β,j=1,2,
we get the iteration matrix(163)M(q)=∏p=1q(I+τpA2)-1(I-τpA1)(I+τpA1)-1(I-τpA2).
The eigenvalues of M(q) are(164)∏p=1q‍1-τpλi(1)1+τpλi(1)⋅1-τpλi(2)1+τpλi(2).
In the same way as above, ρ(M(q)) is minimized when(165)d(α,β,q)=max⁡α≤x≤β∏p=1q‍|1-τpx1+τpx|.


### 4.6. A Preconditioning Method for Complex Valued Matrices

Complex valued systems of equations arise in many applications. A commonly occuring case is the solution of a matrix polynomial equation(166)Qm(A)x=b,
where A is a real square matrix and Qm is a polynomial of degree m that has no zeroes at the eigenvalues of A. Here Qm can be factored in the product of second degree, and possibly some factors of first degree polynomials with real coefficients.The second degree polynomials can be factored in products of first degree polynomials with complex coefficients.Consider then a linear system(167)Az=b.in the form(168)(R+iS)(x+iy)=c+id,
where ℛ,   S are real matrices of order n and x,y,c,d∈ℛn.The system can be solved in complex arithmetic. However, complex arithmetic leads to heavier computational complexity and it is in general difficult to precondition complex valued matrices, as the eigenvalues of the given matrix or the preconditioned matrix can be spread in the whole complex plane and the iterative solution method may then converge too slowly.One can alternatively apply a preconditioned conjugate gradient method to the Hermitian positive definite normal matrix systemAHAu=AHb for which the eigenvalues are real. At any rate, this involves complex arithmetic that costs typically three to four times as much as corresponding real arithmetic.Complex arithmetic can be avoided by rewriting (168) in real valued form, such as(169)A(1)[xy]=[R-SSR][xy]=[cd],or(170)A(2)[x-y]=[RSS-R][x-y]=[cd].The block matrices are here real but, in general, nonsymmetric and/or indefinite. For the solution, one can use a generalized conjugate gradient method such as GMRES [40] or GCG [65].ForA(1) it holds that any eigenvalue λ appears in complex conjugate pairs λ,  λ¯. For A(2), which is real symmetric, for any eigenvalue λ≠0,-λ is also an eigenvalue. Thus, the spectrum σ(A(1)) is symmetric with respect to the real axes and the spectrum σ(A(2)) is symmetric with respect to the imaginary axes, that is, in both cases the spectrum embraces the origin. From best polynomial approximation properties it is known that such point distributions leads to polynomials of essentially a square degree as for the same approximation accuracy compared to the case with a one-sided  spectrum.In [21], one finds further explanations why Krylov subspace methods can be inefficient for solving complex valued systems, represented in the above real forms. Several iterative solution methods, such as the QMR [22] have been developped and proven to be efficient for these types of problems. However, it is difficult to precondition complex valued matrices and unpreconditioned methods converge in general very slowly.Following [66], we consider here instead an approach based on rewriting the equation in the form (170).Instead of solving the full block matrix system we apply a Schur complement approach by the elimination of one component, which results in the following reduced system(171)Cx=f,
where (172)C=R+SR-1S,f=c+SR-1d.
As an introduction, assume first that R is symmetric and positive definite and S is symmetric and positive semidefinite. As a preconditioner to the matrix C in (171) we take R+S.For the generalized eigen value problem,(173)μ(R+S)z=(R+SR-1S)z,
it holds then(174)μ(I+H)y=(I+H2)y,
where μ=R1/2z and H=R-1/2SR-1/2.IfλRz=Sz,  z≠0, or, equivalently, Hy=λy,   y≠0, it follows from (174) that(175)μ=1+λ2(1+λ)2,
that is, (176)μ=11+(2λ/(1+λ2)).
Since, by assumption λ≥0, it follows(177)12≤μ≤1,
and the condition number satisfies the bound 𝒦((R+S)-1C)≤2.The correspondingly preconditioned conjugate gradient method to solve (174) converges therefore rapidlyThere exists an even more efficient form of the iteration matrix that also shows that we can weaken the assumptions made onR and S. Hence, consider(178)Rx-Sy=c,Sx+Ry-d
and assume that α is a real parameter such that R+αS is nonsingular. Such a parameter exists if ker(R)∩ker(S)=0.Multiplying the first equation by-α(R+αS)-1, the second by (R+αS)-1 and adding yields(179)(R+αS)-1(S-αR)x+y=(R+αS)-1(d-αc).
Now multiplying this equation by S, using Sy=Rx-c, and rewriting the equation properly, we find(180)r≡Rx-c+S(R+αS)-1((S-αR)x+αc-d)=0.When solving the system by iteration, such as by Chebyshev iterations,r will be the residual and we observe that r can be written in the form (see (179))(181)y=(R+αS)-1((S-αR)x+αc-d),r=Rx-Sy-c.In this form there is no need to compute the right hand side vectorf initially as if (168) is used and the vector y is found during the iteration process. This saves two solutions with the matrix R+αS. Since we need few iterations, such a saving can be important to decrease the total expense of the method.To solve (179) we use R+αS as a preconditioner. The resulting preconditioned matrix takes the form(182)Mα=(R+αS)-1[R+S(R+αS)-1(S-αR)]=(R+αS)-1[(R+αS)-αS](R+αS)-1R+(R+αS)-1S(R+αS)-1S=((R+αS)-1R)2+((R+αS)-1S)2.This form can also be used to derive eigenvalue estimates in more general cases than was done above. IfR is nonsingular, we find(183)Mα=(I+αR-1S)-2(I+R-1S)2.
Therefore, the preconditioned matrix is a rational function in the matrix R-1S. It follows that the eigenvalues μ of Eα satisfy(184)μ=1+λ2(1+αλ)2,
were λ is an eigenvalue of R-1S.We want to chooseα to minimize the spectral condition number(185)K(Mα)=μmax⁡μmin⁡.Theorem 11.
Assume thatR is s.p.d and S is s.p.s-d. Then, the extreme eigenvalues of the preconditioned matrix Mα, defined in (182) satisfy
(186)μmin⁡={11+α2if0≤α≤λ̂,1+λ̂2(1+αλ̂)2ifα≥λ̂,μmax⁡={1ifα̂≤α,1+λ̂2(1+αλ̂)2if0≤α≤α̂,
where λ̂ is the maximal eigenvalue of R-1S,
(187)R-1S≤λ̂I,α̂=λ̂1+1+λ̂2.
The spectral condition number is minimized when α=α̂, in which case
(188)μmin⁡=11+α̂2,μmax⁡=1,K(Mα)=1+α̂2=21+λ̂21+1+λ̂2.Proof.
The bounds of the extreme eigenvalues follow by elementary computations ofμ=(1+λ2)/(1+αλ)2,   0≤λ≤λ̂. Similarly, it is readily seen that μmax/μmax is minimized for some α in the interval α̂≤α≤λ̂, where μmax=1. Hence, it is minimized for α=argmaxα̂≤α(1+α2)-1, that is, for α=α̂.For applications, see [66]. An important application arises when one uses Padé type approximations, and related implicit Runge-Kutta methods (see [67]), to solve initial value problems.


### 4.7. Historical Remarks

Because incomplete factorization methods has had a strong influence on the development of preconditioning methods we give here some historical remarks.The idea of an incomplete factorization method goes back to early papers by Buleev [68], Varga [10], Oliphant [69], Dupont et al. [70], Dupont [71], and Woźnički [72], where it was presented for matrices of a type arising from difference approximations of elliptic problems. The first more general form (unmodified methods for pointwise matrices) was studied for M-matrices by Meijerink and van der Vorst [46]. For a review and general formalism for describing such methods, see Axelsson [18], Birkhoff et al. [63], Beauwens [12], and Il'in [60]. For a similar but more involved type of methods for difference matrices, which allowed for variable parameters from one iteration to the next, see Stone [73].A modified form of the method, where a certain row sum criterion was imposed, was studied by Gustafsson [47]. Actually, as is readily seen, the method of Dupont et al. [70] and as further discussed in Axelsson  [59], using a perturbation technique, can be seen as a modified version of the general incomplete factorization method when applied to the five-point elliptic difference matrices, assuming that no fill-in is accepted outside the sparsity structure of A itself and assuming a natural ordering of the grid points. The advantage of modified versions is that they can give condition numbers of the iteration matrices that are of an order of magnitude smaller than for the original matrix.The incomplete factorization method can be readily generalized to matrices partitioned in block matrix form. This was done first for matrices partitioned in block tridiagonal form in Axelsson et al. [49] and Concus et al. [50], the latter being based on earlier work by Underwood  [74]. A general form was presented in Axelsson  [67] and Beauwens and Ben Bouzid  [52], where existence of the method was proven for M-matrices.The existence, that is, the existence of nonzero pivot entries of pointwise incomplete factorization methods forM-matrices was first shown by Meijerink and van der Vorst  [46] and, for pointwise H-matrices, by Varga et al. [75]. The existence of incomplete factorization methods for M-matrices in block form was shown in Axelsson et al. [49] and Concus et al. [50] for block tridiagonal matrices; in Axelsson  [76] and Beauwens and Ben Bouzid [52], for general block matrices; and in Axelsson and Polman  [51] for relaxted versions of such methods.Kolotilina  [77] shows the existence of convergent splittings for block H-matrices, and Axelsson  [78] shows the existence of general incomplete factorizations for block H-matrices.


## 4.1. The Exact Inverse on Additive Form

Let thenA(0)=A and consider a matrix(87)A(k)=[A1(k)A12(k)A21(k)A2(k)],
in a sequence defined by(88)A(k+1)≡S2(k)=A2(k)-A21(k)A1(k)-1A21(k),k=0,1,…,k0,
where each A1(k) in nonsingular, being a block diagonal of a symmetric positive definite matrix. Hence, the following recursion holds(89)A(k)-1=[A1(k)-1000]+[-A1(k)-1A12(k)I2(k+1)]A(k+1)-1[A21(k)A1(k)-1,I2(k+1)],k=0,  1,…,k0.Here, I2(k+1) is the identity matrix on level k+1. Note that in this example the dimensions decrease with increasing level number and the final matrix (i.e., Schur complement) in the sequence is A(k0)=S2(k0+1). The above recursion can be rewritten in compact form(90)A-1=[A1(0)-1000]+UD-1L,
where the kth column of the block upper triangular matrix U equals [-A1(k)-1A12(k)I2(k+1)] and L=UT. Further,(91)D=[A1(k)0A(2)⋱0A(k0)].Hence, this is the (block matrix) Gaussian elimination method applied directly to form the inverse matrix. In this way, there is no need to first form the factorizationA=L̃DŨ and then A-1=Ũ-1D-1L̃-1. As is wellknown and readily seen, the columns of Ũ-1 and L̃-1 are formed directly with no additional computation, from those of Ũ and L̃, respectively. Note that Ũ-1 is upper (block) triangular and L̃-1 is lower (block) triangular).The matrixD contains the (block) pivot matrices which arise during the factorization. Permutations to increase the stability by finding dominating pivots can be done by replacing A(k+1) with P(k)TÃ(k+1)P(k) where Ã(k+1)=P(k)A(k+1)P(k)T is the permuted matrix on which the next elimination step takes place.An incomplete factorization method for approximate inverses can be defined by approximating each arising Schur complement matrix with some sparse matrixB̃(k+1) and possibly also approximating A1(k)-1 with some matrix B1(k), to yield the approximate inverse(92)B(k)=[B1(k)000]+V(k)B̃(k+1)V(k)T,
where V(k)=[-B1(k)A12(k)I2(k+1)].In forming the approximate Schur complement one can use a simpler matrixD1(k)than B1(k), often a diagonal matrix suffices. The intermediate Schur complement matrix S̃2(k)=A2(k)-A21(k)D1(k)A12(k) can be possibly further approximated by deleting certain off-diagonal entries to preserve sparsity. These entries can be compensated for by modifying the diagonal of S̃2(k) to form the final approximation B̃(k+1). Thereby, it can be important to make the approximate Schur complement exact on some particular vector or vector space. (We do not go into these aspects further here, see [43, 56] for details.)The eigenvectors for the smallest eigenvalues ofA provide efficient column vectors for the matrix V to reduce significantly the condition number of BA as compared to that of A. However, in general the eigenvectors are not known, and even if they are known it would be costly to apply the corresponding preconditioner as V would be a full matrix. A more viable choice is to let V be defined by the basis functions {φi(H)} of a coarse mesh (or coarsened matrix graph) so that(93)Im⁡V=span⁡{φi(H)}.V,   VT acts then, respectively, as prolongation and restriction operators and(94)V=[J12I2],
where J12 is the interpolation matrix from the coarse mesh (ΩH) to the fine mesh (Ωh), and we assume   ΩH⊂Ωh. Further, letting the matrices be variationally defined, as in a finite element method, we have(95)AH=VTAV,
where A is the finite element matrix on the fine mesh.Now, the eigenvectors for the smaller eigenvalues ofAH are normally accurate approximations of the corresponding eigenvectors for A. Furthermore, the eigenvectors of AH are members of ImV. Therefore, the matrix V in (94) acts nearly as well as the eigenvector matrix but, in addition, is sparse. Hence the approximate inverse takes the form(96)B=G+σVAH-1VT,
where σ=λmax(GA),  AH=VTAV.Here, the projection matrix(97)P=VAH-1VTA,
projects vectors on the subspace ImV, containing normally good approximations of the eigenvectors for the smallest eigenvalues of A, that is, those who may cause severe ill-conditioning. Clearly, the approximation is more accurate as closer ΩH is to Ωh. However, the cost of the action of P (mainly the coarse mesh solver for the action of AH-1) increases when ΩH expands. One can balance ΩH to Ωh in order to let the action of P involve the same order of computational complexity as an action of A, that is, O(h-2) for a sparse matrix A. Assuming that the cost of an action of AH-1 is O(H-2.5) in a 2D diffusion problem (e.g., using a modified incomplete factorization method as preconditioner for the conjugate gradient method), we find H=h4/5. The number of outer iterations with preconditioner B depends also on the choice of G. We refer the discussion of how G can be chosen properly to [56].As an example, for a model diffusion problem with constant coefficients on a regular mesh, say for the Laplacian operator on unit square, the eigenvectors for the Laplacian(-Δ) on Ωh with Dirichlet boundary conditions are(98)vk,l(h)=sin⁡kπsin⁡lπy,x,y∈Ωh,
where k,  l=1,2,…,h-1-1, for the eigenvalues(99)λk,l(h)=(2sin⁡kπh2)2+(2sin⁡lπh2)2.
For the coarse mesh, it holds(100)vk,l(H)=sin⁡kπxsin⁡lπy,x,y∈ΩH,
where k,  l,=1,2,…,H-1, and here Vvk,l(H) are good approximations (interpolants) of the eigenvectors vk,l(h) on Ωh for the smallest eigenvalues.An alternative choice of matrixV is to take eigenvectors from a nearby problem, normally defined by taking limit values of some problem parameter, see [56].Multigrid, algebraic multilevel and algebraic multigrid methods have been presented thoroughly in, for  example [29, 43, 57, 58]. Because of space limitations, they can not be presented here.


## 4.2. Symmetrization of Preconditioners; the SSOR and ADI Methods

As we have seen, the incomplete factorization methods require first a factorization step. There exists simpler preconditioning methods that require no factorization but have a form similar to the incomplete factorization methods. We will present two methods of this type. As an introduction, consider first an iterative method of the form(101)M(xl+1-xl)=b-Axl,l=0,1,…,
to solve Ax=b, where A and M are nonsingular. As we saw in Section 2, the asymptotic rate of convergence is determined by the spectral radius of the iteration matrix(102)B=I-M-1A.
For a method such as the SOR method (which also requires no factorization), with optimal overrelaxation parameter ω (assuming that A has property Aπ or A is s.p.d., see Section 2), the eigenvalues of the corresponding iteration matrix B are situated on a circle. No further acceleration is then possible.There is, however, a simple remedy to this, based on taking a step in the forward direction of the chosen ordering, followed by a backward step, that is, a step in the opposite order to the vector components. This method is said to have its origin in the early days of computers when programs were stored on tapes that had to be rewound before a new forward SOR step could begin. It was found that this otherwise useless computer time for the rewinding could be better used for a backward SOR sweep!As we will see, for symmetric and positive definite matrices the combined forward and backward sweeps correspond to a s.p.d. matrix which, contrary to the SOR method, has the advantage that it can be used as a preconditioning matrix in an iterative acceleration method. This method, called the SSOR method, will be defined later.For an early discussion of the SSOR method, used as a preconditioner, see [59]. For discussions about symmetrization of preconditioners, see [6, 60, 61]. More generally, if A is s.p.d, we consider the symmetrization of an iterative method in the form(103)xl+1=xl+M-1(b-Axl).
For the analysis only, we consider the transformed form of (103),(104)yl+1=(I-A1/2M-1A1/2)yl+b̃,
where(105)yl=A1/2xl,b̃=A1/2M-1b̃.
If M is unsymmetric, the iteration matrix I-A1/2M-1A1/2 is also unsymmetric. We will now consider a method using M and another preconditioner chosen so that the iteration matrix for the combined method becomes symmetric. We call this the symmetrization of the method.LetM1,  M2 be two such preconditioning matrices. Let(106)Bi=I-M̃i-1,M̃i=A-1/2MiA-1/2,
and consider the combined iteration matrix B2B1. As we will now see, it arises as an iteration matrix for the combined method(107)M1(xl+1/2-xl)=b-Axl,M2(xl+1-xl+1/2)=b-Axl+1/2,l=0,1,….
For the analysis only, we transform this to the form(108)yl+1/2-yl=b̃(1)-M̃1-1yl,yl+1-yl+1/2=b̃(2)-M̃2-1yl+1/2,
where(109)b̃(i)=A(1/2)Mi-1b.
This iteration takes the form(110)yl+1=b̃(2)+(I-M̃2-1)b̃(1)+(I-M̃1-1yl),
that is,(111)yl+1=b̃(2)+(I-M̃2-1)b̃(1)+(I-M̃2-1)(I-M̃1-1)yl,
or(112)yl+1=b̂+B2B1yl,l=0,1,…,
where(113)b̂=b̃(2)+(I-M̃2-1)b(1).
For the following we need a lemma.Lemma 2.
IfA,   B, and C are Hermitian positive definite and each pair of them commute, then ABC is Hermitian positive definite.Proof.
We have(ABC)*=CBA and use commutativity to find
(114)CBA=BCA+BAC=ABC.
Hence, ABC is Hermitian. Next, we show that the product of two s.p.d matrices that commute is positive definite. We have
(115)A-1/2ABA1/2ABA1/2=A1/2BA1/2,
which is Hermitian positive definite. Hence, by similarity, the eigenvalues of AB are positive and, since
(116)(AB)*=AB,AB is Hermitian positive definite. In the same way, (AB)C is Hermitian positive definite.Lemma 3.
LetA be s.p.d. and assume either of the following additional conditions:(a)
M2*=M1.(b)
M1,  M2 are s.p.d. ρ(A1/2Mi-1A1/2)<1,   i=1,2, and the pair of matrices M1, M2, commutes.Then, the combined iteration method (107) converges if and only if M1+M2-A is s.p.d.Proof.
It is readily seen thaty=b̂+B2B1y (i.e., the iteration method is consistent with Ax=b), where y=A1/2x and x=A-1b. Hence,
(117)y-yl+1=B2B1(y-yl),
and the iteration method (112), and hence (107) converges for any initial vector if and only if ρ(B2B1)<1, where ρ(·) denotes the spectral radius. But
(118)B2B1=I-M̃1-1-M̃2-1+M̃1-1M̃2-1=I-A1/2M1-1(M1+M2-A)M2-1A1/2.
It is readily seen that under either of the given conditions (a) or (b),
(119)M1-1(M1+M2-A)M2-1=M1-1+M2-1-M1-1AM2-1
is symmetric. Further, it is positive definite if and only if M1+M2-A is positive definite. Hence, I-B2B1 is s.p.d. Further, B2B1=(I-M̃2-1)(I-M̃1-1) is symmetric, and a similarity transformation shows that B2B1 is similar to (I-M̃2-1)1/2(I-M̃1-1)(I-M̃2-1)1/2, which is a congruence transformation of I-M̃1-1, whose eigenvalues are positive. Hence, B2B1 has positive eigenvalues, so the eigenvalues of B2B1 are contained in the interval (0,1) and, in particular, ρ(B2B1)<1.The proof of Lemma3 shows that B2B1 is symmetric, so the combined iteration method is a symmetrized version of either of the simple methods.Let us now consider a special class of symmetrized methods. We letA be split as A=D+L+U, where we assume that D is s.p.d., and let(120)V=(1-1ω)D+L,H=(1-1ω)D+U,D̂=(2/ω-1)D, where ω is a parameter, 0<ω<2. (Here, L and U are not necessarily the lower and upper triangular parts of A.) Note that(121)D̂+V+H=A,
so this is also a splitting of A. As an example of a combined, or symmetrized, iteration method, we consider the preconditioning matrix(122)C=(D̂+V)D̂-1(D̂+H),
and show that this leads to a convergent iteration method(123)C(xl+1-xl)=b-Axl,l=0,1,….
This corresponds to choosing M1=D̂-1/2(D̂+H) and M2=(D̂+V)D̂-1/2, and it can be seen that the conditions of Lemma 3 hold if the conditions in the next theorem hold.Theorem 8.
LetA=D+L+U, where D is s.p.d. Let V,H,D̂ be defined by (120), and assume that either (a) or (b) holds, where(a)
U=L*(b)
L,U are s.p.d. and each pair of matrices L, U, D commute. Then the eigenvalues λ of the matrix C-1A, where C is defined in (122), are contained in the interval 0<λ≤1.Proof.
This can be shown either by verifying the conditions in Lemma3 or more directly as follows. As in the proof of Lemma 3, it follows that C is s.p.d. Hence, the eigenvalues of C-1A are positive. Further,
(124)C=D̂+V+H+VD̂-1H,
so, by (121)
(125)C=A+VD̂-1H.
Under either condition (a) or (b), C=VD̂-1H is symmetric and positive semidefinite.
This shows thatx*Cx≥x*Ax for all  x, so the eigenvalues of C-1A are bounded above by 1.We will now show that the matrixC can also efficiently be used as a preconditioning matrix, which for a proper value of the parameter ω, and under an additional condition, can even reduce the order of magnitude of the condition number. In this respect, note that when C is used as a preconditioning matrix for the Chebyshev iterative method, it is not necessary to have C scaled so that λ(C-1A)≤1, because it is suffices then that 0<m≤λ(C-1A)≤M, for some numbers m, M. Hence, the factor 2/ω-1 in D̂-1 can be neglected.Theorem 9.
LetA=D+L+U be a splitting of A, where A and D are s.p.d. and either (a) U=L* or (b) L,  U are s.p.d. and each pair of D,  L,  U commute. Then, the eigenvalues of matrix C-1A, where
(126)C=(1ωD+L)D̂-1(1ωD+U)
and 0<ω<2,  D̂=(2/ω-1)D, are contained in the interval
(127)[(2-ω){1+ω(1/ω-1/2)2δ-1+ωγ},1],
where
(128)δ=min⁡x≠0xTAxxTDxγ=max⁡x≠0xT(LD-1U-1/4D)xxTAx.
Further, if there exists a vector for whichxT(L+U)xT(L+U)x≤0, then γ≥-1/4, and if(129)ρ(L̃Ũ)≤14,
then γ≤0, and if
(130)ρ(L̃Ũ)≤14+O(δ),thenγ≤O(1),δ⟶0.
Here, L̃=D-1/2LD-1/2.Proof.
It is readily seen that(131)C=12-ω(1ωD+L)(1ωD)-1(1ωD+U)=12-ω[A+(1ω-1)D+ωLD-1U]=12-ω[A+ω(1ω-12)2D+ω(LD-1U-14D)].
This shows the lower bound in (127); the upper bound follows by Theorem 8. By choosing a vector for which xT(L+U)x≤0, it follows that
(132)xT(LD-1U-(1/4)D)xxTAx≥xT(LD-1x-(1/4)Dx)xTDx≥-14,
which shows γ≥-1/4. The remainder of the theorem is immediate.


### 4.2.1. The Condition Number

Theorem9 shows that the optimal value of ω to minimize the upper bound of the condition number of C-1A is the value that minimizes the real-valued  function (133)f(ω)=1+ω(1/ω-1/2)2δ-1+ωγ2-ω.
It is readily seen (see Axelsson and Barker, 1984 [30]), that f(ω) is minimized for(134)ω*=21+2(1/2+γ)δ,min⁡ωf(ω)=f(ω*)=(12+γ)δ-1+12.
In general, δ is not known, but we may know that δ=O(h2), for some problem parameter, h→0 (such as for the step length in second-order elliptic problems). Then, if γ=O(1),   h→0, we let ω=2/(1+ξh) for some ξ>0, in which case(135)f(ω)=O(h-1)=O(δ-1),h⟶0.This means thatC-1A has an order of magnitude smaller condition number than A itself, which latter is O(δ-1).We consider now two applications of Theorem9.


### 4.2.2. The SSOR Method

In the first case,L is the lower triangular part of A or the lower block triangular part, if A is partitioned in block matrix form and U=L*. Then,(136)C=12-ω(1ωD+L)(1ωD)-1(1ωD+L*),
is a symmetrized version of the SOR method and is called the SSOR (symmetric successive overrelaxation) method.As an example, for an elliptic differential equation of second order it can be seen that the conditionρ(L̃L̃T)≤1/4 holds for problems with Dirichlet boundary conditions and constant coefficients. For extensions of this, see Axelsson and Barker [30]. For the model difference equation on a square domain with side π, we have(137)δ=2(sin⁡h2)2,γ≤0,
and we find(138)ω*=21+2sin⁡h/2~21+h,f(ω*)=12δ+12~h-1+12,h⟶0.


## 4.2.1. The Condition Number

Theorem9 shows that the optimal value of ω to minimize the upper bound of the condition number of C-1A is the value that minimizes the real-valued  function (133)f(ω)=1+ω(1/ω-1/2)2δ-1+ωγ2-ω.
It is readily seen (see Axelsson and Barker, 1984 [30]), that f(ω) is minimized for(134)ω*=21+2(1/2+γ)δ,min⁡ωf(ω)=f(ω*)=(12+γ)δ-1+12.
In general, δ is not known, but we may know that δ=O(h2), for some problem parameter, h→0 (such as for the step length in second-order elliptic problems). Then, if γ=O(1),   h→0, we let ω=2/(1+ξh) for some ξ>0, in which case(135)f(ω)=O(h-1)=O(δ-1),h⟶0.This means thatC-1A has an order of magnitude smaller condition number than A itself, which latter is O(δ-1).We consider now two applications of Theorem9.


## 4.2.2. The SSOR Method

In the first case,L is the lower triangular part of A or the lower block triangular part, if A is partitioned in block matrix form and U=L*. Then,(136)C=12-ω(1ωD+L)(1ωD)-1(1ωD+L*),
is a symmetrized version of the SOR method and is called the SSOR (symmetric successive overrelaxation) method.As an example, for an elliptic differential equation of second order it can be seen that the conditionρ(L̃L̃T)≤1/4 holds for problems with Dirichlet boundary conditions and constant coefficients. For extensions of this, see Axelsson and Barker [30]. For the model difference equation on a square domain with side π, we have(137)δ=2(sin⁡h2)2,γ≤0,
and we find(138)ω*=21+2sin⁡h/2~21+h,f(ω*)=12δ+12~h-1+12,h⟶0.


## 4.3. The ADI Method

In the second case of methods of (101), we let L denote the off-diagonal part of the difference operator working in the x-direction and U off-diagonal part of the difference operator in the y-direction. D is its diagonal part. Then, the matrix(139)Ĉ=(1ωD+L)(1ωD)-1(1ωD+U),
is called an alternating direction preconditioning matrix and the corresponding iteration method is called the ADI (alternating direction iteration) method. In this method, we solve alternately one-dimensional difference equations in x-and y-directions. Much has been written on the ADI-method which was originally presented in Peaceman and Rachford [62]; see Varga [10], for an early influential presentation and Birkhoff et al. [63] and Wachspress [64], for instance.As we will see, for the model difference equations we get the same optimal value ofω as in (138). The condition γ=O(1) may be less restrictive, for the ADI-method, but the condition of commutativity is much more restrictive, as the following lemma shows.Lemma 4.
LetA,  B be two Hermitian matrices of order n. Then AB=BA if and only if A and B have a common set of orthonormal eigenvectors.Proof.
If such a common set of eigenvectors{vi} exists, then Avi=σivi,  Bvi=τivi and
(140)ABvi=σiτivii=1,2,…,n.
Since the eigenvector space of an Hermitian matrix is complete, we therefore have
(141)ABx=BAx,∀x∈Cn,
which shows that AB=BA. Conversely, suppose that AB=BA. As A is Hermitian, take U to be a unitary matrix that diagonalizes A, that is
(142)Ã=UAU*=[γ1I10γ2I2⋱0γrIr],
where γ1<γ2<⋯<γr are the distinct eigenvalues of A and Ij is the identity matrix of order nj, the multiplicity of γj. (Here A is posssibly permuted accordingly.) Let B̃=UBU* and partition B̃ corresponding to the partitioning of A, that is,
(143)B̃=[B11B12⋯B1r⋮Br1Br2⋯Brr].
Since AB=BA, we have
(144)ÃB̃=UABU*=UBAU*=B̃Ã.
Carying out the block multiplication ÃB̃=B̃Ã, we find that this, in turn, implies Bij=0,  i≠j, since γi≠γj,  i≠j. Simply stated, a (block) matrix commutes with a (block) diagonal matrix if and only if it is itself (block) diagonal. Hence, B̃ is block diagonal and each Hermitian submatrix Bi,i has ni orthonormal eigenvectors that are also eigenvectors of the submatrix γiIi of Ã. Since ∑i=1rni=n and all eigenvectors are orthonormal, A and B must have the same set of eigenvectors.For the second-order elliptic difference equation in two space dimensions, it turns out that the commutativity ofL and U essentially corresponds to the property that the original problem is separable, that is, that solutions of ℒu=f can be written in the form u=φ(x)ψ(y). This means that the coefficients a(x,y) and b(x,y) in the differential operator ∂/∂x[a(x,y)∂u/∂x]+∂/∂y[b(x,y)∂u/∂x]+c(x,y)u must satisfy a(x,y)=a(x),b(x,y)=b(y), and c(x,y)=c, a constant. Hence, if a(x,y)=b(x,y), then a(x,y)=b(x,y)=a, a constant. Furthermore, the convex closure of the meshpoints must be a rectangle with sides parallel to the coordinate axes (Varga, [10] 1962). If A=A1+A2, the ADI-method can be written in the form(145)(I+τ1A1)xl+1/2=(I-τ1A2)xl+τ1b,(I+τ2A2)xl+1=(I-τ2A1)xl+1/2+τ2b,l=0,1,….This is thePeaceman-Rachford[62]iteration method. The iteration matrix M is similar to  (146)(I-τ2A1)(I+τ1A1)-1(I-τ1A2)(I+τ2A2)-1.
When A1,  A2 are Hermitian positive definite, their eigenvalues λi(1),  λi(2) are positive, and(147)‖(I-τ2A1)(I+τ1A1)-1‖2=ρ((I-τ2A1)(I+τ1A1)-1)=max⁡i|1-τ2λi(1)1+τ1λi(1)|.
Thus,(148)ρ(M)=ρ((I-τ2A1)(I+τ1A1)-1(I-τ1A2)-1(I+τ2A2)-1)≤‖(I-τ2A1)(I+τ1A1)-1(I-τ1A2)(I+τ2A2)-1‖2≤‖(I-τ2A1)(I+τ1A1)-1‖2×‖(I-τ1A2)(I+τ2A2)-1‖2=μ(τ1,τ2),ρ(M)≤μ(τ1,τ2)=max⁡i|1-τ2λi(1)1+τ1λi(1)|max⁡i|1-τ1λi(2)1+τ2λi(2)|.Note that forτ1=τ2=τ>0,  μ(τ1,τ2)=μ(τ,τ)<1, so we have ρ(M)<1, that is convergence for any τ>0. This holds even if A1,A2 do not commute. Note also that when A1 and A2 commute, we have(149)ρ(M)=μ(τ1,τ2).Let us continue the analyses for the general case whereA1,  A2 do not necessarily commute. We want to compute the optimal values of τ1 and τ2 such that μ(τ1,τ2) is minimized. For simplicity, we assume that α, β are the same lower and upper bounds of the eigenvalues of A1 and A2, that is, 0<α≤λi(j)≤β, j=1,2. We have(150)μ(τ1,τ2)≤max⁡{|1-τ2α1+τ1α|,|1-τ2β1+τ1β|}×max⁡{|1-τ1α1+τ2α|,|1-τ1β1+τ2β|}.
We want to choose τ1,  τ2>0 such that this bound is as small as possible. Note, then, that for such values of τ1,  τ2 we must have 1-τiα>0 and 1-τiβ<0. Next note that each factor in the bound (150) is minimized when(151)1-τ2α1+τ1α=τ2β-1τ1β+1,1-τ1α1+τ2α=τ1β-1τ2β+1,
respectively, that is, when(152)τ1τ2-α+β2αβ(τ1-τ2)-1αβ=0,τ1τ2+α+β2αβ(τ1-τ2)-1αβ=0,
respectively. Thus, both factors are simultaneously minimized when τ1=τ2, and then τ1=τ2=1/αβ.Theorem 10.
LetA=A1+A2, where A1,  A2, are s.p.d., and consider the Peaceman-Rachford ADI method (145) to solve Ax=b with τ1=τ2=1/αβ. The spectral radius of the corresponding iteration matrix M satisfies
(153)ρ(M)≤min⁡τ1,τ2μ(τ1,τ2)≤(1-α/β1+α/β)2~1-4αβ,
if α/β→0.Proof.
Forτ1=τ2=1/αβ, we have
(154)μ(τ1,τ2)=(1-α/β1+α/β)2.Remark 2.
For a model difference equation for a second-order elliptic differential equation problem on a square with sideπ, we have with stepsize h,
(155)α=(sin⁡(h/2)h/2)2~1,β=(cos⁡(h/2)h/2)2~4h2,h⟶0.
Then,
(156)μ(τ1,τ2)=(1-tan(h/2)1+tan(h/2))2=1-sin⁡(h)1+sin⁡(h)~1-2h,h⟶0.
Note that this is just the convergence factor we get for the SOR method with an optimal overrelaxation parameter.
Sinceρ(M)≤μ(τ1,τ2), this means that the ADI method with parameters (chosen as above) converges at least as fast as the SOR method. Note, however, that in the ADI-method we must solve two systems of equations with tridiagonal coefficient matrices (I-τAi) on each step, while the pointwise SOR method requires no solution of such systems.


## 4.4. The Commutative Case

Assume now thatA1,  A2 commute. Then, as we have seen, M is symmetric and has real eigenvalues, and we can apply the Chebyshev acceleration method. The eigenvalues of the corresponding preconditioned matrix C̃ are related to the eigenvalues of M by(157)λ(C̃)=1-λ(M).
Since -ρ(M)≤λ(M)≤ρ(M), where(158)ρ(M)=(1-α/β1+α/β)2~1-4αβ,
we have(159)4αβ~1-ρ(M)≤λ(C̃)=1+ρ(M)~2-4αβ~2,αβ⟶0.
The asymptotic rate of convergence of the Chebyshev accelerated method therefore is(160)21-ρ(M)1+ρ(M)~22(αβ)1/4,α/β⟶0.
For the model difference equation, we have the asymptotic rate of convergence(161)~2h1/2,h⟶0.


## 4.5. The Cyclically Repeated ADI Method

The real power of the ADI method is brought forth when we use a sequence of parametrsτl. Assume that A1,  A2 commute, then we choose the parameters τl cyclically. With a cycle of q parameters and with the assumption(162)0<α≤λi(j)≤β,j=1,2,
we get the iteration matrix(163)M(q)=∏p=1q(I+τpA2)-1(I-τpA1)(I+τpA1)-1(I-τpA2).
The eigenvalues of M(q) are(164)∏p=1q‍1-τpλi(1)1+τpλi(1)⋅1-τpλi(2)1+τpλi(2).
In the same way as above, ρ(M(q)) is minimized when(165)d(α,β,q)=max⁡α≤x≤β∏p=1q‍|1-τpx1+τpx|.


## 4.6. A Preconditioning Method for Complex Valued Matrices

Complex valued systems of equations arise in many applications. A commonly occuring case is the solution of a matrix polynomial equation(166)Qm(A)x=b,
where A is a real square matrix and Qm is a polynomial of degree m that has no zeroes at the eigenvalues of A. Here Qm can be factored in the product of second degree, and possibly some factors of first degree polynomials with real coefficients.The second degree polynomials can be factored in products of first degree polynomials with complex coefficients.Consider then a linear system(167)Az=b.in the form(168)(R+iS)(x+iy)=c+id,
where ℛ,   S are real matrices of order n and x,y,c,d∈ℛn.The system can be solved in complex arithmetic. However, complex arithmetic leads to heavier computational complexity and it is in general difficult to precondition complex valued matrices, as the eigenvalues of the given matrix or the preconditioned matrix can be spread in the whole complex plane and the iterative solution method may then converge too slowly.One can alternatively apply a preconditioned conjugate gradient method to the Hermitian positive definite normal matrix systemAHAu=AHb for which the eigenvalues are real. At any rate, this involves complex arithmetic that costs typically three to four times as much as corresponding real arithmetic.Complex arithmetic can be avoided by rewriting (168) in real valued form, such as(169)A(1)[xy]=[R-SSR][xy]=[cd],or(170)A(2)[x-y]=[RSS-R][x-y]=[cd].The block matrices are here real but, in general, nonsymmetric and/or indefinite. For the solution, one can use a generalized conjugate gradient method such as GMRES [40] or GCG [65].ForA(1) it holds that any eigenvalue λ appears in complex conjugate pairs λ,  λ¯. For A(2), which is real symmetric, for any eigenvalue λ≠0,-λ is also an eigenvalue. Thus, the spectrum σ(A(1)) is symmetric with respect to the real axes and the spectrum σ(A(2)) is symmetric with respect to the imaginary axes, that is, in both cases the spectrum embraces the origin. From best polynomial approximation properties it is known that such point distributions leads to polynomials of essentially a square degree as for the same approximation accuracy compared to the case with a one-sided  spectrum.In [21], one finds further explanations why Krylov subspace methods can be inefficient for solving complex valued systems, represented in the above real forms. Several iterative solution methods, such as the QMR [22] have been developped and proven to be efficient for these types of problems. However, it is difficult to precondition complex valued matrices and unpreconditioned methods converge in general very slowly.Following [66], we consider here instead an approach based on rewriting the equation in the form (170).Instead of solving the full block matrix system we apply a Schur complement approach by the elimination of one component, which results in the following reduced system(171)Cx=f,
where (172)C=R+SR-1S,f=c+SR-1d.
As an introduction, assume first that R is symmetric and positive definite and S is symmetric and positive semidefinite. As a preconditioner to the matrix C in (171) we take R+S.For the generalized eigen value problem,(173)μ(R+S)z=(R+SR-1S)z,
it holds then(174)μ(I+H)y=(I+H2)y,
where μ=R1/2z and H=R-1/2SR-1/2.IfλRz=Sz,  z≠0, or, equivalently, Hy=λy,   y≠0, it follows from (174) that(175)μ=1+λ2(1+λ)2,
that is, (176)μ=11+(2λ/(1+λ2)).
Since, by assumption λ≥0, it follows(177)12≤μ≤1,
and the condition number satisfies the bound 𝒦((R+S)-1C)≤2.The correspondingly preconditioned conjugate gradient method to solve (174) converges therefore rapidlyThere exists an even more efficient form of the iteration matrix that also shows that we can weaken the assumptions made onR and S. Hence, consider(178)Rx-Sy=c,Sx+Ry-d
and assume that α is a real parameter such that R+αS is nonsingular. Such a parameter exists if ker(R)∩ker(S)=0.Multiplying the first equation by-α(R+αS)-1, the second by (R+αS)-1 and adding yields(179)(R+αS)-1(S-αR)x+y=(R+αS)-1(d-αc).
Now multiplying this equation by S, using Sy=Rx-c, and rewriting the equation properly, we find(180)r≡Rx-c+S(R+αS)-1((S-αR)x+αc-d)=0.When solving the system by iteration, such as by Chebyshev iterations,r will be the residual and we observe that r can be written in the form (see (179))(181)y=(R+αS)-1((S-αR)x+αc-d),r=Rx-Sy-c.In this form there is no need to compute the right hand side vectorf initially as if (168) is used and the vector y is found during the iteration process. This saves two solutions with the matrix R+αS. Since we need few iterations, such a saving can be important to decrease the total expense of the method.To solve (179) we use R+αS as a preconditioner. The resulting preconditioned matrix takes the form(182)Mα=(R+αS)-1[R+S(R+αS)-1(S-αR)]=(R+αS)-1[(R+αS)-αS](R+αS)-1R+(R+αS)-1S(R+αS)-1S=((R+αS)-1R)2+((R+αS)-1S)2.This form can also be used to derive eigenvalue estimates in more general cases than was done above. IfR is nonsingular, we find(183)Mα=(I+αR-1S)-2(I+R-1S)2.
Therefore, the preconditioned matrix is a rational function in the matrix R-1S. It follows that the eigenvalues μ of Eα satisfy(184)μ=1+λ2(1+αλ)2,
were λ is an eigenvalue of R-1S.We want to chooseα to minimize the spectral condition number(185)K(Mα)=μmax⁡μmin⁡.Theorem 11.
Assume thatR is s.p.d and S is s.p.s-d. Then, the extreme eigenvalues of the preconditioned matrix Mα, defined in (182) satisfy
(186)μmin⁡={11+α2if0≤α≤λ̂,1+λ̂2(1+αλ̂)2ifα≥λ̂,μmax⁡={1ifα̂≤α,1+λ̂2(1+αλ̂)2if0≤α≤α̂,
where λ̂ is the maximal eigenvalue of R-1S,
(187)R-1S≤λ̂I,α̂=λ̂1+1+λ̂2.
The spectral condition number is minimized when α=α̂, in which case
(188)μmin⁡=11+α̂2,μmax⁡=1,K(Mα)=1+α̂2=21+λ̂21+1+λ̂2.Proof.
The bounds of the extreme eigenvalues follow by elementary computations ofμ=(1+λ2)/(1+αλ)2,   0≤λ≤λ̂. Similarly, it is readily seen that μmax/μmax is minimized for some α in the interval α̂≤α≤λ̂, where μmax=1. Hence, it is minimized for α=argmaxα̂≤α(1+α2)-1, that is, for α=α̂.For applications, see [66]. An important application arises when one uses Padé type approximations, and related implicit Runge-Kutta methods (see [67]), to solve initial value problems.


## 4.7. Historical Remarks

Because incomplete factorization methods has had a strong influence on the development of preconditioning methods we give here some historical remarks.The idea of an incomplete factorization method goes back to early papers by Buleev [68], Varga [10], Oliphant [69], Dupont et al. [70], Dupont [71], and Woźnički [72], where it was presented for matrices of a type arising from difference approximations of elliptic problems. The first more general form (unmodified methods for pointwise matrices) was studied for M-matrices by Meijerink and van der Vorst [46]. For a review and general formalism for describing such methods, see Axelsson [18], Birkhoff et al. [63], Beauwens [12], and Il'in [60]. For a similar but more involved type of methods for difference matrices, which allowed for variable parameters from one iteration to the next, see Stone [73].A modified form of the method, where a certain row sum criterion was imposed, was studied by Gustafsson [47]. Actually, as is readily seen, the method of Dupont et al. [70] and as further discussed in Axelsson  [59], using a perturbation technique, can be seen as a modified version of the general incomplete factorization method when applied to the five-point elliptic difference matrices, assuming that no fill-in is accepted outside the sparsity structure of A itself and assuming a natural ordering of the grid points. The advantage of modified versions is that they can give condition numbers of the iteration matrices that are of an order of magnitude smaller than for the original matrix.The incomplete factorization method can be readily generalized to matrices partitioned in block matrix form. This was done first for matrices partitioned in block tridiagonal form in Axelsson et al. [49] and Concus et al. [50], the latter being based on earlier work by Underwood  [74]. A general form was presented in Axelsson  [67] and Beauwens and Ben Bouzid  [52], where existence of the method was proven for M-matrices.The existence, that is, the existence of nonzero pivot entries of pointwise incomplete factorization methods forM-matrices was first shown by Meijerink and van der Vorst  [46] and, for pointwise H-matrices, by Varga et al. [75]. The existence of incomplete factorization methods for M-matrices in block form was shown in Axelsson et al. [49] and Concus et al. [50] for block tridiagonal matrices; in Axelsson  [76] and Beauwens and Ben Bouzid [52], for general block matrices; and in Axelsson and Polman  [51] for relaxted versions of such methods.Kolotilina  [77] shows the existence of convergent splittings for block H-matrices, and Axelsson  [78] shows the existence of general incomplete factorizations for block H-matrices.


## 5. Approximate Inverses Methods

In many applications, it is of interest to compute approximations of the inverse(A-1) of a given matrix A, such that these approximations can be readily used in various iterative methods.LetG denote an approximation of A-1.Following [6], first we present an example of an explicit and an implicit method, which is followed by a general framework for computing approximate inverses. At the end, we present an efficient way to construct symmetric and positive definite approximate inverses.An approximate inverse to a given operator may be constructed in several ways. The simplest way is to use a Neumann expansion, that is letD-1A=I-E, where D is the diagonal of A, for instance.Assuming that∥E∥<1, then the expansion(189)A-1=(I-E)-1D-1=(I+E+E2+⋯)D-1is convergent and any truncated part of this series provides an approximate inverse. However, this will normally give poore approximations. As we will see, more accurate approximate inverses can be constructed as best, possibly weighted, Frobenius norm approximations.In many applications the matrixA is sparse, but the exact inverse will be just a full matrix. A natural condition on G then arises: we can impose that G has some a priori chosen sparsity pattern (the same as A or different) which will make the calculations with G easy and cheap, and also will provide a sufficient accuracy.LetA have order n and S̲={(i,j),1≤i≤n;1≤i≤j≤n}. Any proper subset S of S̲ will be referred to as a sparsity pattern S⊂S̲.SL denotes the corresponding sparsity pattern for the lower triangular matrix and SL̃ denotes the corresponding sparsity pattern for the strictly lower triangular matrix.For simplicity, we use the same notationS for matrices having sparsity pattern S. Thus, A∈S if aij≠0⇔(i,j)∈S.


### 5.1. Explicit Methods

In these methods, an approximation of the inverseA-1 of a given nonsingular matrix A is computed explicitly, that is, without solving a linear globally coupled system of equations.LetS be a sparsity pattern. We want to compute G∈S, such that(190)(GA)ij=δij,(i,j)∈S,that is(191)∑k:(i,k)∈Sgikakj=δij,(i,j)∈S.
Some observations can be made from (191):(i)
the elements in each row ofG can be computed independently;(ii)
even ifA is symmetric, G in not necessarily symmetric, because gi,j,  j≠i,  and gj,i are, in general, not equal.


### 5.2. Implicit Methods

These methods require thatA is factored first. In practice, they are used mainly for band or “envelope" matrices. The algorithm was presented in [79]. It is based on an idea in [80]; see also [81].Suppose thatA=LD-1U is a triangular matrix factorization of A. If A is a band matrix then L and U are also band matrices.Let(192)L=I-L̃,U=I-Ũ,where L̃ and Ũ are strictly lower and upper triangular matrices correspondingly.The following lemma holds.Lemma 5.
Using the above notations it holds that(i)
A-1=DL-1+ŨA-1,(ii)
A-1=U-1D+A-1L̃.Proof.
Consider the following(193)A=LD-1U⟹A-1=U-1DL-1⟹(I-Ũ)A-1=DL-1⟹A-1=DL-1+ŨA-1.
Also,
(194)A-1(I-L̃)=U-1D⟹A-1=U-1D+A-1L̃.SinceDL-1 is lower triangular and Ũ is upper triangular, using (i) we can compute entries in the upper triangular part of A-1 with no need to use entries of L-1. Similarly, using (ii) we can compute entries of the lower triangular part A-1 without computing U-1.Suppose now thatA is a block banded matrix with a semibandwidth p, and we want to form A-1 also as block banded with a semibandwidth q:q≥p. The identities (i) and (ii) can be used then for the computation of the upper and lower parts of A-1.Remark 3.
(i) The algorithm involves onlymatrix×matrix operations.
(ii) There is no need to compute any entries outside the bands.
(iii) IfA is symmetric then it suffices executing only (i) or (ii).
(iv) It can be seen that(A-1)nn=Dnn-1.
There are two drawbacks with the above algorithm. It requires first the factorizationA=LD-1U and even if A is s.p.d, the band matrix part of A-1, which is computed, need not be s.p.d. The next example illustrates this.Example 2.
Consider an s.p.d. matrix(195)G=[1-21-25-31-34],Gband=[1-20-25-30-34],
is indefinite.


### 5.3. A General Framework for Computing Approximate Inverses

It turns out that both the explicit and implicit method can be characterized as methods to compute best approximations ofA-1 of all matrices having a given sparsity pattern, in some norm. The basic idea is due to Kolotilina and Yeremin  [38, 79], see [6]. Recall that the trace function is defined by tr(A)=∑i=1naii, which also equals ∑i=1nλi(A). Let a sparsity pattern S be given. Consider the functional(196)FW(G)≡‖I-GA‖W2=tr⁡((I-GA)W(I-GA)T),where the weight matrix W is s.p.d. If W≡I then ∥I-GA∥I is the Frobenius norm of I-GA.ClearlyFW(G)≥0. If G=A-1 then FW(G)=0. We want to compute the entries of G in order to minimize FW(G), that is, to find Ĝ∈S, such that(197)‖I-ĜA‖W≤‖I-GA‖W,∀G∈S.The following properties of the trace function will be used(198)tr⁡A=tr⁡AT,tr⁡(A+B)=tr⁡A+tr⁡B.Then,(199)FW(G)=tr⁡(I-GA)W(I-GA)T=tr⁡(W-GAW-W(GA)T+GAW(GA)T)=tr⁡W-tr⁡GAW-tr⁡(GAW)T+tr⁡GAWATGT.Further, as we are interested in minimizingFW with respect to G∈S, we consider the entries gi,j as variables. The necessary condition for a minimizing point are then(200)∂FW(G)∂gij=0,(i,j)∈S.
From (199) and (200), we get(201)-2(WAT)ij+2(GAWAT)ij=0,or(202)(GAWAT)ij=(WAT)ij,(i,j)∈S.
Depending on the particular matrix A and the choice of S and W,  (202) may or may not have a solution. We give some examples where a solution exists.Example 3.
LetA be s.p.d. Choose W=A-1 which is also s.p.d. Then, (202) implies
(203)(GA)ij=δij,(i,j)∈S,
which is the formula for the previously presented explicit method which, hence, is a special case of the more general framework for computing approximate inverses using weighted Frobenius norm.Example 4.
LetW=(ATA)-1. Then (202) implies
(204)(G)ij=(A-1)ij,(i,j)∈S,
which is the relation for the previously presented implicit method. In this case the entries of G are the corresponding entries of the exact inverse.Example 5.
LetW=I. Then,
(205)FW(G)=n-tr⁡(GA),(GAAT)ij=(AT)ij,(i,j)∈S.
This method is also explicit.
We can expect that such methods will be accurate only if all elements ofA which are not used in the computations are zero or are relatively small. In some cases the quality of the computed approximation G to A-1 can be significantly improved using diagonal compensation of the entries of A which are outside S. The best approximation G to A-1 in a (weighted) Frobenius norm is in general not symmetric and, as we have seen, not always positive definite. For this reason, the next, alternate method, is considered.


### 5.4. Constructing a Symmetric and Positive Definite Approximate Inverse

For some methods (as in the preconditioned Chebyshev and the conjugate gradient iteration methods) it is of importance to use s.p.d. preconditioners. As we have seen, the methods described till now do not guarantee thatG will be such a matrix.In order to compute an s.p.d. approximate inverse of an s.p.d. matrix, we can proceed as follows. It will be shown that this approximation gives a best approximation to minimize theK-condition number of the correspondingly preconditioned matrix.A Symmetric and Positive Definite Approximate Factorized Inverse
Seek an approximate inverse in the formG=LLT, where L∈SL,
(206)SL={(i,j)∈S,i≥j},
and let
(207)SL̃={(i,j)∈SL,i>j},
that is, denote by SL and SL̃ the lower and strictly lower triangular part of the sparsity set S.Theorem 12.
LetA be s.p.d. and consider matrices L with sparsity pattern Sl. Let the matrix L̂ be computed by the following steps.(i)
Compute firstL̃ such that
(208)(L̃A)ij=Aij,(i,j)∈SL̃,
and L̃ij=0,  (i,j)∈SL̃c (the complement set).(ii)
LetL̂=D(I-L̃), where
(209)D=diag⁡(d1,d2,…,dn),di=1[(I-L̃)A(I-L̃T)]i,i1/2.
Then, L̂∈SL and minimizes the K-condition number of L̂AL̂T, that is,
(210)((1/n)tr⁡(L̂AL̂T))ndet⁡(L̂AL̂T)=inf⁡L∈SL((1/n)tr⁡(LALT))ndet⁡(LALT).Proof.
LetD=diag(d1,…,dn) denote the diagonal part of a matrix X∈SL and let X̃=I-D-1X, that is, X̃∈SL̃. Then,
(211)((1/n)tr⁡(XAXT))ndet⁡(XAXT)=((1/n)∑i(XAXT)ii)n(det⁡(X))2det⁡(A)=((1/n)∑i[D(I-X̃)A(I-X̃T)D]i,i)n(det⁡(X))2det⁡(A)=((1/n)∑idi2[(I-X̃)A(I-X̃T)]i,i)n(Πidi2)det⁡(A)=((1/n)∑iαi2)nΠiαi2⋅Πi[(I-X̃)A(I-X̃T)]i,idet⁡(A),
where αi2=di2[(I-X̃)A(I-X̃T)]i,i.
Note now that(212)Πi[(I-X̃)A(I-X̃T)]i,i,
does not depend on di, so we can minimize this factor independently of di.
Consider then the general weighted Frobenius norm minimization problem(213)min⁡G∈Str⁡(I-GB)W(I-GB)T.
As we have seen, its solution G satisfies the relation
(214)(GBWBT)i,j=(WBT)i,j,∀(i,j)∈S.
Let nowG=X̃,   W=A,   B=I,   S=SL̃. Then,
(215)(GBWBT)ij=(WBT)ij
takes the form
(216)(X̃A)ij=Aij∀i,j∈SL̃.
This is an explicit method and since the minimization is done rowwise it follows from (213), with the chosen matrices G,  B and W, that each of
(217)[(I-X̃)A(I-X̃)T]i,ii=1,…,n
is minimized separately. By construction L̃ satisfies (216), so the minimization problem is has the solution X̃=L̃. Hence,
(218)min⁡X̃Πi[(I-X̃)A(I-X̃)T]i,i=Πi[(I-L̃)A(I-L̃)T]i,i.
Consider next the first factor in (211). Here,
(219)((1/n)∑jαj2)nΠjαj2≥1,
since a geometric average is less or equal to an arithmetic average. Equality is taken if and only if all αj are equal and with no limitation we can take αj=1,   j=1,…,n. Hence,
(220)di2=1[(I-L̃)A(I-L̃)T]i,i
which completes the proof.The method above provides a simple and cheap method to compute approximate inverses on factorized form. The proof of the theorem shows that theK-condition number is reduced in a way as follows from the next corollary.Corollary 2.
LetL̂,  D be defined as in Theorem 12. Then,
(221)K(L̂AL̂T)≡((1/n)tr⁡(L̂AL̂T))ndet⁡(L̂AL̂T)=Π1n[(I-L̃)A(I-L̃T)]iidet⁡(A),
where L̃=I-D-1L.Hence, the trace is replaced by a product, that is then′th power of the arithmetic average is replaced the n′th power of a geometric average. This is illustrated in the next example.Example 6.
LetSL={(1,1),(2,2),…,(n,n)}, that is, let L be a diagonal matrix. Then, we find di2=aii and Corollary 2 shows that
(222)K(LALT)=min⁡L∈SL((1/n)tr⁡(LALT))ndet⁡(LALT)=Π1naiidet⁡(A),
which is to be compared with
(223)K(A)=((1/n)tr⁡(A))ndet⁡(A)=((1/n)∑1naii)ndet⁡(A),
that is, we have
(224)K(LALT)=(g̲a̲)nK(A).
Hence, the K-condition number K(LALT) of the diagonally scaled matrix LALT is substantially smaller than K(A) if the geometric average g̲ of the diagonal entries aii of A are much smaller than their arithmetic average a¯. This holds when the entries aii vary significantly. Note that it always holds that g̲≤a¯.
We conclude this section by mentioning that theK-condition number can be take large values even for seemingly harmless eigenvalue distributions.Example 7 (Arithmetic distribution).
Let0<λ1<λ2<⋯<λn, the eigenvalues of B be distributed uniformly as an arithmetic sequence in the interval [a,b],a=λ1, b=λn. For simplicity, assume that n/2 is even. Then,
(225)K(B)=(b+a)/2)n∏1nλi.
On the other hand, (57) shows k≤(1/2)b/aln2/ɛ, which is asymptotically smaller than n if (b/a)n-2=o(1). In particular, if b/a does not depend on n then we have k≤O(ln1/ɛ). Therefore, the estimate in Theorem 6 inferior even to the simple estimate in (56). For other distributions, however, Theorem 6 can give a smaller upper bound.


## 5.1. Explicit Methods

In these methods, an approximation of the inverseA-1 of a given nonsingular matrix A is computed explicitly, that is, without solving a linear globally coupled system of equations.LetS be a sparsity pattern. We want to compute G∈S, such that(190)(GA)ij=δij,(i,j)∈S,that is(191)∑k:(i,k)∈Sgikakj=δij,(i,j)∈S.
Some observations can be made from (191):(i)
the elements in each row ofG can be computed independently;(ii)
even ifA is symmetric, G in not necessarily symmetric, because gi,j,  j≠i,  and gj,i are, in general, not equal.


## 5.2. Implicit Methods

These methods require thatA is factored first. In practice, they are used mainly for band or “envelope" matrices. The algorithm was presented in [79]. It is based on an idea in [80]; see also [81].Suppose thatA=LD-1U is a triangular matrix factorization of A. If A is a band matrix then L and U are also band matrices.Let(192)L=I-L̃,U=I-Ũ,where L̃ and Ũ are strictly lower and upper triangular matrices correspondingly.The following lemma holds.Lemma 5.
Using the above notations it holds that(i)
A-1=DL-1+ŨA-1,(ii)
A-1=U-1D+A-1L̃.Proof.
Consider the following(193)A=LD-1U⟹A-1=U-1DL-1⟹(I-Ũ)A-1=DL-1⟹A-1=DL-1+ŨA-1.
Also,
(194)A-1(I-L̃)=U-1D⟹A-1=U-1D+A-1L̃.SinceDL-1 is lower triangular and Ũ is upper triangular, using (i) we can compute entries in the upper triangular part of A-1 with no need to use entries of L-1. Similarly, using (ii) we can compute entries of the lower triangular part A-1 without computing U-1.Suppose now thatA is a block banded matrix with a semibandwidth p, and we want to form A-1 also as block banded with a semibandwidth q:q≥p. The identities (i) and (ii) can be used then for the computation of the upper and lower parts of A-1.Remark 3.
(i) The algorithm involves onlymatrix×matrix operations.
(ii) There is no need to compute any entries outside the bands.
(iii) IfA is symmetric then it suffices executing only (i) or (ii).
(iv) It can be seen that(A-1)nn=Dnn-1.
There are two drawbacks with the above algorithm. It requires first the factorizationA=LD-1U and even if A is s.p.d, the band matrix part of A-1, which is computed, need not be s.p.d. The next example illustrates this.Example 2.
Consider an s.p.d. matrix(195)G=[1-21-25-31-34],Gband=[1-20-25-30-34],
is indefinite.


## 5.3. A General Framework for Computing Approximate Inverses

It turns out that both the explicit and implicit method can be characterized as methods to compute best approximations ofA-1 of all matrices having a given sparsity pattern, in some norm. The basic idea is due to Kolotilina and Yeremin  [38, 79], see [6]. Recall that the trace function is defined by tr(A)=∑i=1naii, which also equals ∑i=1nλi(A). Let a sparsity pattern S be given. Consider the functional(196)FW(G)≡‖I-GA‖W2=tr⁡((I-GA)W(I-GA)T),where the weight matrix W is s.p.d. If W≡I then ∥I-GA∥I is the Frobenius norm of I-GA.ClearlyFW(G)≥0. If G=A-1 then FW(G)=0. We want to compute the entries of G in order to minimize FW(G), that is, to find Ĝ∈S, such that(197)‖I-ĜA‖W≤‖I-GA‖W,∀G∈S.The following properties of the trace function will be used(198)tr⁡A=tr⁡AT,tr⁡(A+B)=tr⁡A+tr⁡B.Then,(199)FW(G)=tr⁡(I-GA)W(I-GA)T=tr⁡(W-GAW-W(GA)T+GAW(GA)T)=tr⁡W-tr⁡GAW-tr⁡(GAW)T+tr⁡GAWATGT.Further, as we are interested in minimizingFW with respect to G∈S, we consider the entries gi,j as variables. The necessary condition for a minimizing point are then(200)∂FW(G)∂gij=0,(i,j)∈S.
From (199) and (200), we get(201)-2(WAT)ij+2(GAWAT)ij=0,or(202)(GAWAT)ij=(WAT)ij,(i,j)∈S.
Depending on the particular matrix A and the choice of S and W,  (202) may or may not have a solution. We give some examples where a solution exists.Example 3.
LetA be s.p.d. Choose W=A-1 which is also s.p.d. Then, (202) implies
(203)(GA)ij=δij,(i,j)∈S,
which is the formula for the previously presented explicit method which, hence, is a special case of the more general framework for computing approximate inverses using weighted Frobenius norm.Example 4.
LetW=(ATA)-1. Then (202) implies
(204)(G)ij=(A-1)ij,(i,j)∈S,
which is the relation for the previously presented implicit method. In this case the entries of G are the corresponding entries of the exact inverse.Example 5.
LetW=I. Then,
(205)FW(G)=n-tr⁡(GA),(GAAT)ij=(AT)ij,(i,j)∈S.
This method is also explicit.
We can expect that such methods will be accurate only if all elements ofA which are not used in the computations are zero or are relatively small. In some cases the quality of the computed approximation G to A-1 can be significantly improved using diagonal compensation of the entries of A which are outside S. The best approximation G to A-1 in a (weighted) Frobenius norm is in general not symmetric and, as we have seen, not always positive definite. For this reason, the next, alternate method, is considered.


## 5.4. Constructing a Symmetric and Positive Definite Approximate Inverse

For some methods (as in the preconditioned Chebyshev and the conjugate gradient iteration methods) it is of importance to use s.p.d. preconditioners. As we have seen, the methods described till now do not guarantee thatG will be such a matrix.In order to compute an s.p.d. approximate inverse of an s.p.d. matrix, we can proceed as follows. It will be shown that this approximation gives a best approximation to minimize theK-condition number of the correspondingly preconditioned matrix.A Symmetric and Positive Definite Approximate Factorized Inverse
Seek an approximate inverse in the formG=LLT, where L∈SL,
(206)SL={(i,j)∈S,i≥j},
and let
(207)SL̃={(i,j)∈SL,i>j},
that is, denote by SL and SL̃ the lower and strictly lower triangular part of the sparsity set S.Theorem 12.
LetA be s.p.d. and consider matrices L with sparsity pattern Sl. Let the matrix L̂ be computed by the following steps.(i)
Compute firstL̃ such that
(208)(L̃A)ij=Aij,(i,j)∈SL̃,
and L̃ij=0,  (i,j)∈SL̃c (the complement set).(ii)
LetL̂=D(I-L̃), where
(209)D=diag⁡(d1,d2,…,dn),di=1[(I-L̃)A(I-L̃T)]i,i1/2.
Then, L̂∈SL and minimizes the K-condition number of L̂AL̂T, that is,
(210)((1/n)tr⁡(L̂AL̂T))ndet⁡(L̂AL̂T)=inf⁡L∈SL((1/n)tr⁡(LALT))ndet⁡(LALT).Proof.
LetD=diag(d1,…,dn) denote the diagonal part of a matrix X∈SL and let X̃=I-D-1X, that is, X̃∈SL̃. Then,
(211)((1/n)tr⁡(XAXT))ndet⁡(XAXT)=((1/n)∑i(XAXT)ii)n(det⁡(X))2det⁡(A)=((1/n)∑i[D(I-X̃)A(I-X̃T)D]i,i)n(det⁡(X))2det⁡(A)=((1/n)∑idi2[(I-X̃)A(I-X̃T)]i,i)n(Πidi2)det⁡(A)=((1/n)∑iαi2)nΠiαi2⋅Πi[(I-X̃)A(I-X̃T)]i,idet⁡(A),
where αi2=di2[(I-X̃)A(I-X̃T)]i,i.
Note now that(212)Πi[(I-X̃)A(I-X̃T)]i,i,
does not depend on di, so we can minimize this factor independently of di.
Consider then the general weighted Frobenius norm minimization problem(213)min⁡G∈Str⁡(I-GB)W(I-GB)T.
As we have seen, its solution G satisfies the relation
(214)(GBWBT)i,j=(WBT)i,j,∀(i,j)∈S.
Let nowG=X̃,   W=A,   B=I,   S=SL̃. Then,
(215)(GBWBT)ij=(WBT)ij
takes the form
(216)(X̃A)ij=Aij∀i,j∈SL̃.
This is an explicit method and since the minimization is done rowwise it follows from (213), with the chosen matrices G,  B and W, that each of
(217)[(I-X̃)A(I-X̃)T]i,ii=1,…,n
is minimized separately. By construction L̃ satisfies (216), so the minimization problem is has the solution X̃=L̃. Hence,
(218)min⁡X̃Πi[(I-X̃)A(I-X̃)T]i,i=Πi[(I-L̃)A(I-L̃)T]i,i.
Consider next the first factor in (211). Here,
(219)((1/n)∑jαj2)nΠjαj2≥1,
since a geometric average is less or equal to an arithmetic average. Equality is taken if and only if all αj are equal and with no limitation we can take αj=1,   j=1,…,n. Hence,
(220)di2=1[(I-L̃)A(I-L̃)T]i,i
which completes the proof.The method above provides a simple and cheap method to compute approximate inverses on factorized form. The proof of the theorem shows that theK-condition number is reduced in a way as follows from the next corollary.Corollary 2.
LetL̂,  D be defined as in Theorem 12. Then,
(221)K(L̂AL̂T)≡((1/n)tr⁡(L̂AL̂T))ndet⁡(L̂AL̂T)=Π1n[(I-L̃)A(I-L̃T)]iidet⁡(A),
where L̃=I-D-1L.Hence, the trace is replaced by a product, that is then′th power of the arithmetic average is replaced the n′th power of a geometric average. This is illustrated in the next example.Example 6.
LetSL={(1,1),(2,2),…,(n,n)}, that is, let L be a diagonal matrix. Then, we find di2=aii and Corollary 2 shows that
(222)K(LALT)=min⁡L∈SL((1/n)tr⁡(LALT))ndet⁡(LALT)=Π1naiidet⁡(A),
which is to be compared with
(223)K(A)=((1/n)tr⁡(A))ndet⁡(A)=((1/n)∑1naii)ndet⁡(A),
that is, we have
(224)K(LALT)=(g̲a̲)nK(A).
Hence, the K-condition number K(LALT) of the diagonally scaled matrix LALT is substantially smaller than K(A) if the geometric average g̲ of the diagonal entries aii of A are much smaller than their arithmetic average a¯. This holds when the entries aii vary significantly. Note that it always holds that g̲≤a¯.
We conclude this section by mentioning that theK-condition number can be take large values even for seemingly harmless eigenvalue distributions.Example 7 (Arithmetic distribution).
Let0<λ1<λ2<⋯<λn, the eigenvalues of B be distributed uniformly as an arithmetic sequence in the interval [a,b],a=λ1, b=λn. For simplicity, assume that n/2 is even. Then,
(225)K(B)=(b+a)/2)n∏1nλi.
On the other hand, (57) shows k≤(1/2)b/aln2/ɛ, which is asymptotically smaller than n if (b/a)n-2=o(1). In particular, if b/a does not depend on n then we have k≤O(ln1/ɛ). Therefore, the estimate in Theorem 6 inferior even to the simple estimate in (56). For other distributions, however, Theorem 6 can give a smaller upper bound.


## 6. Augmented Subspace  Preconditioning Method

### 6.1. Introduction; Preconditioners for Very Ill-Conditioned Problems

In this section, we consider the solution of systemsAx=b, where Ais an n×n matrix which is symmetric and positive definite (s.p.d) and can have a very large condition number, that is, be ill-conditioned. Such systems arise typically for near-limit values of some problem parameter. (Ratio of material coefficients, aspect ratio of the domain, nearly incompressible materials in elasticity theory, etc.) The condition number can be additionally very large due to the size of the matrix A (a small value of the discretization parameter) and also due to an irregular mesh and/or large aspect ratios of the mesh in partial differential equation (PDE) problems.If the size of the system is not too large one can use direct solution methods, possibly coupled with an iterative refinement method.LetB=LDLT ( or B=L̃L̃T ) be a triangular matrix or the Cholesky factorization of A. Due to finite precision computations (say, in single precision) in general B is only an approximation of A. The iterative refinement method takes the following form.Algorithm 1 (Iterative refinement method).
Givenx(0)=0
fork=0,1,2,…, until convergence
(i)
computer(k)=b-Ax(k),(ii)
solveBd(k)=rk,(iii)
letx(k+1)=x(k)+d(k), and repeat (i)–(iii).Frequently, it suffices with one iterative refinement step. The ability of iterative refinement to produce a more accurate solution vector depends crucially on how the computation of the residual vectorr(k) in (i) is implemented. A safe way is to use double precision for this computation but possibly single precision in (ii) and (iii). However, as described in [30], if one rewrites the computation of Ax(k) as a sum of differences, in some cases it suffices to use single precision in (i) also.The computational labor is normally dominated by the initial factorization ofA. For large systems this cost can become too big as it grows in general fast with problem size. (For an elliptic difference problem on a 3D  N×N×N mesh it grows as O(N7) for certain band-matrix orderings. Furthermore, the demand of memory to store the factor L grows as O(N5). For certain nested dissection and other orderings, the complexity is somewhat reduced, however.)Therefore, iterative solution methods become the ultimate methods of choice. As we have seen in Section1, the basic idea behind the iterative solution technique is to use a cheaper (incomplete) factorization or other approximation B of A and to compensate for this approximation by repeating the steps in the iterative refinement method until the residual is sufficiently small. In addition, to speed up the convergence of the method, one or more acceleration parameters are introduced, for instance, (iii) becomes x(k+1)=x(k)+αkd(k) for certain parameters αk.By a proper choice ofB the number of required iterations may not be too big while the expense in solving systems with matrix B may not be larger than the order of some “work units", for instance, can correspond to a few actions (matrix-vector multiplications) of Aon a vector. In this way, one can gain significantly in computational labor and less demand of memory resources as compared with a direct solver. Actually, the direct solver can be viewed as an approximate factorization with the full amount of fill-in allowed, while as we have seen, in a incomplete factorization method one controls the amount of fill-in either by using a predetermined sparsity pattern in L or by allowing a variable pattern, which depends on some relative drop tolerance. (Such a drop tolerance is to delete a fill-in entry aij if there holds |aij|≤ɛaiiajj,  j≠i for some ɛ,   0<ɛ<1. More details can be found in [82]).A problem with iterative solution methods for ill-conditioned systems is that they may stagnate, that is, there is no further improvement as the method proceeds. This occurs typically for minimum residual or minimumA-norm methods. For other type of methods even divergence may be observed. Another problematic issue is the fact that if the residual norm has taken a small value, this does not necessarily mean that the error norm is sufficiently small, since(226)‖x-x(k)‖2=‖A-1A(x-x(k))‖2≤‖A-1‖2‖r(k)‖2=1λmin⁡(A)‖r(k)‖2,and here λmin(A) takes very small values for ill-conditioned systems. Hence, even if ∥r(k)∥ is small, ∥x-x(k)∥ may still be large. For ill-conditioned systems one sees then typically a reduction of the residual to some limit value while the errors hardly decay at all. This was illustrated in Section 3. For studies on the influence of inexact arithmetics, see for example [83–85].This situation can be significantly improved by using a proper preconditioner. Then,(227)‖x-x(k)‖2=‖(B-1A)-1B-1A(x-x(k))‖2≤1λmin⁡(B-1A)‖r̃(k)‖2,where r̃(k)=B-1A(x-x(k))=B-1r(k) is the so called preconditionedor pseudo-residual. Here λmin(B-1A)≫λmin(A) with a proper preconditioner. Therefore, the importance of choosing a proper preconditioner is twofold:  (1)
to increase the rate of convergence while keeping the expense in solving systems withB low, and(2)
to enable a small error norm when the pseudoresidual is small.Preconditioning methods, such as the modified incomplete factoriztion method, multigrid and multilevel methods, aim at reducing arror components correponding both to the large eigenvalues with rapidly oscillating components and the smaller eigenvalues for smoother eigen functions. In the modified method, this is partly achieved by letting the preconditioner by exact for a particular smooth component of the solution, such as for the constant component vector. It has been shown, see [6, 47] when applied for elliptic difference problems, that under certain conditions the spectral condition number is reduced from O(h-2) to O(h-1). In multigrid methods, one works on two or more levels of meshes where the finer grid component should smooth out the fast, oscillating components in the iteration error, while the coarser mesh should handle the smooth components. Under certain conditions, such methods way reduce the above condition number to optimal order, O(1), as h→0.The multigrid method was first introduced for finite difference methods in the 1960s by Fedorenko [86], and Bakhvalov [87], and further developed and advocated by Brandt in the 1970s, see, for example, Brandt [88]. For finite elements it has been pursued by, for example, Braess [89], Hackbusch [90], Bramble et al. [91], Mandel  et al. [92], McCormick [57], Bramble et al. [93] and Bank et al. [94], among others.As it turns out, such standard preconditioning methods, namely (modified) incomplete factorization ((M)ILU), [46, 47], Multigrid (MG) [90], or Algebraic Multilevel Iteration (AMLI), [95–97], methods may not be efficient in both and in particular, in the second of the above mentioned requirements. This might be due to the fact that the smallest eigenvalue (in the preconditioned system) is caused by some problem parameter which these methods leave unaffected. Therefore there is a demand for new types of preconditioners (or new combinations of already known preconditioners). To satisfy the above need, two types of preconditioners have been constructed:(a)
deflation methods,(b)
augmented matrix methods,which we now describe.


### 6.2. Deflation Methods

The deflation technique is based on a projection matrix. Assume thatA has a number of (very) small eigenvalues, say m̃,  0<λ1≤λ2…≤λm̃, and let 𝒲={w(i)},   i=1,…,m̃ be their corresponding eigenvectors (Awi=λiwi). Let V be a rectangular matrix of order n×m, where m<n (in practice m≪n) of full rank, where the m columns of V span a subspace γ, such that Imγ contains the eigenvectors corresponding to the  “bad" subspace 𝒲. Hence, m≥m̃.Lemma 6.
LetP=AVAV-1VT, where AV=VTAV. Then, the following holds:(a)
P2=P, that is a projector;(b)
P(AV)=AV;(c)
(I-P)b=0 if b∈ℐm(AV);(d)
PTV=V;(e)
(I-P)A is symmetric and positive definite and has a nullspace of dimension m.Proof.
Note first thatAV is nonsingular since V has a full rank (=m). The statements follow now by straightforward computations.Lemma6 shows that P is projection matrix which maps any vector onto AV. Similarly, PT is a projection matrix which maps Vonto itself. We will use the matrix Pin three slightly different ways to solve ill-conditioned systemsWe split first the right-side vectorb in two components:(228)b=Pb+(I-P)b.(These components areA-1 orthogonal, i.e., (Pb)TA-1(I-P)b=0.) The first splits the computation of thesolution vector corresponding.Method 1 (Splitting of the solution vector).
Let(229)x(0)=VAV-1VTb.Then,(230)Ax(0)=Pb.Solve(231)Az=(I-B)b.The solutionx of Ax=bis then(232)x=x(0)+z.
Here, x(0) and z are A-orthogonal.
Note thatAz=b-Ax(0). The matrix AV is normally of small order and the arising system in (229) can be solved with relatively little expense using a direct solution method. Furthermore, the system (231) is well-conditioned on the solution subspace, because, as follows from part (c) of Lemma 6, (I-P)b, and hence z do not contain components of any of the first m  “small" eigenvectors wi,  i=1,2,…,m. Hence, (231) can be solved by the CG method with a rate of convergence determined by the effective condition numberλn/λm+1, which is expected to be substantially smaller than λn/λ1.
However, the method requires exact solution of systems withAV and for some problems m is not that small. Also, it is assumed that the projection Pb is computed exactly (or to a sufficient accuracy), which may be unfeasible in many applications.Method 2 (Defect-correction with projectors).
In the presence of round-off errors,x(0) may not be sufficiently accurate and b-Ax(0) may still contain components in the “bad" subspace. A defect-correction (iterative refinement) procedure may then help.  Let  x(0)=0 for    k=0,1,2,…, until convergence. Compute  r(k)=b-Ax(k).  Solve  Ad(k)=b-Ax(k) as follows:(i)
d(k,0)=VAV-1VTr(k),(ii)
Az(k)=(I-P)r(k), or   Ay(k)=r(k)-Ad(k,0),(iii)
d(k)=d(k,0)+y(k).Let  x(k+1)=x(k)+d(k).
In this method, it normally suffices with few defect-correction steps.
For some extremely ill-conditioned systems, the implementation of the defect-correction method as a preconditioning method may be necessary. Note then that as follows from Lemma6, (I-P)A is symmetric so the standard conjugate gradient method can be used.Method 3 (Preconditioning by a projection matrix).
Letx(0)=VAV-1VTb.  Solve   (I-P)Az=(I-P)b,  or  (I-P)Az=b-Ax(0) by CG iteration.  Let  x=z+x(0).
Note thatx(0) is contained in the null-space of (I-P)A, since (I-P)Ax(0)=(I-P)Pb=(P2-P)b=0. Here, the system (I-P)A is well-conditioned on the orthogonal complement to the null-space and, in addition, the right-hand-side has no or only small components in the bad subspace.
Methods1, 2, and 3 require accurate solution of systems with the matrix AV. It is a viable step for small values of m. However, when the dimension of the “bad" subspace of A is relatively big, it may be too costly. Furthermore, the iteration Method 3 involves two multiplications with A (one involved in P and one required to computeAz(k)) at each iteration step when computing the search direction vectors and is therefore particularly expensive.
Another issue to comment on is that the above methods are assumed to move the components of the eigenvectors for the smallest eigenvalues ofA to become exactly zero. However, this can be sensitive to perturbations and occurs only in exact arithmetic. As we have seen, it is a viable method for small dimensions of the subspace causing the ill-conditioning but it may be inefficient for larger dimensions.
Deflation methods have been used and analysed by [98–100], among others.
In the next section and Section4, we present a method which move the small eigenvalues to the cluster of bigger eigenvalues which is much less dependent on having the right subspace spanned by the columns of V and which do not require exact solution of systems with AV.


### 6.3. Augmented Matrix Preconditioning Methods, the Ideal Case

We now present an alternative method to handle ill-posed problems. In this method the small eigenvalues are moved to the cluster of bigger eigenvalues, instead of being deflated to zero, as in the deflation method. The method is an extension of the method presented in [101]. The presentation here is based in [25, 102]. First, we consider B=I+VVT as a (multiplicative) preconditioner to A.In this case one must scale the column vectors appearing inV properly. A method involving an automatic scaling is based on a projection matrix. Let then(233)B=I+σVAV-1VT,AV=VTAV.Let(λi,vi)i=1m be the eigenpairs of A for the smallest eigenvalues, 0<λ1≤λ2≤⋯λm. If V=[v1,…,vm], we get then(234)λ̃i=λi(BA)={λi+σ,1≤i≤m,λi,m+1≤i≤n.
Hence, σ determines how much the smallest eigenvalues are moved. If λm+1≤λ1+σ and λm+σ≤λn, then λm+1≤λ̃i≤λn, that is, the m smallest eigenvalues have been moved to the cluster [λm+1,λn] of bigger eigenvalues and the spectral condition number of BAis κ(BA)=λn/λm+1, which normally means a significant reduction, compared to κ(A)=λn/λ1.The above illustrates what can be achieved in an ideal case. In practice, the exact eigenvectors{v1,…,vm} (or the subspace spanned by them) are not known. Even if the eigenvectors are known it is not efficient to use them to form the matrix V because they are in general not sparse and the matrix vector multiplications with V will be costly. Hence, in practice other vectors which are sparse but spans about the same subspace as the smoothest eigenvectors must be used; otherwise the expense of the preconditioner would be too high. We consider therefore more general subspaces spanned by the column vectors of V. The next lemma will be useful.Lemma 7.
LetA be s.p.d. Then,
(235)P=A1/2V(VTAV)-1VTA1/2,
is an orthogonal projection, that is, P2=P and P*=P. Therefore, the only eigenvalues of P are 0 and 1.Proof.
Consider the following(236)P2=A1/2V(VTAV)-1VTA1/2A1/2V(VTAV)-1VTA1/2=A1/2V(VTAV)-1VTA1/2=P.The next theorem shows (what can also be expected) that the clustering can never get worse for expanding subspaces spanned by the column-vectors ofV, that is, there holds a monotonicity principle.Theorem 13 (see [103]).
LetA and Â be s.p.d. matrices of order n×n and let   Vk  be rectangular matrices of order n×mk,  k=1,2 such that rankVk=mk,  k=1,2. If ImV1⊆ImV2, then for all i,   1≤i≤n the following inequality holds
(237)λi((I+V2(V2TÂV2)-1V2T)A)≥λi((I+V1(V1TÂV1)-1V1T)A).Proof.
It is readily seen that the proposition holds ifF=V2(V2TÂV2)-1V2T-V1(V1TÂV1)-1V1T is negative definite. But since ImV1⊆ImV2, there exists some matrix Q of order m2×m1 such that V1=V2Q. Then, with DK=VkTÂVk, we have
(238)F=V2(D2-1-QD1-1QT)V2T=V2D2-1/2(I-D2-1/2QD1-1/2QTD2-1/2)D2-1/2V2T,
where
(239)P≡D2-1/2QD1-1/2QTD2-1/2=D2-1/2Q(QTD2Q)-1QTD2-1/2,
is an orthogonal projector (P2=P), whose eigenvalues are 0 and 1.Corollary 3.
IfImV1=ImV2 then I+V2D2-1V2T=I+V1D1-1V1T.Proof.
In this case,Q in V1=V2Q is invertible. Thus, D2-1/2Q(QTD2Q)-1QTD2-1/2=I.Remark 4.
The above corollary shows that the individual eigenvectors ofA are not needed when constructing the matrix V; we are rather interested in the subspace spanned by them.
The most interesting case for us is whenImV⊃span{v1,…,vm}, where vi are the eigenvalues of Afor the smallest eigenvalues λ1,…,λm. Then, the preconditioner B=I+σVAV-1VT moves the smallest eigenvalues λi at least to λi+σ, where σis a scaling parameter.Theorem 14.
LetB=I+σVAV-1VT and let  λ1,…,λm be the smallest eigenvalues of A. Then, for the eigenvalues of BA there holds. 									(240)λ̃i≥{λi+σ,1≤i≤mλi,m+1≤i≤n.Proof.
LetV1=[v1,…,vm], where {vi}1m are the first m eigenvalues of A and let V2=V. Then, Lemma 7 and (234) show the result.It may happen that the eigenvalues are moved too far so that the maximum eigenvalue ofBA is much larger that of A.Theorem 15.
LetA be s.p.d. of order n×n and let the rectangular matrix V of order n×m be defined as V=[v1,v2,…,vm]. Assume that rankV=m. Further, define Ã=(I+σVAV-1VT)A, where AV=VTAV. Then,
(241)λmax⁡(Ã)≤σ+λmax⁡(A).Proof.
The result follows from the following relations:(242)λmax⁡(Ã)≤λmax⁡(A)+σsup⁡xTV(VTAVV)-1VTxxTA-1x=λmax⁡(A)+σsup⁡xTA1/2V(VTAV)-1VTA1/2xxTx=λmax⁡(A)+σ,
where the last equality follows from Lemma 7It follows from Theorems 6.3 and16 that the optimal value of σ=λm+1, in which case κ(BA)≤(λn+σ)/λm+1. In general, λm+1 may not be known. With σ=λn, we obtain κ(BA)≤(2λn)/λm+1.The moral we can draw from the above is that the suggested technique can be a very useful means to reduce the condition number of a given matrixA if we have information about the eigenvectors corresponding to the smallest eigenvalues of A. Since in practice this is hardly ever the case, a natural step to undertake is to consider not the individual eigenvectors but the subspace spanned by some approximation of them.In the next section, we present a generalized form of the augmented matrix preconditioner which allows for both approximate subspaces and the replacement ofAV by a simpler matrix BV.


## 6.1. Introduction; Preconditioners for Very Ill-Conditioned Problems

In this section, we consider the solution of systemsAx=b, where Ais an n×n matrix which is symmetric and positive definite (s.p.d) and can have a very large condition number, that is, be ill-conditioned. Such systems arise typically for near-limit values of some problem parameter. (Ratio of material coefficients, aspect ratio of the domain, nearly incompressible materials in elasticity theory, etc.) The condition number can be additionally very large due to the size of the matrix A (a small value of the discretization parameter) and also due to an irregular mesh and/or large aspect ratios of the mesh in partial differential equation (PDE) problems.If the size of the system is not too large one can use direct solution methods, possibly coupled with an iterative refinement method.LetB=LDLT ( or B=L̃L̃T ) be a triangular matrix or the Cholesky factorization of A. Due to finite precision computations (say, in single precision) in general B is only an approximation of A. The iterative refinement method takes the following form.Algorithm 1 (Iterative refinement method).
Givenx(0)=0
fork=0,1,2,…, until convergence
(i)
computer(k)=b-Ax(k),(ii)
solveBd(k)=rk,(iii)
letx(k+1)=x(k)+d(k), and repeat (i)–(iii).Frequently, it suffices with one iterative refinement step. The ability of iterative refinement to produce a more accurate solution vector depends crucially on how the computation of the residual vectorr(k) in (i) is implemented. A safe way is to use double precision for this computation but possibly single precision in (ii) and (iii). However, as described in [30], if one rewrites the computation of Ax(k) as a sum of differences, in some cases it suffices to use single precision in (i) also.The computational labor is normally dominated by the initial factorization ofA. For large systems this cost can become too big as it grows in general fast with problem size. (For an elliptic difference problem on a 3D  N×N×N mesh it grows as O(N7) for certain band-matrix orderings. Furthermore, the demand of memory to store the factor L grows as O(N5). For certain nested dissection and other orderings, the complexity is somewhat reduced, however.)Therefore, iterative solution methods become the ultimate methods of choice. As we have seen in Section1, the basic idea behind the iterative solution technique is to use a cheaper (incomplete) factorization or other approximation B of A and to compensate for this approximation by repeating the steps in the iterative refinement method until the residual is sufficiently small. In addition, to speed up the convergence of the method, one or more acceleration parameters are introduced, for instance, (iii) becomes x(k+1)=x(k)+αkd(k) for certain parameters αk.By a proper choice ofB the number of required iterations may not be too big while the expense in solving systems with matrix B may not be larger than the order of some “work units", for instance, can correspond to a few actions (matrix-vector multiplications) of Aon a vector. In this way, one can gain significantly in computational labor and less demand of memory resources as compared with a direct solver. Actually, the direct solver can be viewed as an approximate factorization with the full amount of fill-in allowed, while as we have seen, in a incomplete factorization method one controls the amount of fill-in either by using a predetermined sparsity pattern in L or by allowing a variable pattern, which depends on some relative drop tolerance. (Such a drop tolerance is to delete a fill-in entry aij if there holds |aij|≤ɛaiiajj,  j≠i for some ɛ,   0<ɛ<1. More details can be found in [82]).A problem with iterative solution methods for ill-conditioned systems is that they may stagnate, that is, there is no further improvement as the method proceeds. This occurs typically for minimum residual or minimumA-norm methods. For other type of methods even divergence may be observed. Another problematic issue is the fact that if the residual norm has taken a small value, this does not necessarily mean that the error norm is sufficiently small, since(226)‖x-x(k)‖2=‖A-1A(x-x(k))‖2≤‖A-1‖2‖r(k)‖2=1λmin⁡(A)‖r(k)‖2,and here λmin(A) takes very small values for ill-conditioned systems. Hence, even if ∥r(k)∥ is small, ∥x-x(k)∥ may still be large. For ill-conditioned systems one sees then typically a reduction of the residual to some limit value while the errors hardly decay at all. This was illustrated in Section 3. For studies on the influence of inexact arithmetics, see for example [83–85].This situation can be significantly improved by using a proper preconditioner. Then,(227)‖x-x(k)‖2=‖(B-1A)-1B-1A(x-x(k))‖2≤1λmin⁡(B-1A)‖r̃(k)‖2,where r̃(k)=B-1A(x-x(k))=B-1r(k) is the so called preconditionedor pseudo-residual. Here λmin(B-1A)≫λmin(A) with a proper preconditioner. Therefore, the importance of choosing a proper preconditioner is twofold:  (1)
to increase the rate of convergence while keeping the expense in solving systems withB low, and(2)
to enable a small error norm when the pseudoresidual is small.Preconditioning methods, such as the modified incomplete factoriztion method, multigrid and multilevel methods, aim at reducing arror components correponding both to the large eigenvalues with rapidly oscillating components and the smaller eigenvalues for smoother eigen functions. In the modified method, this is partly achieved by letting the preconditioner by exact for a particular smooth component of the solution, such as for the constant component vector. It has been shown, see [6, 47] when applied for elliptic difference problems, that under certain conditions the spectral condition number is reduced from O(h-2) to O(h-1). In multigrid methods, one works on two or more levels of meshes where the finer grid component should smooth out the fast, oscillating components in the iteration error, while the coarser mesh should handle the smooth components. Under certain conditions, such methods way reduce the above condition number to optimal order, O(1), as h→0.The multigrid method was first introduced for finite difference methods in the 1960s by Fedorenko [86], and Bakhvalov [87], and further developed and advocated by Brandt in the 1970s, see, for example, Brandt [88]. For finite elements it has been pursued by, for example, Braess [89], Hackbusch [90], Bramble et al. [91], Mandel  et al. [92], McCormick [57], Bramble et al. [93] and Bank et al. [94], among others.As it turns out, such standard preconditioning methods, namely (modified) incomplete factorization ((M)ILU), [46, 47], Multigrid (MG) [90], or Algebraic Multilevel Iteration (AMLI), [95–97], methods may not be efficient in both and in particular, in the second of the above mentioned requirements. This might be due to the fact that the smallest eigenvalue (in the preconditioned system) is caused by some problem parameter which these methods leave unaffected. Therefore there is a demand for new types of preconditioners (or new combinations of already known preconditioners). To satisfy the above need, two types of preconditioners have been constructed:(a)
deflation methods,(b)
augmented matrix methods,which we now describe.


## 6.2. Deflation Methods

The deflation technique is based on a projection matrix. Assume thatA has a number of (very) small eigenvalues, say m̃,  0<λ1≤λ2…≤λm̃, and let 𝒲={w(i)},   i=1,…,m̃ be their corresponding eigenvectors (Awi=λiwi). Let V be a rectangular matrix of order n×m, where m<n (in practice m≪n) of full rank, where the m columns of V span a subspace γ, such that Imγ contains the eigenvectors corresponding to the  “bad" subspace 𝒲. Hence, m≥m̃.Lemma 6.
LetP=AVAV-1VT, where AV=VTAV. Then, the following holds:(a)
P2=P, that is a projector;(b)
P(AV)=AV;(c)
(I-P)b=0 if b∈ℐm(AV);(d)
PTV=V;(e)
(I-P)A is symmetric and positive definite and has a nullspace of dimension m.Proof.
Note first thatAV is nonsingular since V has a full rank (=m). The statements follow now by straightforward computations.Lemma6 shows that P is projection matrix which maps any vector onto AV. Similarly, PT is a projection matrix which maps Vonto itself. We will use the matrix Pin three slightly different ways to solve ill-conditioned systemsWe split first the right-side vectorb in two components:(228)b=Pb+(I-P)b.(These components areA-1 orthogonal, i.e., (Pb)TA-1(I-P)b=0.) The first splits the computation of thesolution vector corresponding.Method 1 (Splitting of the solution vector).
Let(229)x(0)=VAV-1VTb.Then,(230)Ax(0)=Pb.Solve(231)Az=(I-B)b.The solutionx of Ax=bis then(232)x=x(0)+z.
Here, x(0) and z are A-orthogonal.
Note thatAz=b-Ax(0). The matrix AV is normally of small order and the arising system in (229) can be solved with relatively little expense using a direct solution method. Furthermore, the system (231) is well-conditioned on the solution subspace, because, as follows from part (c) of Lemma 6, (I-P)b, and hence z do not contain components of any of the first m  “small" eigenvectors wi,  i=1,2,…,m. Hence, (231) can be solved by the CG method with a rate of convergence determined by the effective condition numberλn/λm+1, which is expected to be substantially smaller than λn/λ1.
However, the method requires exact solution of systems withAV and for some problems m is not that small. Also, it is assumed that the projection Pb is computed exactly (or to a sufficient accuracy), which may be unfeasible in many applications.Method 2 (Defect-correction with projectors).
In the presence of round-off errors,x(0) may not be sufficiently accurate and b-Ax(0) may still contain components in the “bad" subspace. A defect-correction (iterative refinement) procedure may then help.  Let  x(0)=0 for    k=0,1,2,…, until convergence. Compute  r(k)=b-Ax(k).  Solve  Ad(k)=b-Ax(k) as follows:(i)
d(k,0)=VAV-1VTr(k),(ii)
Az(k)=(I-P)r(k), or   Ay(k)=r(k)-Ad(k,0),(iii)
d(k)=d(k,0)+y(k).Let  x(k+1)=x(k)+d(k).
In this method, it normally suffices with few defect-correction steps.
For some extremely ill-conditioned systems, the implementation of the defect-correction method as a preconditioning method may be necessary. Note then that as follows from Lemma6, (I-P)A is symmetric so the standard conjugate gradient method can be used.Method 3 (Preconditioning by a projection matrix).
Letx(0)=VAV-1VTb.  Solve   (I-P)Az=(I-P)b,  or  (I-P)Az=b-Ax(0) by CG iteration.  Let  x=z+x(0).
Note thatx(0) is contained in the null-space of (I-P)A, since (I-P)Ax(0)=(I-P)Pb=(P2-P)b=0. Here, the system (I-P)A is well-conditioned on the orthogonal complement to the null-space and, in addition, the right-hand-side has no or only small components in the bad subspace.
Methods1, 2, and 3 require accurate solution of systems with the matrix AV. It is a viable step for small values of m. However, when the dimension of the “bad" subspace of A is relatively big, it may be too costly. Furthermore, the iteration Method 3 involves two multiplications with A (one involved in P and one required to computeAz(k)) at each iteration step when computing the search direction vectors and is therefore particularly expensive.
Another issue to comment on is that the above methods are assumed to move the components of the eigenvectors for the smallest eigenvalues ofA to become exactly zero. However, this can be sensitive to perturbations and occurs only in exact arithmetic. As we have seen, it is a viable method for small dimensions of the subspace causing the ill-conditioning but it may be inefficient for larger dimensions.
Deflation methods have been used and analysed by [98–100], among others.
In the next section and Section4, we present a method which move the small eigenvalues to the cluster of bigger eigenvalues which is much less dependent on having the right subspace spanned by the columns of V and which do not require exact solution of systems with AV.


## 6.3. Augmented Matrix Preconditioning Methods, the Ideal Case

We now present an alternative method to handle ill-posed problems. In this method the small eigenvalues are moved to the cluster of bigger eigenvalues, instead of being deflated to zero, as in the deflation method. The method is an extension of the method presented in [101]. The presentation here is based in [25, 102]. First, we consider B=I+VVT as a (multiplicative) preconditioner to A.In this case one must scale the column vectors appearing inV properly. A method involving an automatic scaling is based on a projection matrix. Let then(233)B=I+σVAV-1VT,AV=VTAV.Let(λi,vi)i=1m be the eigenpairs of A for the smallest eigenvalues, 0<λ1≤λ2≤⋯λm. If V=[v1,…,vm], we get then(234)λ̃i=λi(BA)={λi+σ,1≤i≤m,λi,m+1≤i≤n.
Hence, σ determines how much the smallest eigenvalues are moved. If λm+1≤λ1+σ and λm+σ≤λn, then λm+1≤λ̃i≤λn, that is, the m smallest eigenvalues have been moved to the cluster [λm+1,λn] of bigger eigenvalues and the spectral condition number of BAis κ(BA)=λn/λm+1, which normally means a significant reduction, compared to κ(A)=λn/λ1.The above illustrates what can be achieved in an ideal case. In practice, the exact eigenvectors{v1,…,vm} (or the subspace spanned by them) are not known. Even if the eigenvectors are known it is not efficient to use them to form the matrix V because they are in general not sparse and the matrix vector multiplications with V will be costly. Hence, in practice other vectors which are sparse but spans about the same subspace as the smoothest eigenvectors must be used; otherwise the expense of the preconditioner would be too high. We consider therefore more general subspaces spanned by the column vectors of V. The next lemma will be useful.Lemma 7.
LetA be s.p.d. Then,
(235)P=A1/2V(VTAV)-1VTA1/2,
is an orthogonal projection, that is, P2=P and P*=P. Therefore, the only eigenvalues of P are 0 and 1.Proof.
Consider the following(236)P2=A1/2V(VTAV)-1VTA1/2A1/2V(VTAV)-1VTA1/2=A1/2V(VTAV)-1VTA1/2=P.The next theorem shows (what can also be expected) that the clustering can never get worse for expanding subspaces spanned by the column-vectors ofV, that is, there holds a monotonicity principle.Theorem 13 (see [103]).
LetA and Â be s.p.d. matrices of order n×n and let   Vk  be rectangular matrices of order n×mk,  k=1,2 such that rankVk=mk,  k=1,2. If ImV1⊆ImV2, then for all i,   1≤i≤n the following inequality holds
(237)λi((I+V2(V2TÂV2)-1V2T)A)≥λi((I+V1(V1TÂV1)-1V1T)A).Proof.
It is readily seen that the proposition holds ifF=V2(V2TÂV2)-1V2T-V1(V1TÂV1)-1V1T is negative definite. But since ImV1⊆ImV2, there exists some matrix Q of order m2×m1 such that V1=V2Q. Then, with DK=VkTÂVk, we have
(238)F=V2(D2-1-QD1-1QT)V2T=V2D2-1/2(I-D2-1/2QD1-1/2QTD2-1/2)D2-1/2V2T,
where
(239)P≡D2-1/2QD1-1/2QTD2-1/2=D2-1/2Q(QTD2Q)-1QTD2-1/2,
is an orthogonal projector (P2=P), whose eigenvalues are 0 and 1.Corollary 3.
IfImV1=ImV2 then I+V2D2-1V2T=I+V1D1-1V1T.Proof.
In this case,Q in V1=V2Q is invertible. Thus, D2-1/2Q(QTD2Q)-1QTD2-1/2=I.Remark 4.
The above corollary shows that the individual eigenvectors ofA are not needed when constructing the matrix V; we are rather interested in the subspace spanned by them.
The most interesting case for us is whenImV⊃span{v1,…,vm}, where vi are the eigenvalues of Afor the smallest eigenvalues λ1,…,λm. Then, the preconditioner B=I+σVAV-1VT moves the smallest eigenvalues λi at least to λi+σ, where σis a scaling parameter.Theorem 14.
LetB=I+σVAV-1VT and let  λ1,…,λm be the smallest eigenvalues of A. Then, for the eigenvalues of BA there holds. 									(240)λ̃i≥{λi+σ,1≤i≤mλi,m+1≤i≤n.Proof.
LetV1=[v1,…,vm], where {vi}1m are the first m eigenvalues of A and let V2=V. Then, Lemma 7 and (234) show the result.It may happen that the eigenvalues are moved too far so that the maximum eigenvalue ofBA is much larger that of A.Theorem 15.
LetA be s.p.d. of order n×n and let the rectangular matrix V of order n×m be defined as V=[v1,v2,…,vm]. Assume that rankV=m. Further, define Ã=(I+σVAV-1VT)A, where AV=VTAV. Then,
(241)λmax⁡(Ã)≤σ+λmax⁡(A).Proof.
The result follows from the following relations:(242)λmax⁡(Ã)≤λmax⁡(A)+σsup⁡xTV(VTAVV)-1VTxxTA-1x=λmax⁡(A)+σsup⁡xTA1/2V(VTAV)-1VTA1/2xxTx=λmax⁡(A)+σ,
where the last equality follows from Lemma 7It follows from Theorems 6.3 and16 that the optimal value of σ=λm+1, in which case κ(BA)≤(λn+σ)/λm+1. In general, λm+1 may not be known. With σ=λn, we obtain κ(BA)≤(2λn)/λm+1.The moral we can draw from the above is that the suggested technique can be a very useful means to reduce the condition number of a given matrixA if we have information about the eigenvectors corresponding to the smallest eigenvalues of A. Since in practice this is hardly ever the case, a natural step to undertake is to consider not the individual eigenvectors but the subspace spanned by some approximation of them.In the next section, we present a generalized form of the augmented matrix preconditioner which allows for both approximate subspaces and the replacement ofAV by a simpler matrix BV.


## 7. Preconditioners with an Approximate Subspace Correction Term

The preconditioner presented in the previous subsection will now be extended to include an approximate subspace correction term.We replace firstAV with a possibly simpler matrix BV. The resulting eigenvalue bounds are found in the next theorem.Theorem 16.
LetA be s.p.d. Define the preconditioner B as B=I+σVBV-1VT,where,σ>0,ImV⊇span{v1,…,vm}, wherevi are the eigenvectors of A for the smallest eigenvalues and BV is an m×m s.p.d. approximation of AV. Then, the eigenvalues λ(BA) of BA are bounded as follows:(a)
(243)min⁡{σλmin⁡(BV-1AV)+λ1,λm+1}≤λ(BA)≤σλmax⁡(BV-1AV)+λmax⁡(A).(b)
Withσ=λmax(A)/λmax(BV-1AV), we have(244)min⁡{λmax⁡(A)/κ(BV-1AV)+λ1,λm+1}≤λ(BA)≤2λmax⁡(A).Proof.
The minimal eigenvalue ofBA can be estimated as
(245)λmin⁡(BA)=inf⁡x⁡{xT(I+σVBV-1VT)xxTA-1x}=inf⁡x⁡{xTxxTA-1x+σxTVBV-1VTxxTVAV-1VTx⋅xTVAV-1VTxxTA-1x}=inf⁡x⁡{xTxxTA-1x+σλmin⁡(BV-1AV)xTVAV-1VTxxTA-1x}.
Here, infxxTx/xTA-1x=λ1 and infx;VTx=0xTx/xTA-1x=λm+1.
The lower bound equals the minimal eigenvalue of the matrixB̂A, where
(246)B̂=B̂(V)≡I+σ̂VAV-1VTσ̂=σλmin⁡(NV-1AV).
By Theorem14, we have with V1=span{v1,…,v+m} and V2 a matrix satisfying ImV2⊇ImV1, that
(247)λi(B̂(V2)A)≥λi(B̂(V1)A)≥min⁡{λ1+σ̂,λm+1}
and, in particular, for V2=V(248)λmin⁡(BA)≥min⁡(λ1+σ̂,λm+1)=min⁡{λ1+σλmin⁡(BV-1AV),λm+1}
which is the lower bound in part (a). The upper bound follows in a similar way. Since there is no upper bound assumed on rank V, and since supx,VTx=0xTx/xTA-1x≤λmax(A), we obtain
(249)λmax⁡(BA)≤σλmax⁡(BV-1AV)+λmax⁡(A).
If we let σ=λmax(A)/λmax(BV-1AV), we get the stated lower bound in (b) and λmax(BA)≤2λmax(A).Since normallyλmax(A) and λmax(BV-1AV) are readily estimated, the given choice of σ is viable. It may increase the maximal eigenvalue with a factor 2, which is acceptable.Corollary 4.
Ifκ(BV-1AV)≤λmax(A)λm+1, then
(250)κ(BA)≤2λmax⁡(A)λm+1,
that is, the upper bound, coincides with the bound where BV=AV. In particular, if κ(AV)≤λmax(A)/λm+1, then one can simply let BV=I or
(251)BV=diag⁡(AV).Hence, seen that the matrixAV in the preconditioner can be replaced with a simpler matrix BV where the action of BV-1 is cheap, without deteriorating the eigenvalue bounds.It still remains to weaken the assumption(252)Im⁡V⊇span⁡{v1,…vm},
as this is not easy to satisfy in practice. Due to space limitations this will not be presented here. Instead, refer to [56].


## 8. Krylov Subspace Methods for  Singular Systems

Singular systems, that is, with a nontrivial kernel, arise in various applications, such as in boundary value problems with pure Neumann type boundary conditions imposed on the whole boundary. Nullspaces of large dimension may arise in finite element methods using edge element methods, see, for example, [104] and in the analysis of Markov chains when stationary probability vectors of stochastic matrices are computed, see [105, 106], see also [107].A singular system does not always have a solution and it is more appropriate to consider the least squares problem: findx∈Rn such that ∥b-Ax∥≤∥b-Ax∥ for all x∈Rn. We recall that a basic iterative solution method to solve a linear system, either has the form(253)xk+1=xk+τkrk,rk=b-Axk,k=0,1,…,or the preconditioned form(254)solveBδk=τkrk,and update(255)xk+1=xk+δk,k=0,1,….Here,B is a preconditioner to A. Similarly, as we have seen, more involved methods, such as (generalized) CG-methods (GCG, GMRES, GCR, etc.) are based on approximations taken from a Krylov subspace(256)K(A,r0,k)={r0,Ar0,…,Akr0}orK(B-1A,r0,k).In general, they are based on a minimum residual approach where, at each iteration step, we compute an updated solution that satisfies the best approximation property(257)min⁡x∈K(A,r0,k)‖b-Ax‖ormin⁡x∈K(B-1A,r0,k)‖b-Ax‖.We will show that the convergence of such iterative solution methods can stall, or suffer a breakdown, when applied to certain singular systems. For the analysis, we will use the following properties of relevant subspaces for matrixA. In this generality, they can be stated for rectangular matrices.Definition 2.
LetA∈ℛm×n. Then R(A) of dimension ≤n, called the rangeof A, is the subspace spanned by the column vectors a.j that is
(258)y∈R(A)⟺y=∑j=1nαja.j.Definition 3.
𝒩(A), of dimension ≤n, is the nullspace of A, that is, the subspace of vectors v̲∈Rn s.t. Av̲=0̲.
By a classical result (see e.g., [6]), it holds R(A)⊥=𝒩(AT)Definition 4.
A linear systemAx=bis called consistent if b∈R(A) and inconsistentotherwise. If Ax=b is consistent, there exists a solution. Cle arly, any system with a nonsingular matrix A is consistent.
To provide a general method, applicable for all types of systems, we will use a least squares type of method, that is determine an approximate solution tox s.t. minx∥b-Ax∥ (or, similarly, a preconditioned form).
In practice, the approximations are mostly computed by a Krylov subspace method, where at each step a solutionxk is computed such that
(259)‖b-Axk‖≤‖b-Ax‖,∀x∈K(A,r0,k),
or a preconditioned form of the method. As the next Theorem shows even if the system is consistent, a breakdown can occur.Theorem 17.
Any minimum residual Krylov subspace method may suffer a breakdown for some initial vector if and only ifR(A)∩𝒩(A) contains a nontrivial vector.Proof.
The sufficiency follows since ifR(A)∩𝒩(A)≠{0̲}, there exists a nonzero vector x∈𝒩(A) which is also in R(A). Then, at some stage k, there exists a vector rk=y, where Ay=x, x∈𝒩(A). Hence, Ark=Ay=x, which implies A2rk=Ax=0, so a zero vector arises in the Krylov subspace at some stage. This means that convergence stalls. On the other hand, R(A)∩𝒩(A)={0̲} implies the existence of nonzero vectors Akr0 of any order in the Krylov subspace, which implies that there is an improved approximate solution for each higher stage k+1.Example 8.
Let(260)A=[01-10-12-101/21/2-100001].Here, Ae̲=0,e̲T=(1,1,,1,0), that is,  e̲∈𝒩(A).
Furthermore, there is a solution ofAy=e̲, for instance
(261)y=1/2[3310].
Hence, e̲∈R(A)∩𝒩(A), where e̲≠0̲. Since A2r0=0, with yas initial vector, the Krylov subspace for the system Ax=0 stalls at the second step.Corollary 5.
(a) If𝒩(A)=𝒩(AT), in particular if A is symmetric, then minimal residual type methods-based the Krylov subspace converges.
(b) More generally, this holds ifA is a (H-) normal matrix.Proof.
(a) SinceR(A)⊥=𝒩(AT), it holds
(262)R(A)⊥=N(A)and,henceR(A)∩N(A)=R(A)∩R(A)⊥={0̲}.
(b) For a normal matrix, there exists a unitary matrixU that diagonalizes A, that is
(263)UTAU=D.
Hence,
(264)UTATU=DT=D.
Therefore, if AUx=Dx=0 for some x≠0, then also ATUx=Dx=0, so
(265)Ux∈N(A),Ux∈N(AT),
for any such vector x. Since U is nonsingular, this implies 𝒩(A)=𝒩(AT).Remark 5.
Corollary5 can be extended to H-normal matrices, that is matrices for which A commutes with its H-adjoint,
(266)A′=H-1A*H,
for some Hermitian matrix H see, for example, [6].
A remedy to avoid breakdowns for matricesA for which the vector space 𝒱=R(A)∩𝒩(A) is nontrivial, is to work in a subspace orthogonal to 𝒱. This can be achieved by use of the augmented subspace projection method in Section 6. This method works also to avoid situations, where R(A)∩𝒩(A) contains eigenvectors to A corresponding to nearly zero eigenvalues, causing a near breakdown or, in finite precision computations, an actual breakdown. For further comments on near breakdowns, see, for example, [83–85].


## 9. Concluding Remarks

Some milestones in the development iterative solutions methods have been presented. By the combination of improved methods and the developments of computer hardware one can presently solve problems with a degree of freeadoms nearly billionfold compared to that in the early ages of the computer age.There remains still, however, very difficult problems such as in multiphysics and heterogeneous media problems and various forms of inverse problems, which need further improvement of solution methods.Some problems, such as those arising in constrained optimization and mixed finite element methods, lead to matrices on saddle pointform. Due to space limitations, they have not been discussed in this paper, however, see, for example, [108]. In the last centuries, much work has been devoted to multigrid, algebraic multigrid and multilevel iteration methods which have shown an optimal order of performance for many types of problems, for example see [58]. Also, domain decomposition methods which go back to the Schwarz alternating decomposition method, have shown developments, see, for example, [109–112]. For an early survey of domain decomposition methods, see [113]. For the same reason, they could not be discussed in this paper. Much work has also been developed to parallelization aspects of solution methods. This topic deserves a separate survey article and has also not been discussed in this paper.



---

*Source: 972794-2010-12-09.xml*