examine the accuracy of the result. Explain your ¬ndings. Compare

the execution times on your computing environment (using the cputime

command in MATLAB, for instance).

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

31

CONJUGATE GRADIENT ITERATION

2.8.12. Use the Jacobi and symmetric Gauss“Seidel iterations from Chapter 1

to solve the elliptic boundary value problem in § 2.7. How does the

performance compare to CG and PCG?

2.8.13. Implement Jacobi (1.17) and symmetric Gauss“Seidel (1.18) precondi-

tioners for the elliptic boundary value problem in § 2.7. Compare the

performance with respect to both computer time and number of itera-

tions to preconditioning with the Poisson solver.

2.8.14. Modify pcgsol so that φ(x) is computed and stored at each iterate

and returned on output. Plot φ(xn ) as a function of n for each of the

examples.

2.8.15. Apply CG and PCG to solve the ¬ve-point discretization of

’uxx (x, y) ’ uyy (x, y) + ex+y u(x, y) = 1, 0 < x, y, < 1,

subject to the inhomogeneous Dirichlet boundary conditions

u(x, 0) = u(x, 1) = u(1, y) = 0, u(0, y) = 1, 0 < x, y < 1.

Experiment with di¬erent mesh sizes and preconditioners (Fast Poisson

solver, Jacobi, and symmetric Gauss“Seidel).

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

32 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

Chapter 3

GMRES Iteration

3.1. The minimization property and its consequences

The GMRES (Generalized Minimum RESidual) was proposed in 1986 in [167]

as a Krylov subspace method for nonsymmetric systems. Unlike CGNR,

GMRES does not require computation of the action of AT on a vector. This is

a signi¬cant advantage in many cases. The use of residual polynomials is made

more complicated because we cannot use the spectral theorem to decompose

A. Moreover, one must store a basis for Kk , and therefore storage requirements

increase as the iteration progresses.

The kth (k ≥ 1) iteration of GMRES is the solution to the least squares

problem

(3.1) minimizex∈x0 +Kk b ’ Ax 2 .

The beginning of this section is much like the analysis for CG. Note that

if x ∈ x0 + Kk then

k’1

γj Aj r0

x = x0 +

j=0

and so

k’1 k

j+1

γj’1 Aj r0 .

b ’ Ax = b ’ Ax0 ’ γj A r0 = r0 ’

j=0 j=1

Hence if x ∈ x0 + Kk then r = p(A)r0 where p ∈ Pk is a residual polynomial.

¯ ¯

We have just proved the following result.

Theorem 3.1.1. Let A be nonsingular and let xk be the kth GMRES

iteration. Then for all pk ∈ Pk

¯

(3.2) rk = min p(A)r0

¯ ¤ pk (A)r0 2 .

¯

2 2

p∈Pk

From this we have the following corollary.

Corollary 3.1.1. Let A be nonsingular and let xk be the kth GMRES

iteration. Then for all pk ∈ Pk

¯

rk 2

(3.3) ¤ pk (A) 2 .

¯

r0 2

33

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

34 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

We can apply the corollary to prove ¬nite termination of the GMRES

iteration.

Theorem 3.1.2. Let A be nonsingular. Then the GMRES algorithm will

¬nd the solution within N iterations.

Proof. The characteristic polynomial of A is p(z) = det(A ’ zI). p has

degree N , p(0) = det(A) = 0 since A is nonsingular, and so

pN (z) = p(z)/p(0) ∈ PN

¯

is a residual polynomial. It is well known [141] that p(A) = pN (A) = 0. By

¯

(3.3), rN = b ’ AxN = 0 and hence xN is the solution.

In Chapter 2 we applied the spectral theorem to obtain more precise infor-

mation on convergence rates. This is not an option for general nonsymmetric

matrices. However, if A is diagonalizable we may use (3.2) to get information

from clustering of the spectrum just like we did for CG. We pay a price in that

we must use complex arithmetic for the only time in this book. Recall that

A is diagonalizable if there is a nonsingular (possibly complex!) matrix V such

that

A = V ΛV ’1 .

Here Λ is a (possibly complex!) diagonal matrix with the eigenvalues of A on

the diagonal. If A is a diagonalizable matrix and p is a polynomial then

p(A) = V p(Λ)V ’1

A is normal if the diagonalizing transformation V is orthogonal. In that case

the columns of V are the eigenvectors of A and V ’1 = V H . Here V H is the

complex conjugate transpose of V . In the remainder of this section we must

use complex arithmetic to analyze the convergence. Hence we will switch to

complex matrices and vectors. Recall that the scalar product in C N , the space

of complex N -vectors, is xH y. In particular, we will use the l2 norm in C N .

Our use of complex arithmetic will be implicit for the most part and is needed

only so that we may admit the possibility of complex eigenvalues of A.

We can use the structure of a diagonalizable matrix to prove the following

result.

Theorem 3.1.3. Let A = V ΛV ’1 be a nonsingular diagonalizable matrix.

Let xk be the kth GMRES iterate. Then for all pk ∈ Pk

¯

rk 2

(3.4) ¤ κ2 (V ) max |¯k (z)|.

p

r0 z∈σ(A)

2

Proof. Let pk ∈ Pk . We can easily estimate pk (A)

¯ ¯ by

2

V ’1

pk (A)

¯ ¤V pk (Λ)

¯ ¤ κ2 (V ) max |¯k (z)|,

p

2 2 2 2

z∈σ(A)

as asserted.

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

35

GMRES ITERATION

It is not clear how one should estimate the condition number of the

diagonalizing transformation if it exists. If A is normal, of course, κ2 (V ) = 1.

As we did for CG, we look at some easy consequences of (3.3) and (3.4).

Theorem 3.1.4. Let A be a nonsingular diagonalizable matrix. Assume

that A has only k distinct eigenvalues. Then GMRES will terminate in at most

k iterations.

Theorem 3.1.5. Let A be a nonsingular normal matrix. Let b be a linear

combination of k of the eigenvectors of A

k

b= γl uil .

l=1

Then the GMRES iteration, with x0 = 0, for Ax = b will terminate in at most

k iterations.

3.2. Termination

As is the case with CG, GMRES is best thought of as an iterative method.

The convergence rate estimates for the diagonalizable case will involve κ2 (V ),

but will otherwise resemble those for CG. If A is not diagonalizable, rate

estimates have been derived in [139], [134], [192], [33], and [34]. As the set of

nondiagonalizable matrices has measure zero in the space of N — N matrices,

the chances are very high that a computed matrix will be diagonalizable. This

is particularly so for the ¬nite di¬erence Jacobian matrices we consider in

Chapters 6 and 8. Hence we con¬ne our attention to diagonalizable matrices.

As was the case with CG, we terminate the iteration when

(3.5) rk ¤· b

2 2

for the purposes of this example. We can use (3.3) and (3.4) directly to estimate

the ¬rst k such that (3.5) holds without requiring a lemma like Lemma 2.3.2.

Again we look at examples. Assume that A = V ΛV ’1 is diagonalizable,

that the eigenvalues of A lie in the interval (9, 11), and that κ2 (V ) = 100.

We assume that x0 = 0 and hence r0 = b. Using the residual polynomial

pk (z) = (10 ’ z)k /10k we ¬nd

¯

rk 2

¤ (100)10’k = 102’k .

r0 2

Hence (3.5) holds when 102’k < · or when

k > 2 + log10 (·).

¤ ρ < 1. Let pk (z) = (1 ’ z)k . It is a direct

Assume that I ’ A ¯

2

consequence of (3.2) that

¤ ρk r 0 2 .

(3.6) rk 2

The estimate (3.6) illustrates the potential bene¬ts of a good approximate

inverse preconditioner.

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

36 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

The convergence estimates for GMRES in the nonnormal case are much

less satisfying that those for CG, CGNR, CGNE, or GMRES in the normal

case. This is a very active area of research and we refer to [134], [33], [120],

[34], and [36] for discussions of and pointers to additional references to several

questions related to nonnormal matrices.

3.3. Preconditioning

Preconditioning for GMRES and other methods for nonsymmetric problems is

di¬erent from that for CG. There is no concern that the preconditioned system

be spd and hence (3.6) essentially tells the whole story. However there are two

di¬erent ways to view preconditioning. If one can ¬nd M such that

I ’ MA = ρ < 1,

2

then applying GMRES to M Ax = M b allows one to apply (3.6) to the

preconditioned system. Preconditioning done in this way is called left

preconditioning. If r = M Ax ’ M b is the residual for the preconditioned

system, we have, if the product M A can be formed without error,

ek rk

2 2

¤ κ2 (M A) ,

e0 r0

2 2

by Lemma 1.1.1. Hence, if M A has a smaller condition number than A, we

might expect the relative residual of the preconditioned system to be a better

indicator of the relative error than the relative residual of the original system.

If

I ’ AM 2 = ρ < 1,

one can solve the system AM y = b with GMRES and then set x = M y. This is

called right preconditioning. The residual of the preconditioned problem is the

same as that of the unpreconditioned problem. Hence, the value of the relative

residuals as estimators of the relative error is unchanged. Right preconditioning

has been used as the basis for a method that changes the preconditioner as the

iteration progresses [166].

One important aspect of implementation is that, unlike PCG, one can

apply the algorithm directly to the system M Ax = M b (or AM y = b). Hence,

one can write a single matrix-vector product routine for M A (or AM ) that

includes both the application of A to a vector and that of the preconditioner.

Most of the preconditioning ideas mentioned in § 2.5 are useful for GMRES

as well. In the examples in § 3.7 we use the Poisson solver preconditioner for

nonsymmetric partial di¬erential equations. Multigrid [99] and alternating

direction [8], [182] methods have similar performance and may be more

generally applicable. Incomplete factorization (LU in this case) preconditioners

can be used [165] as can polynomial preconditioners. Some hybrid algorithms

use the GMRES/Arnoldi process itself to construct polynomial preconditioners

for GMRES or for Richardson iteration [135], [72], [164], [183]. Again we

mention [8] and [12] as a good general references for preconditioning.

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

37

GMRES ITERATION

3.4. GMRES implementation: Basic ideas

Recall that the least squares problem de¬ning the kth GMRES iterate is

minimizex∈x0 +Kk b ’ Ax 2 .

Suppose one had an orthogonal projector Vk onto Kk . Then any z ∈ Kk can

be written as

k

yl vlk ,

z=

l=1

where vlk is the lth column of Vk . Hence we can convert (3.1) to a least squares

problem in Rk for the coe¬cient vector y of z = x ’ x0 . Since

x ’ x0 = Vk y

for some y ∈ Rk , we must have xk = x0 + Vk y where y minimizes

b ’ A(x0 + Vk y) = r0 ’ AVk y 2 .

2

Hence, our least squares problem in Rk is

(3.7) minimizey∈Rk r0 ’ AVk y 2 .

This is a standard linear least squares problem that could be solved by a QR

factorization, say. The problem with such a direct approach is that the matrix

vector product of A with the basis matrix Vk must be taken at each iteration.

If one uses Gram“Schmidt orthogonalization, however, one can represent

(3.7) very e¬ciently and the resulting least squares problem requires no extra

multiplications of A with vectors. The Gram“Schmidt procedure for formation

of an orthonormal basis for Kk is called the Arnoldi [4] process. The data are

vectors x0 and b, a map that computes the action of A on a vector, and a

dimension k. The algorithm computes an orthonormal basis for Kk and stores

it in the columns of V .

Algorithm 3.4.1. arnoldi(x0 , b, A, k, V )

1. De¬ne r0 = b ’ Ax0 and v1 = r0 / r0 2 .

2. For i = 1, . . . , k ’ 1

i T

Avi ’ j=1 ((Avi ) vj )vj

vi+1 = i T

Avi ’ j=1 ((Avi ) vj )vj 2

If there is never a division by zero in step 2 of Algorithm arnoldi, then

the columns of the matrix Vk are an orthonormal basis for Kk . A division by

zero is referred to as breakdown and happens only if the solution to Ax = b is

in x0 + Kk’1 .

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.

Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.

38 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Lemma 3.4.1. Let A be nonsingular, let the vectors vj be generated by

Algorithm arnoldi, and let i be the smallest integer for which

i

((Avi )T vj )vj = 0.

Avi ’

j=1

Then x = A’1 b ∈ x0 + Ki .

Proof. By hypothesis Avi ∈ Ki and hence AKi ‚ Ki . Since the columns of

Vi are an orthonormal basis for Ki , we have

AVi = Vi H,

where H is an i — i matrix. H is nonsingular since A is. Setting β = r0 2

and e1 = (1, 0, . . . , 0)T ∈ Ri we have

ri = b ’ Axi = r0 ’ A(xi ’ x0 ) 2 .

2 2

Since xi ’ x0 ∈ Ki there is y ∈ Ri such that xi ’ x0 = Vi y. Since r0 = βVi e1

and Vi is an orthogonal matrix

ri = Vi (βe1 ’ Hy) = βe1 ’ Hy Ri+1 ,

2 2

where · Rk+1 denotes the Euclidean norm in Rk+1 .

Setting y = βH ’1 e1 proves that ri = 0 by the minimization property.

The upper Hessenberg structure can be exploited to make the solution of

the least squares problems very e¬cient [167].

If the Arnoldi process does not break down, we can use it to implement

GMRES in an e¬cient way. Set hji = (Avj )T vi . By the Gram“Schmidt

construction, the k +1—k matrix Hk whose entries are hij is upper Hessenberg,

i.e., hij = 0 if i > j + 1. The Arnoldi process (unless it terminates prematurely