<< . .

. 6
( : 26)



. . >>

and accuracy a¬ected by changes in a(x, y)? Try a(x, y) = .1 + x and
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

<< . .

. 6
( : 26)



. . >>