± ’ 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