<<

. 9
( 16)



>>

¬’1 0 0 4 ’1 0 ’1 0 0 ·
¬ ·
¬ 0 ’1 0 ’1 4 ’1 0 ’1 0 ·
¬ ·
¬ 0 0 ’1 0 ’1 4 0 0 ’1 ·
¬ ·
¬ 0 0 0 ’1 0 0 4 ’1 0 ·
¬ ·
 0 0 0 0 ’1 0 ’1 4 ’1 
0 0 0 0 0 ’1 0 ’1 4


m = 3 — 3 : rowwise ordering.


« 
’1 ’1 0 0
4 0 0 0 0
¬ ·
’1 0 ’1 0
0 4 0 0 0
¬ ·
¬ ·
’1 ’1 ’1 ’1
0 0 4 0 0
¬ ·
¬ ·
0 ’1 0 ’1
0 0 0 4 0
¬ ·
¬ ·
0 0 ’1 ’1
0 0 0 0 4
¬ ·
¬ ·
¬’1 ’1 ’1 0 0 0·
¬ ·
4 0 0
¬’1 0 ’1 ’1 0 0·
¬ ·
0 4 0
 0 ’1 ’1 0 ’1 0
0 0 4
0 0 ’1 ’1 ’1 0 0 0 4


red-black ordering:
red: node 1, 3, 5, 7, 9 from rowwise ordering
black: node 2, 4, 6, 8 from rowwise ordering

Figure 5.2. Comparison of orderings.



that contribute to the corresponding row of the discretization matrix. In
the example of the ¬ve-point stencil, starting with rowwise numbering, one
can combine all odd indices to a block S1 (the “red nodes”) and all even
indices to a block S2 (the “black” nodes). Here we have p = 2. We call this
a red-black ordering (see Figure 5.2). If two “colours” are not su¬cient, one
can choose p > 2.
We return to the SOR method and its convergence: In the following the
iteration matrix will be denoted by MSOR(ω) with the relaxation parameter
ω. Likewise, MJ and MGS are the iteration matrices of Jacobi™s and the
Gauss“Seidel method, respectively. General propositions are summarized
in the following theorem:
5.1. Linear Stationary Iterative Methods 213

Theorem 5.6 (of Kahan; Ostrowski and Reich)


MSOR(ω) ≥ |1 ’ ω| for ω = 0.
(1)
(2) If A is symmetric and positive de¬nite, then
for ω ∈ (0, 2) .
MSOR(ω) < 1

2
Proof: See [16, pp. 91 f.].

Therefore, we use only ω ∈ (0, 2). For a useful procedure we need more
information about the optimal relaxation parameter ωopt , given by
MSOR(ωopt ) = min MSOR(ω) ,
0<ω<2

and about the size of the contraction number. This is possible only if the
ordering of equations and unknowns has certain properties:
De¬nition 5.7 A matrix A ∈ Rm,m is consistently ordered if for the
partition (5.18), D is nonsingular and
C(±) := ±’1 D’1 L + ±D’1 R
has eigenvalues independent of ± for ± ∈ C\{0}.
There is a connection to the possibility of a multi-colour ordering, because
a matrix in the block form (5.39) is consistently ordered if it is block-
tridiagonal (i.e., Aij = 0 for |i ’ j| > 1) and the diagonal blocks Aii are
nonsingular diagonal matrices (see [28, pp. 114 f.]).
In the case of a consistently ordered matrix one can prove a relation
between the eigenvalues of MJ , MGS , and MSOR(ω) . From this we can see
how much faster the Gauss“Seidel method converges than Jacobi™s method:
Theorem 5.8 If A is consistently ordered, then
(MJ )2 = (MGS ) .

2
Proof: For a special case see Remark 5.5.2 in [16].

Due to (5.4) we can expect a halving of the number of iteration steps,
but this does not change the asymptotic statement (5.27).
Finally, in the case that Jacobi™s method converges the following theorem
holds:
Theorem 5.9 Let A be consistently ordered with nonsingular diagonal ma-
trix D, the eigenvalues of MJ being real and β := (MJ ) < 1. Then we have
for the SOR method:
2
(1) ωopt = ,
1 + (1 ’ β 2 )1/2
214 5. Iterative Methods for Systems of Linear Equations
±
 2 2 1/2

 1 ’ ω + 1 ω 2 β 2 + ωβ 1 ’ ω + ω β

 2 4
(2) (MSOR(ω) ) = for 0 < ω < ωopt





ω’1 for ωopt ¤ ω < 2 ,

β2
(3) MSOR(ωopt ) = .
(1 + (1 ’ β 2 )1/2 )2

Proof: See [18, p. 216].
2

ρ( M SOR(ω) )
1




0 ω opt 2ω
1

Figure 5.3. Dependence of MSOR(ω) on ω.


If (MJ ) is known for Jacobi™s method, then ωopt can be calculated. This
is the case in the example of the ¬ve-point stencil discretization on a square:
From (5.26) and Theorem 5.9 it follows that
π2
2
π
+ O(n’4 ) ;
=1’
(MGS ) = cos 2
n n
hence
π
ωopt = 2/ 1 + sin n ,

ωopt ’ 1 = 1 ’ 2 π + O(n’2 ) .
MSOR(ωopt ) = n
Therefore, the optimal SOR method has a lower complexity than all
methods described up to now.
Correspondingly, the number of operations to reach the relative er-
ror level µ > 0 is reduced to ln 1 O(m3/2 ) operations in comparison to
µ
ln 1 O(m2 ) operations for the previous procedures.
µ
Table 5.1 gives an impression of the convergence for the model problem.
It displays the theoretically to be expected values for the numbers of iter-
ations of the Gauss“Seidel method (mGS ), as well as for the SOR method
5.1. Linear Stationary Iterative Methods 215

n mGS mSOR
8 43 8
16 178 17
32 715 35
64 2865 70
128 11466 140
256 45867 281

Table 5.1. Gauss“Seidel and optimal SOR method for the model problem.

with optimal relaxation parameter (mSOR ). Here we use the very moderate
termination criterion µ = 10’3 measured in the Euclidean norm.
The optimal SOR method is superior, even if we take into account the
almost doubled e¬ort per iteration step. But generally, ωopt is not known
explicitly. Figure 5.3 shows that it is probably better to overestimate ωopt
instead of underestimating. More generally, one can try to improve the
relaxation parameter during the iteration:
If (MJ ) is a simple eigenvalue, then this also holds true for the spectral
radius (MSOR(ω) ). The spectral radius can thus be approximated by the
power method on the basis of the iterates. By Theorem 5.9 (3) one can
approximate (MJ ), and by Theorem 5.9 (1) then also ωopt .
This basic principle can be extended to an algorithm (see, for example,
[18, Section 9.5]), but the upcoming overall procedure is no longer a linear
stationary method.


5.1.5 Extrapolation Methods
Another possibility for an extension of the linear stationary methods, re-
lated to the adaption of the relaxation parameter, is the following: Starting
with a linear stationary basic iteration xk+1 := ¦ xk we de¬ne a new
˜ ˜
iteration by
x(k+1) := ωk ¦ x(k) + (1 ’ ωk )x(k) , (5.41)
with extrapolation factors ωk to be chosen. A generalization of this de¬-
nition is to start with the iterates of the basic iteration x(0) , x(1) , . . .. The
˜ ˜
iterates of the new method are to be determined by
k
(k)
±kj x(j) ,
x := ˜
j=0

with ±kj de¬ned by a polynomial pk ∈ Pk , with the property pk (t) =
k j
j=0 ±kj t and pk (1) = 1. For an appropriate de¬nition of such extrapola-
tion or semi-iterative methods we need to know the spectrum of the basic
iteration matrix M , since the error e(k) = x(k) ’ x satis¬es
e(k) = pk (M )e(0) ,
216 5. Iterative Methods for Systems of Linear Equations

where M is the iteration matrix of the basic iteration. This matrix should
be normal, for example, such that
pk (M ) = (pk (M ))
2

holds. Then we have the obvious estimation
¤ pk (M )e(0) ¤ pk (M ) ¤ (pk (M )) e(0) 2 .
e(k) e(0) (5.42)
2 2 2 2

If the method is to be de¬ned in such a way that
(pk (M )) = max |pk (»)| » ∈ σ(M )
is minimized by choosing pk , then the knowledge of the spectrum σ(M ) is
necessary. Generally, instead of this, we assume that suitable supersets are
known: If σ(M ) is real and
a¤»¤b for all » ∈ σ(M ) ,
then, due to
¤ max pk (») e(0) 2 ,
e(k) 2 »∈[a,b]

it makes sense to determine the polynomials pk as a solution of the
minimization problem on [a, b],
max |pk (»)| ’ min p ∈ Pk
for all with p(1) = 1 . (5.43)
»∈[a,b]

In the following sections we will introduce methods with an analogous
convergence behaviour, without control parameters necessary for their
de¬nition.
For further information on semi-iterative methods see, for example, [16,
Chapter 7].



Exercises
5.1 Investigate Jacobi™s method and the Gauss“Seidel method for solving
the linear system of equations Ax = b with respect to their convergence if
we have the following system matrices:
«  « 
1 2 ’2 2 ’1 1
1
(a) A =  1 1 1  , (b) A =  2 2 2  .
2
’1 ’1 2
221

5.2 Prove the consistency of the SOR method.

5.3 Prove Theorem 5.6, (1).
5.2. Gradient and Conjugate Gradient Methods 217

5.2 Gradient and Conjugate Gradient Methods
In this section let A ∈ Rm,m be symmetric and positive de¬nite. Then the
system of equations Ax = b is equivalent to the problem
1
Minimize f (x) := xT Ax ’ bT x for x ∈ Rm , (5.44)
2
since for such a functional the minima and stationary points coincide, where
a stationary point is an x satisfying
0 = ∇f (x) = Ax ’ b . (5.45)
In contrast to the notation x · y for the “short” space vectors x, y ∈ Rd we
write here the Euclidean scalar product as matrix product xT y.
For the ¬nite element discretization this corresponds to the equivalence
of the Galerkin method (2.23) with the Ritz method (2.24) if A is the
sti¬ness matrix and b the load vector (see (2.34) and (2.35)). More generally,
Lemma 2.3 implies the equivalence of (5.44) and (5.45), if as bilinear form
the so-called energy scalar product
:= xT Ay
x, y (5.46)
A

is chosen.
A general iterative method to solve (5.44) has the following structure:
De¬ne a search direction d(k) .
˜
± ’ f (±) := f x(k) + ±d(k)
Minimize (5.47)
exactly or approximately, with the solution ±k .
x(k+1) := x(k) + ±k d(k) .
De¬ne (5.48)
If f is de¬ned as in (5.44), the exact ±k can be computed from the condition
˜
f (±) = 0 and
T (k)
˜
f (±) = ∇f x(k) + ±d(k) d
as
T
g (k) d(k)
±k = ’ , (5.49)
T
d(k) Ad(k)
where
g (k) := Ax(k) ’ b = ∇f x(k) . (5.50)
The error of the kth iterate is denoted by e(k) :
e(k) := x(k) ’ x .
Some relations that are valid in this general fromework are the following:
Due to the one-dimensional minimization of f , we have
T
g (k+1) d(k) = 0 , (5.51)
218 5. Iterative Methods for Systems of Linear Equations

and from (5.50) we can conclude immediately that
Ae(k) = g (k) , e(k+1) = e(k) + ±k d(k) , (5.52)
(k+1) (k) (k)
g =g + ±k Ad . (5.53)
We consider the energy norm
1/2
:= xT Ax
x (5.54)
A

induced by the energy scalar product. For a ¬nite element sti¬ness matrix
A with a bilinear form a we have the correspondence
= a(u, u)1/2 = u
x A a
m
for u = i=1 xi •i if the •i are the underlying basis functions. Comparing
the solution x = A’1 b with an arbitrary y ∈ Rm leads to
1
f (y) = f (x) + y ’ x 2 , (5.55)
A
2
so that condition (5.44) also minimizes the distance to x in · A . The
energy norm will therefore have a special importance. Measured in the
energy norm we have, due to (5.52),
T T
2
= e(k) g (k) = g (k) A’1 g (k) ,
e(k) A
and therefore due to (5.52) and (5.51),
T
2
e(k+1) = g (k+1) e(k) .
A

The vector ’∇f x(k) in x(k) points in the direction of the locally steepest
descent, which motivates the gradient method, i.e.,
d(k) := ’g (k) , (5.56)
and thus
T
d(k) d(k)
±k = . (5.57)
T
d(k) Ad(k)
The above identities imply for the gradient method
T
d(k) d(k)
(k+1) 2 (k) T (k)
1 ’ ±k
(k)
e(k) 2
e =g + ±k Ad e = A T
d(k) A’1 d(k)
and thus by means of the de¬nition of ±k from (5.57)
± 
T (k) 2
 
 
d(k) d
2 2
’x A = x ’x A 1’
(k+1) (k)
x .
 (k) T Ad(k) d(k) T A’1 d(k) 
 
d

With the inequality of Kantorovich (see, for example, [28, p. 132]),
xT Ax xT A’1 x
2
1 1/2 1 ’1/2
¤ κ +κ ,
2 2 2
(xT x)
5.2. Gradient and Conjugate Gradient Methods 219

where κ := κ(A) is the spectral condition number, and the relation

(a ’ 1)2
4
1’ = for a > 0 ,
2 2
a1/2 + a’1/2 (a + 1)

we obtain the following theorem:

Theorem 5.10 For the gradient method we have
k
κ’1
’x ¤ x(0) ’ x
(k)
x . (5.58)
A A
κ+1
This is the same estimate as for the optimally relaxed Richardson method
(with the sharper estimate M A ¤ κ’1 instead of (M ) ¤ κ’1 ). The
κ+1 κ+1
essential di¬erence lies in the fact that this is possible without knowledge
of the spectrum of A.
Nevertheless, for ¬nite element discretizations we have the same poor
convergence rate as for Jacobi™s or similar methods. The reason for this
T
de¬ciency lies in the fact that due to (5.51), we have g (k+1) g (k) = 0, but
T
in general not g (k+2) g (k) = 0. On the contrary, these search directions are
very often almost parallel, as can be seen from Figure 5.4.


.
m = 2:
x (0)

f = constant
(contour lines)


Figure 5.4. Zigzag behaviour of the gradient method.


The reason for this problem is the fact that for large κ the search di-
rections g (k) and g (k+1) can be almost parallel with respect to the scalar
products ·, · A (see Exercise 5.4), but with respect to · A the distance
to the solution will be minimized (see (5.55)).
The search directions d(k) should be orthogonal with respect to ·, · A ,
which we call conjugate.

De¬nition 5.11 Vectors d(0) , . . . , d(l) ∈ Rm are conjugate if they satisfy

d(i) , d(j) =0 for i, j = 0, . . . , l , i = j .
A

If the search directions of a method de¬ned according to (5.48), (5.49) are
chosen as conjugate, it is called a method of conjugate directions.
Let d(0) , . . . , d(m’1) be conjugate directions. Then they are also linearly
independent and thus form a basis in which the solution x of (5.1) can be
220 5. Iterative Methods for Systems of Linear Equations

represented, say by the coe¬cients γk :
m’1
γk d(k) .
x=
k=0

Since the d(k) are conjugate and Ax = b holds, we have
T
d(k) b
γk = , (5.59)
T
d(k) Ad(k)
and the γk can be calculated without knowledge of x. If the d(k) would by
given a priori, for example by orthogonalization of a basis with respect to
·, · A , then x would be determined by (5.59).
If we apply (5.59) to determine the coe¬cients for x ’ x(0) in the form
m’1
x’x (0)
γk d(k) ,
=
k=0

which means replacing b with b ’ Ax(0) in (5.59), then we get
T
g (0) d(k)
γk = ’ .
T
d(k) Ad(k)
For the kth iterate we have, according to (5.48);
k’1
(k) (0)
±i d(i)
x =x +
i=0

and therefore (see (5.50))
k’1
(k) (0)
±i Ad(i) .
g =g +
i=0

For a method of conjugate directions this implies
T T
g (k) d(k) = g (0) d(k)
and therefore
T
g (k) d(k)
γk = ’ = ±k ,
T
d(k) Ad(k)
which means that x = x(m) . A method of conjugate directions therefore
is exact after at most m steps. Under certain conditions such a method
may terminate before reaching this step number with g (k) = 0 and the
¬nal iterate x(k) = x. If m is very large, this exactness of a method of
conjugate directions is less important than the fact that the iterates can
be interpreted as the solution of a minimization problem approximating
(5.44):
5.2. Gradient and Conjugate Gradient Methods 221

Theorem 5.12 The iterates x(k) that are determined by a method of con-
jugate directions minimize the functional f from (5.44) as well as the error
x(k) ’ x A on x(0) + Kk (A; g (0) ), where

Kk (A; g (0) ) := span d(0) , . . . , d(k’1) .
This is due to
T
for i = 0, . . . , k ’ 1 .
g (k) d(i) = 0 (5.60)

Proof: It is su¬cent to prove (5.60). Due to the one-dimensional mini-
mization this holds for k = 1 and for i = k ’ 1 (see (5.51) applied to k ’ 1).
To conclude the assertion for k from its knowledge for k ’ 1, we note that
(5.53) implies, for 0 ¤ i < k ’ 1,
T T
g (k) ’ g (k’1) = ±k’1 d(i) Ad(k’1) = 0 .
d(i)
2

In the method of conjugate gradients, or CG method, the d(k) are
determined during the iteration by the ansatz
d(k+1) := ’g (k+1) + βk d(k) . (5.61)
Then we have to clarify whether
d(k) , d(i) = 0 for k > i
A

can be obtained. The necessary requirement d(k+1) , d(k) = 0 leads to
A

’ g (k+1) , d(k) ⇐’
+ βk d(k) , d(k) =0
A A
T
g (k+1) Ad(k)
βk = . (5.62)
T
d(k) Ad(k)
In applying the method it is recommended not to calculate g (k+1) directly
but to use (5.53) instead, because Ad(k) is already necessary to determine
±k and βk .
The following equivalences hold:
Theorem 5.13 In case the CG method does not terminate prematurely
with x(k’1) being the solution of (5.1), then we have for 1 ¤ k ¤ m
Kk (A; g (0) ) = span g (0) , Ag (0) , . . . , Ak’1 g (0)
(5.63)
= span g (0) , . . . , g (k’1) .
Furthermore,
T
0 for i = 0, . . . , k ’ 1, and
g (k) g (i) = (5.64)
dim Kk (A; g (0) ) = k.
222 5. Iterative Methods for Systems of Linear Equations

The space Kk (A; g (0) ) = span g (0) , Ag (0) , . . . , Ak’1 g (0) is called the
Krylov (sub)space of dimension k of A with respect to g (0) .

Proof: The identities (5.64) are immediate consequences of (5.63) and
Theorem 5.12. The proof of (5.63) is given by induction:
For k = 1 the assertion is trivial. Let us assume that for k ≥ 1 the identity
(5.63) holds and therefore also (5.64) does. Due to (5.53) (applied to k ’ 1)
it follows that
g (k) ∈ A Kk A; g (0) ‚ span g (0) , . . . , Ak g (0)
and thus
span g (0) , . . . , g (k) = span g (0) , . . . , Ak g (0) ,
because the left space is contained in the right one and the dimension of
the left subspace is maximal (= k + 1) due to (5.64) and g (i) = 0 for all
i = 0, . . . , k. The identity
span d(0) , . . . , d(k) = span g (0) , . . . , g (k)
2
follows from the induction hypothesis and (5.61).

The number of operations per iteration can be reduced to one matrix
vector, two scalar products, and three SAXPY operations, if the following
equivalent terms are used:
T T
g (k) g (k) g (k+1) g (k+1)
±k = , βk = . (5.65)
T T
d(k) Ad(k) g (k) g (k)
Here a SAXPY operation is of the form
z := x + ±y
for vectors x, y, z and a scalar ±.
The identities (5.65) can be seen as follows: Concerning ±k we note that
because of (5.51) and (5.61),
T T T
’g (k) d(k) = ’g (k) ’ g (k) + βk’1 d(k’1) = g (k) g (k) ,
and concerning βk , because of (5.53), (5.64), (5.62), and the identity (5.49)
for ±k , we have
T T T T
g (k+1) g (k+1) = g (k+1) g (k) + ±k Ad(k) = ±k g (k+1) Ad(k) = βk g (k) g (k)
and hence the assumption. The algorithm is summarized in Table 5.2.
Indeed, the algorithm de¬nes conjugate directions:
Theorem 5.14 If g (k’1) = 0, then d(k’1) = 0 and the d(0) , . . . , d(k’1) are
conjugate.
5.2. Gradient and Conjugate Gradient Methods 223


Choose any x(0) ∈ Rm and calculate
d(0) := ’g (0) = b ’ Ax(0) .
For k = 0, 1, . . . put
T
g (k) g (k)
±k = ,
T
d(k) Ad(k)
x(k+1) x(k) + ±k d(k) ,
=
g (k+1) g (k) + ±k Ad(k) ,
=
T
g (k+1) g (k+1)
βk = ,
T
g (k) g (k)
d(k+1) = ’g (k+1) + βk d(k) ,
until the termination criterion (“|g (k+1) |2 = 0”) is ful¬lled.

Table 5.2. CG method.


Proof: The proof is done by induction:
The case k = 1 is clear. Assume that d(0) , . . . , d(k’1) are all nonzero and
conjugate. Thus according to Theorem 5.12 and Theorem 5.13 the identities
(5.60)“(5.64) hold up to index k. Let us ¬rst prove that d(k) = 0:
Due to g (k) + d(k) = βk’1 d(k’1) ∈ Kk (A; g (0) ) the assertion d(k) = 0
would imply directly g (k) ∈ Kk (A; g (0) ). But relations (5.63) and (5.64)
imply for the index k,
T
g (k) x = 0 for all x ∈ Kk (A; g (0) ) ,

which contradicts g (k) = 0.
T
In order to prove d(k) Ad(i) = 0 for i = 0, . . . , k ’ 1, according to (5.62)
we have to prove only the case i ¤ k ’ 2. We have
T T T
d(i) Ad(k) = ’d(i) Ag (k) + βk’1 d(i) Ad(k’1) .

The ¬rst term disappears due to Ad(i) ∈ A Kk’1 A; g (0) ‚ Kk A; g (0) ,
which means that Ad(i) ∈ span d(0) , . . . , d(k’1) , and (5.60). The second
2
term disappears because of the induction hypothesis.

Methods that aim at minimizing the error or residual on Kk A; g (0)
with respect to a norm · are called Krylov subspace methods. Here the
error will be minimized in the energy norm · = · A according to (5.55)
and Theorem 5.12.
Due to the representation of the Krylov space in Theorem 5.13 the
elements y ∈ x(0) + Kk A; g (0) are exactly the vectors of the form
y = x(0) + q(A)g (0) , for any q ∈ Pk’1 (for the notation q(A) see Appendix
224 5. Iterative Methods for Systems of Linear Equations

A.3). Hence it follows that
y ’ x = x(0) ’ x + q(A)A x(0) ’ x = p(A) x(0) ’ x ,
with p(z) = 1 + q(z)z, i.e., p ∈ Pk and p(0) = 1. On the other hand,
any such polynomial can be represented in the given form (de¬ne q by
q(z) = (p(z) ’ 1) /z). Thus Theorem 5.12 implies
x(k) ’ x ¤ y’x = p(A) x(0) ’ x (5.66)
A
A A

for any p ∈ Pk with p(0) = 1.
Let z1 , . . . , zm be an orthonormal basis of eigenvectors, that is,
T
Azj = »j zj and zi zj = δij for i, j = 1, . . . , m . (5.67)
m
Then we have x(0) ’ x = for certain cj ∈ R, and hence
j=1 cj zj
m
’x =
(0)
p(A) x p (»j ) cj zj
j=1

and therefore
m m
2 T 2
’x ’x ’x = »j |cj |
(0) (0) (0) T
x =x Ax ci cj zi Azj =
A
i,j=1 j=1

and analogously
m 2
2 2
2
’x »j |cj p(»j )| ¤ max |p(»i )| x(0) ’ x
(0)
p(A) x = .
A A
i=1,...,m
j=1
(5.68)
Relations (5.66), (5.68) imply the following theorem:
Theorem 5.15 For the CG method and any p ∈ Pk satisfying p(0) = 1,
we have
x(k) ’ x ¤ max |p(»i )| x(0) ’ x ,
A A
i=1,...,m

with the eigenvalues »1 , . . . , »m of A.
If the eigenvalues of A are not known, but their location is, i.e., if one
knows a, b ∈ R such that
a ¤ »1 , . . . , »m ¤ b , (5.69)
then only the following weaker estimate can be used:
x(k) ’ x ¤ max |p(»)| x(0) ’ x . (5.70)
A A
»∈[a,b]

Therefore, we have to ¬nd p ∈ Pm with p(0) = 1 that minimizes
max |p(»)| » ∈ [a, b] .
5.2. Gradient and Conjugate Gradient Methods 225

This approximation problem in the maximum norm appeared already in
(5.43), because there is a bijection between the sets p ∈ Pk p(1) = 1
and p ∈ Pk p(0) = 1 through
p’p , p(t) := p(1 ’ t) .
˜ ˜ (5.71)
Its solution can represented by using the Chebyshev polynomials of the
¬rst kind (see, for example, [38, p. 302]). They are recursively de¬ned by
Tk+1 (x) := 2xTk (x) ’ Tk’1 (x) for x ∈ R
T0 (x) := 1 , T1 (x) := x ,
and have the representation
Tk (x) = cos(k arccos(x))
for |x| ¤ 1. This immediately implies
|Tk (x)| ¤ 1 for |x| ¤ 1 .
A further representation, valid for x ∈ R, is
1/2 k 1/2 k
1
x + x2 ’ 1 + x ’ x2 ’ 1
Tk (x) = . (5.72)
2
The optimal polynomial in (5.70) is then de¬ned by
Tk ((b + a ’ 2z)/(b ’ a))
for z ∈ R .
p(z) :=
Tk ((b + a)/(b ’ a))
This implies the following result:
Theorem 5.16 Let κ be the spectral condition number of A and assume
κ > 1. Then
k
κ1/2 ’ 1
1
’x ¤ ’x ¤2 x(0) ’ x
(k) (0)
x x . (5.73)
A A A
κ1/2 + 1
κ+1
Tk κ’1


Proof: Choose a as the smallest eigenvalue »min and b as the largest one
»max .
The ¬rst inequality follows immediately from (5.70) and κ = b/a. For
the second inequality note that due to (κ + 1)/(κ ’ 1) = 1 + 2/(κ ’ 1) =:
1 + 2· ≥ 1, (5.72) implies
1/2 k
κ+1 1
≥ 1 + 2· + (1 + 2·)2 ’ 1
Tk
κ’1 2
k
1 1/2
= 1 + 2· + 2 (·(· + 1)) .
2
Finally,
(· + 1)1/2 + · 1/2
2
1/2 1/2 1/2
1 + 2· + 2 (·(· + 1)) = · + (· + 1) =
(· + 1)1/2 ’ · 1/2
226 5. Iterative Methods for Systems of Linear Equations

(1 + 1/·)1/2 + 1
= ,
1/2
’1
(1 + 1/·)
2
which concludes the proof because of 1 + 1/· = κ.

For large κ we have again
κ1/2 ’ 1 2
≈ 1 ’ 1/2 .
κ1/2 + 1 κ
Compared with (5.58), κ has been improved to κ1/2 .
From (5.4) and (5.34) the complexity of the ¬ve-point stencil discretiza-
tion of the Poisson equation on the square results in
1
O κ1/2 O(m) = O(n) O(m) = O m3/2 .
ln
µ
This is the same behaviour as that of the SOR method with optimal re-
laxation parameter. The advantage of the above method lies in the fact
that the determination of parameters is not necessary for applying the
CG method. For quasi-uniform triangulations, Theorem 3.45 implies an
analogous general statement.
A relation to the semi-iterative methods follows from (5.71): The estimate
(5.66) can also be expressed as
¤ p(I ’ A)e(0)
e(k) (5.74)
A A

for any p ∈ Pk with p(1) = 1.
This is the same estimate as (5.42) for the Richardson iteration (5.28) as
basis method, with the Euclidean norm |·|2 replaced by the energy norm ·
A . While the semi-iterative methods are de¬ned by minimization of upper
bounds in (5.42), the CG method is optimal in the sense of (5.74), without
knowledge of the spectrum σ(I ’ A). In this manner the CG method can
be seen as an (optimal) acceleration method for the Richardson iteration.



Exercises
5.4 Let A ∈ Rm,m be a symmetric positive de¬nite matrix.
(a) Show that for x, y with xT y = 0 we have
κ’1
x, y
¤
A
,
x y κ+1
A A

where κ denotes the spectral condition number of A.
Hint: Represent x, y in terms of an orthonormal basis consisting of
eigenvectors of A.
5.3. Preconditioned Conjugate Gradient Method 227

(b) Show using the example m = 2 that this estimate is sharp. To this
end, look for a positive de¬nite symmetric matrix A ∈ R2,2 as well
as vectors x, y ∈ R2 with xT y = 0 and
κ’1
x, y A
= .
x y κ+1
A A


5.5 Prove that the computation of the conjugate direction in the CG
method in the general step k ≥ 2 is equivalent to the three-term recursion
formula
d(k+1) = [±k A + (βk + 1)I] d(k) ’ βk’1 d(k’1) .

5.6 Let A ∈ Rm,m be a symmetric positive de¬nite matrix with spectral
condition number κ. Suppose that the spectrum σ(A) of the matrix A
satis¬es a0 ∈ σ(A) as well as σ(A) \ {a0 } ‚ [a, b] with 0 < a0 < a ¤ b.
Show that this yields the following convergence estimate for the CG
method:
√ k’1
b ’ a0 κ’1
ˆ

x(k) ’ x A ¤ 2 x(0) ’ x A ,
a0 κ+1
ˆ
where κ := b/a ( < κ ).
ˆ



5.3 Preconditioned Conjugate Gradient Method
Due to Theorem 5.16, κ(A) should be small or only weakly growing in m,
which is not true for a ¬nite element sti¬ness matrix.
The technique of preconditioning is used ” as already discussed in Sec-
tion 5.1 ” to transform the system of equations in such a way that the
condition number of the system matrix is reduced without increasing the
e¬ort in the evaluation of the matrix vector product too much.
In a preconditioning from the left the system of equations is transformed
to
C ’1 Ax = C ’1 b
with a preconditioner C; in a preconditioning from the right it is transformed
to
AC ’1 y = b ,
such that x = C ’1 y is the solution of (5.1). Since the matrices are generally
sparse, this always has to be interpreted as a solution of the system of
equations Cx = y.
If A is symmetric and positive de¬nite, then this property is generally
violated by the transformed matrix for both variants, even for a symmetric
228 5. Iterative Methods for Systems of Linear Equations

positive de¬nite C. We assume for a moment to have a decomposition of
C with a nonsingular matrix W as
C = WWT .
Then Ax = b can be transformed to W ’1 AW ’T W T x = W ’1 b, i.e., to
B = W ’1 AW ’T , c = W ’1 b .
By = c with (5.75)
The matrix B is symmetric and positive de¬nite. The solution x is then
given by x = W ’T y. This procedure is called split preconditioning.
Due to W ’T BW T = C ’1 A and W BW ’1 = AC ’1 , B, C ’1 A and AC ’1
have the same eigenvalues, and therefore also the same spectral condition
number κ. Therefore, C should be “close” to A in order to reduce the
condition number. The CG method, applied to (5.75) and then back trans-
formed, leads to the preconditioned conjugate gradient method (PCG):
The terms of the CG method applied to (5.75) will all be marked by ˜,
with the exception of ±k and βk .
Due to the back transformation
x = W ’T x
˜
the algorithm has the search direction
d(k) := W ’T d(k)
˜

for the transformed iterate
x(k) := W ’T x(k) .
˜ (5.76)
The gradient g (k) of (5.44) in x(k) is given by
g (k) := Ax(k) ’ b = W B x(k) ’ c = W g (k) ,
˜ ˜
and hence
˜
g (k+1) = g (k) + ±k W B d(k) = g (k) + ±k Ad(k) ,
so that this formula remains unchanged compared with the CG method
with a new interpretation of the search direction. The search directions are
updated by
d(k+1) = ’W ’T W ’1 g (k+1) + βk d(k) = ’C ’1 g (k+1) + βk d(k) ,
so that in each iteration step additionally the system of equations
Ch(k+1) = g (k+1) has to be solved.
Finally, we have
T T
g (k) g (k) = g (k) C ’1 g (k) = g (k) h(k)
T
˜ ˜
and
T
˜T ˜
d(k) B d(k) = d(k) Ad(k) ,
so that the algorithm takes the form of Table 5.3.
5.3. Preconditioned Conjugate Gradient Method 229


Choose any x(0) ∈ Rm and calculate
d(0) := ’h(0) := ’C ’1 g (0) .
g (0) = Ax(0) ’ b ,
For k = 0, 1, . . . put
T
g (k) h(k)
±k = ,
T
d(k) Ad(k)
x(k+1) x(k) + ±k d(k) ,
=
g (k+1) g (k) + ±k Ad(k) ,
=
C ’1 g (k+1) ,
h(k+1) =
T
g (k+1) h(k+1)
βk = ,
T
g (k) h(k)
’h(k+1) + βk d(k) ,
d(k+1) =
up to the termination criterion (“|g (k+1) |2 = 0”) .

Table 5.3. PCG method.


The solution of the additional systems of equations for sparse matrices
should have the complexity O(m), in order not to worsen the complexity
for an iteration. It is not necessary to know a decomposition C = W W T .
Alternatively, the PCG method can be established by noting that C ’1 A
is self-adjoint and positive de¬nite with respect to the energy scalar product
·, · C de¬ned by C:
T
C ’1 Ax, y = C ’1 Ax Cy = xT Ay = xT C(C ’1 Ay) = x, C ’1 Ay
C C

and hence also C ’1 Ax, x C > 0 for x = 0.
Choosing the CG method for (5.75) with respect to ·, · C , we obtain
precisely the above method.
In case the termination criterion “ g (k+1) 2 = 0” is used for the iteration,
the scalar product must be additionally calculated. Alternatively, we may
T
use “ g (k+1) h(k+1) = 0”. Then the residual is measured in the norm
· C ’1 .
Following the reasoning at the end of Section 5.2, the PCG method can be
interpreted as an acceleration of a linear stationary method with iteration
matrix

M = I ’ C ’1 A .

For a consistent method, we have N = C ’1 or, in the formulation (5.10),
W = C. This observation can be extended in such a way that the CG
method can be used for the acceleration of iteration methods, for example
also for the multigrid method, which will be introduced in Section 5.5. Due
230 5. Iterative Methods for Systems of Linear Equations

to the deduction of the preconditioned CG method and the identity
x(k) ’ x = x(k) ’ x
˜ ˜ ,
A B

which results from the transformation (5.76), the approximation properties
for the CG method also hold for the PCG method if the spectral condition
number κ(A) is replaced by κ(B) = κ(C ’1 A). Therefore,
k
κ1/2 ’ 1
’x ¤2 x(0) ’ x
(k)
x A A
κ1/2 + 1
with κ = κ(C ’1 A).
There is a close relation between those preconditioning matrices C, which
keep κ(C ’1 A) small, and well-convergent linear stationary iteration meth-
ods with N = C ’1 (and M = I ’ C ’1 A) if N is symmetric and positive
de¬nite. Indeed,
κ(C ’1 A) ¤ (1 + (M ))/(1 ’ (M ))
if the method de¬ned by M and N is convergent and N is symmetric for
symmetric A (see Exercise 5.7).
From the considered linear stationary methods because of the required
symmetry we may take
• Jacobi™s method:
This corresponds exactly to the diagonal scaling, which means the division
of each equation by its diagonal element. Indeed, from the decomposition
(5.18) we have C = N ’1 = D, and the PCG method is equivalent to the
preconditioning from the left by the matrix C ’1 in combination with the
usage of the energy scalar product ·, · C .
• The SSOR method:
According to (5.38) we have
C = ω ’1 (2 ’ ω)’1 (D + ωL)D’1 (D + ωLT ) .
Hence C is symmetric and positive de¬nite. The solution of the auxiliary
systems of equations needs only forward and backward substitutions with
the same structure of the matrix as for the system matrix, so that the
requirement of lower complexity is also ful¬lled. An exact estimation of
κ(C ’1 A) shows (see [3, pp. 328 ¬.]) that under certain requirements for A,
which re¬‚ect properties of the boundary value problem and the discretiza-
tion, we ¬nd a considerable improvement of the conditioning by using an
estimate of the type
κ(C ’1 A) ¤ const(κ(A)1/2 + 1) .
The choice of the relaxation parameter ω is not critical. Instead of try-
ing to choose an optimal one for the contraction number of the SSOR
5.3. Preconditioned Conjugate Gradient Method 231

method, we can minimize an estimation for κ(C ’1 A) (see [3, p. 337]),
which recommends a choice of ω in [1.2, 1.6].
For the ¬ve-point stencil discretization of the Poisson equation on the
square we have, according to (5.34), κ(A) = O(n2 ), and the above con-
ditions are ful¬lled (see [3, pp. 330 f.]). By SSOR preconditioning this is
improved to κ(C ’1 A) = O(n), and therefore the complexity of the method
is
1 1
O κ1/2 O(m) = ln O n1/2 O(m) = O m5/4 .
ln (5.77)
µ µ
As discussed in Section 2.5, direct elimination methods are not suitable in
conjunction with the discretization of boundary value problems with large
node numbers, because in general ¬ll-in occurs. As discussed in Section 2.5,
L = (lij ) describes a lower triangular matrix with lii = 1 for all i = 1, . . . , m
(the dimension is described there with the number of degrees of freedom
M ) and U = (uij ) an upper triangular matrix. The idea of the incomplete
LU factorization, or ILU factorization, is to allow only certain patterns
E ∈ {1, . . . , m}2 for the entries of L and U , and instead of A = LU , in
general we can require only
A = LU ’ R.
Here the remainder R = (rij ) ∈ Rm,m has to satisfy
for (i, j) ∈ E .
rij = 0 (5.78)
The requirements
m
for (i, j) ∈ E
aij = lik ukj (5.79)
k=1

mean |E| equations for the |E| entries of the matrices L and U . (Notice that
lii = 1 for all i.) The existence of such factorizations will be discussed later.
Analogously to the close connection between the LU factorization and
an LDLT or LLT factorization for symmetric or symmetric positive def-
inite matrices, as de¬ned in Section 2.5, we can use the IC factorization
(incomplete Cholesky factorization) for such matrices. The IC factorization
needs a representation in the following form:
A = LLT ’ R .
Based on an ILU factorization a linear stationary method is de¬ned by
N = (LU )’1 (and M = I ’ N A), the ILU iteration. We thus have an
expansion of the old method of iterative re¬nement.
Using C = N ’1 = LU for the preconditioning, the complexity of the
auxiliary systems depends on the choice of the matrix pattern E. In general,
the following is required:
E := (i, j) aij = 0 , i, j = 1, . . . , m ‚ E , (i, i) i = 1, . . . , m ‚ E .
(5.80)
232 5. Iterative Methods for Systems of Linear Equations

The requirement of equality E = E is most often used. Then, and also in the
case of ¬xed expansions of E , it is ensured that for a sequence of systems
of equations with a matrix A that is sparse in the strict sense, this will also
hold for L and U . All in all, only O(m) operations are necessary, including
the calculation of L and U , as in the case of the SSOR preconditioning
for the auxiliary system of equations. On the other hand, the remainder R
should be rather small in order to ensure a good convergence of the ILU
iteration and also to ensure a small spectral condition number κ(C ’1 A).
Possible matrix patterns E are shown, for example, in [28, pp. 275 ¬.], where
a more speci¬c structure of L and U is discussed if the matrix A is created
by a discretization on a structured grid, for example by a ¬nite di¬erence
method.
The question of the existence (and stability) of an ILU factorization
remains to be discussed. It is known from (2.56) that also for the existence
of an LU factorization certain conditions are necessary, as for example the
M-matrix property. This is even su¬cient for an ILU factorization.

Theorem 5.17 Let A ∈ Rm,m be an M-matrix. Then for a given pat-
tern E that satis¬es (5.80), an ILU factorization exists. The hereby de¬ned
decomposition of A as A = LU ’ R is regular in the following sense:

(LU )’1 ≥ 0, (R)ij ≥ 0 for all i, j = 1, . . . , m .
ij



2
Proof: See [16, p. 235].


An ILU (or IC) factorization can be de¬ned by solving the equations
(5.78) for lij and uij in an appropriate order. Alternatively, the elimination
or Cholesky method can be used in its original form on the pattern E.
An improvement of the eigenvalue distribution of C ’1 A is sometimes
possible by using an MIC factorization (modi¬ed incomplete Cholesky fac-
torization) instead of an IC factorization. In contrast to (5.79) the updates
in the elimination method for positions outside the pattern are not ignored
here but have to be performed for the corresponding diagonal element.
Concerning the reduction of the condition number by the ILU or IC
preconditioning for the model problem, we have the same situation as for
the SSOR preconditioning. In particular (5.77) holds, too.
The auxiliary system of equations with C = N ’1 , which means that

h(k+1) = N g (k+1) ,

can also be interpreted as an iteration step of the iteration method de¬ned
by N with initial value z (0) = 0 and right-hand side g (k+1) . An expansion
of the discussed possibilities for preconditioning is therefore obtained by
using a ¬xed number of iteration steps instead of only one.
5.4. Krylov Subspace Methods for Nonsymmetric Systems of Equations 233

Exercises
5.7 Let A1 , A2 , . . . , Ak , C1 , C2 , . . . , Ck ∈ Rm,m be symmetric positive
semide¬nite matrices with the property
axT Ci x ¤ xT Ai x ¤ bxT Ci x for x ∈ Rm , i = 1, . . . , k and 0 < a ¤ b .
Prove: If A := k Ai and C := k Ci are positive de¬nite, then the
i=1 i=1
spectral condition number κ of C ’1 A satis¬es
b
κ(C ’1 A) ¤ .
a

5.8 Show that the matrix
« 
2 1 1
1 1
2
A :=
1 1 2
is positive de¬nite and its spectral condition number is 4.
Hint: Consider the associated quadratic form.

5.9 Investigate the convergence of the (P)CG method on the basis of
Theorem 3.45 and distinguish between d = 2 and d = 3.



5.4 Krylov Subspace Methods
for Nonsymmetric Systems of Equations
With the di¬erent variants of the PCG method we have methods that
are quite appropriate ” regarding their complexity ” for those systems
of equations that arise from the discretization of boundary value prob-
lems. However, this holds only under the assumption that the system
matrix is symmetric and positive de¬nite, reducing the possibilities of ap-
plication, for example to ¬nite element discretizations of purely di¬usive
processes without convective transport mechanism (see (3.23)). Exceptions
for time-dependent problems are only the (semi-)explicit time discretization
(compare (7.72)) and the Lagrange“Galerkin method (see Section 9.4). For
all other cases the systems of equations that arise are always nonsymmetric
and positive real, which means that the system matrix A satis¬es
A + AT is positive de¬nite.
It is desirable to generalize the (P)CG methods for such matrices. The CG
method is characterized by two properties:
• The iterate x(k) minimizes f (·) = · ’x on x(0) + Kk A; g (0) ,
A
where x = A’1 b.
234 5. Iterative Methods for Systems of Linear Equations

• The basis vectors d(i) , i = 0, . . . , k ’ 1, of Kk A; g (0) do not have to
be calculated in advance (and stored in the computer), but will be
calculated by a three-term recursion (5.61) during the iteration. An
analogous relation holds by de¬nition for x(k) (see (5.48)).
The ¬rst property can be preserved in the following, whereby the norm
of the error or residual minimization varies in each method. The second
property is partially lost, because generally all basis vectors d(0) , . . . , d(k’1)
are necessary for the calculation of x(k) . This will result in memory space
problems for large k. As for the CG methods, preconditioning will be nec-
essary for an acceptable convergence of the methods. The conditions for
the preconditioning matrices are the same as for the CG method with the
exception of symmetry and positive de¬niteness. All three methods of pre-
conditioning are in principle possible. Therefore, preconditioning will not
be discussed in the following; we refer to Section 5.3.
The simplest approach is the application of the CG method to a system
of equations with symmetric positive de¬nite matrix equivalent to (5.1).
This is the case for the normal equations
AT Ax = AT b . (5.81)
The approach is called CGNR (Conjugate Gradient Normal Residual),
because here the iterate x(k) minimizes the Euclidean norm of the residual
on x(0) + Kk AT A; g (0) with g (0) = AT Ax(0) ’ b . This follows from the
equation
y’x = (Ay ’ b)T (Ay ’ b) = |Ay ’ b|2
2
(5.82)
AT A 2

for any y ∈ Rm and the solution x = A’1 b.
All advantages of the CG method are preserved, although in (5.53) and
(5.65) Ad(k) is to be replaced by AT Ad(k) . Additionally to the doubling of
the number of operations this may be a disadvantage if κ2 (A) is large, since
κ2 (AT A) = κ2 (A)2 can lead to problems of stability and convergence. Due
to (5.34) this is to be expected for a large number of degrees of freedom.
Furthermore, in the case of list-based storage one of the operations Ay
and AT y is always very expensive due to searching. It is even possible
that we do not explicitly know the matrix A but that only the mapping
y ’ Ay can be evaluated, which then disquali¬es this method completely
(see Exercise 8.6).
The same drawback occurs if
AAT x = b
˜ (5.83)
with the solution x = A’T x taken instead of (5.81). If x(k) is the kth iterate
˜ ˜
of the CG method applied to (5.83), then the x := AT x(k) minimizes the
(k)
˜
T T (0)
residual in the Euclidean norm on x0 + A Kk AA ; g : Note that
T 2
2
y’x = AT y ’ x AT y ’ x = AT y ’ x
˜˜ ˜ ˜ ˜
AAT 2
5.4. Krylov Subspace Methods for Nonsymmetric Equations 235


Let g (0) ∈ Rm , g (0) = 0 be given, Set
v1 := g (0) / |g (0) |2 .
For j = 1, . . . , k calculate
T
hij := vi Avj for i = 1, . . . , j ,
j
Avj ’
wj := hij vi ,
i=1
|wj |2 .
hj+1,j :=
If hj+1,j = 0, termination; otherwise, set
vj+1 := wj /hj+1,j .

Table 5.4. Arnoldi algorithm.

holds for any y ∈ Rm and g (0) = Ax(0) ’ b. This explains the terminology
˜
CGNE (with E for Error).
Whether a method minimizes the error of the residual obviously depends
on the norm used. For a symmetric positive de¬nite B ∈ Rm,m , any y ∈ Rm ,
and x = A’1 b, we have
Ay ’ b = y’x .
AT BA
B

For B = A’T and a symmetric positive de¬nite A we get the situation of
the CG method:
Ay ’ b = y’x .
A’T A

For B = I we get again (5.82):
|Ay ’ b|2 = y ’ x AT A .

The minimization of this functional on x(0) + Kk A; g (0) (not
Kk AT A; g (0) ) leads to the GMRES method (Generalized Minimum
RESidual).
This (and other) methods are founded algorithmically on the recur-
sive construction of orthonormal bases of Kk A; g (0) by Arnoldi™s method.
This method combines the generation of a basis according to (5.61) and
Schmidt™s orthonormalization (see Table 5.4).
If Arnoldi™s method can be performed up to the index k, then
hij := 0 for j = 1, . . . , k, i = j + 2, . . . , k + 1 ,
(hij )ij ∈ Rk,k ,
Hk :=
(hij )ij ∈ Rk+1,k ,
¯
Hk :=
(v1 , . . . , vk+1 ) ∈ Rm,k+1 .
Vk+1 :=
The matrix Hk is an upper Hessenberg matrix (see Appendix A.3). The
basis for the GMRES method is the following theorem:
236 5. Iterative Methods for Systems of Linear Equations

Theorem 5.18 If Arnoldi™s method can be performed up to the index k,
then
(1) v1 , . . . , vk+1 form an orthonormal basis of Kk+1 (A; g (0) ).
(2)
¯
AVk = Vk Hk + wk eT = Vk+1 Hk , (5.84)
k

with ek = (0, . . . , 0, 1)T ∈ Rk ,
VkT AVk = Hk . (5.85)
(3) The problem
|Ay ’ b|2 y ∈ x(0) + Kk (A; g (0) )
Minimize for
with minimum x(k) is equivalent to
Hk ξ ’ βe1 ξ ∈ Rk
¯
Minimize for (5.86)
2

with β := ’ g (0) and minimum ξ (k) , and we have
2

x(k) = x(0) + Vk ξ (k) .
If Arnoldi™s method terminates at the index k, then
x(k) = x = A’1 b .

Proof: (1): The vectors v1 , . . . , vk+1 are orthonormal by construction;
hence we have only to prove vi ∈ Kk+1 A; g (0) for i = 1, . . . , k + 1. This
follows from the representation
with polynomials qi’1 ∈ Pi’1 .
vi = qi’1 (A)v1
In this form we can prove the statement by induction with respect to k.
For k = 0 the assertion is trivial. Let the statement hold for k ’ 1. The
validity for k then follows from
k k
hk+1,k vk+1 = Avk ’ Aqk’1 (A) ’
hik vi = hik qi’1 (A) v1 .
i=1 i=1

(2): Relation (5.85) follows from (5.84) by multiplication by VkT , since
VkT Vk = I and VkT wk = hk+1,k VkT vk+1 = 0 due to the orthonormality
of the vi .
The relation in (5.84) is the matrix representation of
j j+1
Avj = hij vi + wj = hij vi for j = 1, . . . , k .
i=1 i=1

(3): Due to (1), the space x(0) + Kk A; g (0) has the parametrisation

ξ ∈ Rk .
y = x(0) + Vk ξ with (5.87)
5.4. Krylov Subspace Methods for Nonsymmetric Equations 237

The assertion is a consequence of the identity
Ay ’ b = A x(0) + Vk ξ ’ b = AVk ξ + g (0)
Vk+1 Hk ξ ’ βv1 = Vk+1 Hk ξ ’ βe1 ,
¯ ¯
=
which follows from (2), since it implies
|Ay ’ b|2 = Vk+1 (Hk ξ ’ βe1 ) = Hk ξ ’ βe1
¯ ¯
2 2

due to the orthogonality of Vk+1 . The last assertion ¬nally can be seen in
this way: If Arnoldi™s method breaks down at the index k, then relation (2)
becomes
AVk = Vk Hk ,
and
¯
AVk = Vk+1 Hk
will further hold with vk+1 chosen arbitrarily (due to hk+1,k = 0). Since A
is nonsingular, this also holds for Hk . Hence the choice
’1
ξ := Hk (βe1 ) ,
which satis¬es
Hk ξ ’ βe1 = |Hk ξ ’ βe1 |2 = 0 ,
¯
2

is possible. Hence the corresponding y ∈ Rm de¬ned by (5.87) ful¬lls y =
2
x(k) = x.

One problem of Arnoldi™s method is that the orthogonality of the vi is
easily lost due to rounding errors. If one substitutes the assignment
j
wj := Avj ’ hij vi
i=1

in Table 5.4 by the operations
wj := Avj ,
for i = 1, . . . , j calculate
T
hij := wj vi ,
wj := wj ’ hij vi ,
which de¬ne the same vector, one obtains the modi¬ed Arnoldi™s method.
From this relation and from (5.86) the GMRES method is constructed in
its basic form. Alternatively, Schmidt™s orthonormalization can be replaced
by the Householder method (see [28, pp. 159 ¬.]). With exact arithmetic
the GMRES algorithm terminates only after reaching the exact solution
(with hk+1,k = 0). This is not always the case for alternative methods
of the same class. For an increasing iteration index k and large problem
dimensions m there may be lack of enough memory for the storage of
238 5. Iterative Methods for Systems of Linear Equations

the basis vectors v1 , . . . , vk . A remedy is o¬ered by working with a ¬xed
number n of iterations and then to restart the algorithm with x(0) := x(n)
and g (0) := Ax(0) ’ b, until ¬nally the convergence criterion is ful¬lled
(GMRES method with restart). There is also a truncated version of the
GMRES method, in which only the last n basis vectors are used. The
minimization of the error in the energy norm (on the vector space K)
as with the CG method makes sense only for symmetric positive de¬nite
matrices A. But the variational equation
(Ay ’ b)T z = 0 for all z ∈ K
that characterizes this minimum in general can be taken as de¬ning con-
dition for y. Further variants of Krylov subspace methods rely on this.
Another large class of such methods is founded on the Lanczos biorthogo-
nalization, in which apart from a basis v1 , . . . , vk of Kk (A; v1 ) another basis
w1 , . . . , wk of Kk (AT ; w1 ) is constructed, such that
T
vj wi = δij for i, j = 1, . . . , k .
The best-known representative of this method is the BICGSTAB method.
For further discussion of this topic see, for example, [28].


Exercises
5.10 Consider the linear system Ax = b, where A = ±Q for some ± ∈
R \ {0} and some orthogonal matrix Q. Show that, for an arbitrary initial
iterate x(0) , the CGNE method terminates after one step with the exact
solution.

5.11 Provided that Arnoldi™s method can be performed up to the index
k, show that it is possible to incorporate a convergence test of the GMRES
method without computing the approximate solution explicitely, i.e., prove
the following formulas:
g (k) := Ax(k) ’ b = hk+1,k eT ξ (k) vk+1 ,
k
|g (k) |2 = hk+1,k |eT ξ (k) | .
k



5.5 The Multigrid Method
5.5.1 The Idea of the Multigrid Method
We discuss again the model problem of the ¬ve-point stencil discretization
for the Poisson equation on the square and use the relaxed Jacobi™s method.
Then due to (5.31) the iteration matrix is
ω
M = ωMJ + (1 ’ ω)I = I ’ A ,
4
5.5. Multigrid Method 239

with A being the sti¬ness matrix according to (1.14). For ω = ω/4 this
˜
coincides with the relaxed Richardson method, which according to (5.35)
has the poor convergence behaviour of Jacobi™s method, even for optimal
choice of the parameter. Nevertheless, for a suitable ω the method has
positive properties. Due to (5.25) the eigenvalues of M are
ω kπ lπ
»k,l = 1 ’ ω + 1 ¤ k, l ¤ n ’ 1 .
cos + cos ,
2 n n
This shows that there is a relation between the size of the eigenvalues and
the position of the frequency of the assigned eigenfunction depending on
the choice of ω: For ω = 1, which is Jacobi™s method, (M ) = »1,1 =
’»n’1,n’1 . Thus the eigenvalues are large if k and l are close to 1 or n.
Hence there are large eigenvalues for eigenfunctions with low frequency as
well as for eigenfunctions with high frequency. For ω = 1 , however, we have
2
(M ) = »1,1 , and the eigenvalues are large only in the case that k and l
are near to 1, which means that the eigenfunctions have low frequency.
In general, if the error of the iterate e(k) had a representation in terms of
orthonormal eigenvectors zν with small eigenvalues, as for example |»ν | ¤
1
2,

e(k) = cν z ν ,
ν:|»ν |¤ 1
2


then according to (5.11) it would follow for the error measured in the
Euclidean vector norm | · |2 that
« 1/2
 »2 c2 

<<

. 9
( 16)



>>