ation parameter, is also called the gradient method with dynamic parameter

(shortly, gradient method), to distinguish it from the stationary Richardson

method (4.23) or gradient method with constant parameter, where ±k = ±

is a constant for any k ≥ 0.

Summarizing, the gradient method can be described as follows:

148 4. Iterative Methods for Solving Linear Systems

given x(0) ∈ Rn , for k = 0, 1, . . . until convergence, compute

r(k) = b ’ Ax(k)

T

r(k) r(k)

±k = T

r(k) Ar(k)

x(k+1) = x(k) + ±k r(k) .

Theorem 4.10 Let A be a symmetric and positive de¬nite matrix; then

the gradient method is convergent for any choice of the initial datum x(0)

and

K2 (A) ’ 1 (k)

¤

e(k+1) e A, k = 0, 1, . . . , (4.38)

A

K2 (A) + 1

·

where is the energy norm de¬ned in (1.28).

A

Proof. Let x(k) be the solution generated by the gradient method at the k-th

(k+1)

step. Then, let xR be the vector generated by taking one step of the non

preconditioned Richardson method with optimal parameter starting from x(k) ,

(k+1)

= x(k) + ±opt r(k) .

i.e., xR

Due to Corollary 4.1 and (4.28), we have

K2 (A) ’ 1 (k)

(k+1)

¤

eR e A,

A

K2 (A) + 1

(k+1) (k+1)

’ x. Moreover, from (4.35) we have that the vector x(k+1) ,

where eR = xR

generated by the gradient method, is the one that minimizes the A-norm of

the error among all vectors of the form x(k) + θr(k) , with θ ∈ R. Therefore,

(k+1)

e(k+1) A ¤ eR 3

A which is the desired result.

We notice that the line through x(k) and x(k+1) is tangent at the point

x(k+1) to the ellipsoidal level surface x ∈ Rn : ¦(x) = ¦(x(k+1) ) (see

also Figure 4.5).

Relation (4.38) shows that convergence of the gradient method can be

quite slow if K2 (A) = »1 /»n is large. A simple geometric interpretation of

this result can be given in the case n = 2. Suppose that A=diag(»1 , »2 ),

with 0 < »2 ¤ »1 and b = (b1 , b2 )T .

In such a case, the curves corresponding to ¦(x1 , x2 ) = c, as c varies

in R+ , form a sequence of concentric ellipses whose semi-axes have length

inversely proportional to the values »1 and »2 . If »1 = »2 , the ellipses

degenerate into circles and the direction of the gradient crosses the center

directly, in such a way that the gradient method converges in one iteration.

Conversely, if »1 »2 , the ellipses become strongly eccentric and the

method converges quite slowly, as shown in Figure 4.5, moving along a

“zig-zag” trajectory.

4.3 Stationary and Nonstationary Iterative Methods 149

1

2

0.5

1

(1)

x (3)

x

0

0

(2)

’1

x

’0.5

(0)

x

’2

’2 0 2 ’1 ’0.5 0 0.5 1

FIGURE 4.5. The ¬rst iterates of the gradient method on the level curves of ¦

Program 19 provides an implementation of the gradient method with dy-

namic parameter. Here and in the programs reported in the remainder of

the section, the input parameters A, x, b, M, maxit and tol respectively

represent the coe¬cient matrix of the linear system, the initial datum x(0) ,

the right side, a possible preconditioner, the maximum number of admis-

sible iterations and a tolerance for the stopping test. This stopping test

checks if the ratio r(k) 2 / b 2 is less than tol. The output parameters of

the code are the the number of iterations niter required to ful¬ll the stop-

ping test, the vector x with the solution computed after niter iterations

and the normalized residual error = r(niter) 2 / b 2 . A null value of the

parameter flag warns the user that the algorithm has actually satis¬ed

the stopping test and it has not terminated due to reaching the maximum

admissible number of iterations.

Program 19 - gradient : Gradient method with dynamic parameter

function [x, error, niter, ¬‚ag] = gradient(A, x, b, M, maxit, tol)

¬‚ag = 0; niter = 0; bnrm2 = norm( b );

if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

r = b - A*x; error = norm( r ) / bnrm2;

if ( error < tol ) return, end

for niter = 1:maxit

z = M \ r; rho = (r™*z);

q = A*z; alpha = rho / (z™*q );

x = x + alpha * z; r = r - alpha*q;

error = norm( r ) / bnrm2;

if ( error <= tol ), break, end

end

if ( error > tol ) ¬‚ag = 1; end

150 4. Iterative Methods for Solving Linear Systems

Example 4.6 Let us solve with the gradient method the linear system with ma-

trix Am ∈ Rm—m generated with the MATLAB commands G=numgrid(™S™,n);

A=delsq(G) where m = (n ’ 2)2 . This matrix is associated with the discretiza-

tion of the di¬erential Laplace operator on the domain [’1, 1]2 . The right-hand

side bm is selected in such a way that the exact solution is the vector 1T ∈ Rm .

The matrix Am is symmetric and positive de¬nite for any m and becomes ill-

conditioned for large values of m. We run Program 19 in the cases m = 16 and

m = 400, with x(0) = 0T , tol=10’10 and maxit=200. If m = 400, the method

fails to satisfy the stopping test within the admissible maximum number of it-

erations and exhibits an extremely slow reduction of the residual (see Figure

4.6). Actually, K2 (A400 ) 258. If, however, we precondition the system with the

matrix P = RT Rin , where Rin is the lower triangular matrix in the Cholesky

in

incomplete factorization of A, the algorithm ful¬lls the convergence within the

maximum admissible number of iterations (indeed, now K2 (P’1 A400 ) 38). •

0

10

(c)

’2

10

’4

10

’6

10

’8

10

(a)

’10

10 (d)

(b)

’12

10

’14

10

0 50 100 150 200 250

FIGURE 4.6. The residual normalized to the starting one, as a function of the

number of iterations, for the gradient method applied to the systems in Example

4.6. The curves labelled (a) and (b) refer to the case m = 16 with the non precon-

ditioned and preconditioned method, respectively, while the curves labelled (c)

and (d) refer to the case m = 400 with the non preconditioned and preconditioned

method, respectively

4.3.4 The Conjugate Gradient Method

The gradient method consists essentially of two phases: choosing a descent

direction (the one of the residual) and picking up a point of local minimum

for ¦ along that direction. The second phase is independent of the ¬rst one

since, for a given direction p(k) , we can determine ±k as being the value

of the parameter ± such that ¦(x(k) + ±p(k) ) is minimized. Di¬erentiating

with respect to ± and setting to zero the derivative at the minimizer, yields

T

p(k) r(k)

±k = , (4.39)

T

p(k) Ap(k)

4.3 Stationary and Nonstationary Iterative Methods 151

instead of (4.37). The question is how to determine p(k) . A di¬erent ap-

proach than the one which led to identify p(k) with r(k) is suggested by the

following de¬nition.

De¬nition 4.4 A direction x(k) is said to be optimal with respect to a

direction p = 0 if

¦(x(k) ) ¤ ¦(x(k) + »p), ∀» ∈ R. (4.40)

If x(k) is optimal with respect to any direction in a vector space V, we say

that x(k) is optimal with respect to V.

From the de¬nition of optimality, it turns out that p must be orthogonal

to the residual r(k) . Indeed, from (4.40) we conclude that ¦ admits a local

minimum along p for » = 0, and thus the partial derivative of ¦ with

respect to » must vanish at » = 0. Since

‚¦ (k)

(x + »p) = pT (Ax(k) ’ b) + »pT Ap,

‚»

we therefore have

‚¦ (k)

pT (r(k) ) = 0,

(x )|»=0 = 0 i¬

‚»

that is, p ⊥ r(k) . Notice that the iterate x(k+1) of the gradient method

is optimal with respect to r(k) since, due to the choice of ±k , we have

r(k+1) ⊥ r(k) , but this property no longer holds for the successive iterate

x(k+2) (see Exercise 12). It is then natural to ask whether there exist descent

directions that maintain the optimality of iterates. Let

x(k+1) = x(k) + q,

and assume that x(k) is optimal with respect to a direction p (thus, r(k) ⊥

p). Let us impose that x(k+1) is still optimal with respect to p, that is,

r(k+1) ⊥ p. We obtain

0 = pT r(k+1) = pT (r(k) ’ Aq) = ’pT Aq.

The conclusion is that, in order to preserve optimality between succes-

sive iterates, the descent directions must be mutually A-orthogonal or A-

conjugate, i.e.

pT Aq = 0.

A method employing A-conjugate descent directions is called conjugate.

The next step is how to generate automatically a sequence of conjugate

152 4. Iterative Methods for Solving Linear Systems

directions. This can be done as follows. Let p(0) = r(0) and search for the

directions of the form

p(k+1) = r(k+1) ’ βk p(k) , k = 0, 1, . . . (4.41)

where βk ∈ R must be determined in such a way that

(Ap(j) )T p(k+1) = 0, j = 0, 1, . . . , k. (4.42)

Requiring that (4.42) is satis¬ed for j = k, we get from (4.41)

(Ap(k) )T r(k+1)

βk = , k = 0, 1, . . .

(Ap(k) )T p(k)

We must now verify that (4.42) holds also for j = 0, 1, . . . , k ’1. To do this,

let us proceed by induction on k. Due to the choice of β0 , relation (4.42)

holds for k = 0; let us thus assume that the directions p(0) , . . . , p(k’1) are

mutually A-orthogonal and, without losing generality, that

(p(j) )T r(k) = 0, j = 0, 1, . . . , k ’ 1, k ≥ 1. (4.43)

Then, from (4.41) it follows that

(Ap(j) )T p(k+1) = (Ap(j) )T r(k+1) , j = 0, 1, . . . , k ’ 1.

Moreover, due to (4.43) and by the assumption of of A-orthogonality we

get

(p(j) )T r(k+1) = (p(j) )T r(k) ’ ±k (p(j) )T Ap(k) = 0, j = 0, . . . , k ’ 1(4.44)

i.e., we conclude that r(k+1) is orthogonal to every vector of the space Vk =

span(p(0) , . . . , p(k’1) ). Since p(0) = r(0) , from (4.41) it follows that Vk is

also equal to span(r(0) , . . . , r(k’1) ). Then, (4.41) implies that Ap(j) ∈ Vj+1

and thus, due to (4.44)

(Ap(j) )T r(k+1) = 0, j = 0, 1, . . . , k ’ 1.

As a consequence, (4.42) holds for j = 0, . . . , k.

The conjugate gradient method (CG) is the method obtained by choosing

the descent directions p(k) given by (4.41) and the acceleration parameter

±k as in (4.39). As a consequence, setting r(0) = b ’ Ax(0) and p(0) = r(0) ,

the k-th iteration of the conjugate gradient method takes the following

4.3 Stationary and Nonstationary Iterative Methods 153

1.4

1.2

1

0.8

G

0.6

CG

0.4

0.2

0

’0.2

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

FIGURE 4.7. Descent directions for the conjugate gradient method (denoted by

CG, dashed line) and the gradient method (denoted by G, solid line). Notice that

the CG method reaches the solution after two iterations

form

T

p(k) r(k)

±k = T

p(k) Ap(k)

x(k+1) = x(k) + ±k p(k)

r(k+1) = r(k) ’ ±k Ap(k)

(Ap(k) )T r(k+1)

βk =

(Ap(k) )T p(k)

p(k+1) = r(k+1) ’ βk p(k) .

It can also be shown (see Exercise 13) that the two parameters ±k and βk

may be alternatively expressed as

r(k) 2

r(k+1) 2

2 2

±k = , βk = . (4.45)

(k) 2

T r

p(k) Ap(k) 2

We ¬nally notice that, eliminating the descent directions from r(k+1) =

r(k) ’ ±k Ap(k) , the following recursive three-terms relation is obtained for

the residuals (see Exercise 14)

1 (k+1) 1 βk’1 βk (k’1)

Ar(k) = ’ ’ r(k) +

r + r . (4.46)

±k ±k ±k’1 ±k’1

As for the convergence of the CG method, we have the following results.

Theorem 4.11 Let A be a symmetric and positive de¬nite matrix. Any

method which employs conjugate directions to solve (3.2) terminates after

at most n steps, yielding the exact solution.

154 4. Iterative Methods for Solving Linear Systems

Proof. The directions p(0) , p(1) , . . . , p(n’1) form an A-orthogonal basis in Rn .

Moreover, since x(k) is optimal with respect to all the directions p(j) , j =

0, . . . , k ’ 1, it follows that r(k) is orthogonal to the space Sk’1 = span(p(0) , p(1) ,

. . . , p(k’1) ). As a consequence, r(n) ⊥ Sn’1 = Rn and thus r(n) = 0 which

3

implies x(n) = x.

Theorem 4.12 Let A be a symmetric and positive de¬nite matrix and

let »1 , »n be its maximum and minimum eigenvalues, respectively. The

conjugate gradient method for solving (3.2) converges after at most n steps.

Moreover, the error e(k) at the k-th iteration (with k < n) is orthogonal to

p(j) , for j = 0, . . . , k ’ 1 and

K2 (A) ’ 1

2ck

¤ e(0)

e(k) A, with c = . (4.47)

A 2k

1+c K2 (A) + 1

Proof. The convergence of the CG method in n steps is a consequence of The-

orem 4.11.

Let us prove the error estimate, assuming for simplicity that x(0) = 0. Notice

¬rst that, for ¬xed k

k

(k+1)

γj Aj b,

x =

j=0

for suitable γj ∈ R. Moreover, by construction, x(k+1) is the vector which min-

imizes the A-norm of the error at step k + 1, among all vectors of the form

z = k δj Aj b = pk (A)b, where pk (ξ) = k δj ξ j is a polynomial of degree

j=0 j=0

k and pk (A) denotes the corresponding matrix polynomial. As a consequence

¤ (x ’ z)T A(x ’ z) = xT qk+1 (A)Aqk+1 (A)x,

e(k+1) 2

(4.48)

A

where qk+1 (ξ) = 1 ’ pk (ξ)ξ ∈ P0,1 , being P0,1 = {q ∈ Pk+1 : q(0) = 1} and

k+1 k+1

qk+1 (A) the associated matrix polynomial. From (4.48) we get

e(k+1) 2

xT qk+1 (A)Aqk+1 (A)x.

= min (4.49)

A

0,1

qk+1 ∈Pk+1

Since A is symmetric positive de¬nite, there exists an orthogonal matrix Q

such that A = QΛQT with Λ = diag(»1 , . . . , »n ). Noticing that qk+1 (A) =

Qqk+1 (Λ)QT , we get from (4.49)

xT Qqk+1 (Λ)QT QΛQT Qqk+1 (Λ)QT x

e(k+1) 2

= min

A 0,1

qk+1 ∈Pk+1

xT Qqk+1 (Λ)Λqk+1 (Λ)QT x

= min

0,1

qk+1 ∈Pk+1

yT diag(qk+1 (»i )»i qk+1 (»i ))y

= min

0,1

qk+1 ∈Pk+1

n

yi »i (qk+1 (»i ))2

2

= min

0,1

qk+1 ∈Pk+1 i=1

4.3 Stationary and Nonstationary Iterative Methods 155

having set y = Qx. Thus, we can conclude that

n

¤

e(k+1) 2 2 2

min max (qk+1 (»i )) yi »i .

A

0,1

qk+1 ∈Pk+1 »i ∈σ(A) i=1

n

yi »i = e(0)

2 2

Recalling that A, we have

i=1

e(k+1) A

¤ max |qk+1 (»i )|.

min

e(0) A 0,1

qk+1 ∈Pk+1 »i ∈σ(A)

Let us now recall the following property

max |q(z)| over the space

Property 4.6 The problem of minimizing

»n ¤z¤»1

P0,1 ([»n , »1 ]) admits a unique solution, given by the polynomial

k+1

»1 + »n ’ 2ξ

ξ ∈ [»n , »1 ],

pk+1 (ξ) = Tk+1 /Ck+1 ,

»1 ’ » n

»1 +»

where Ck+1 = Tk+1 ( »1 ’»n ) and Tk+1 is the Chebyshev polynomial of degree k + 1

n

(see Section 10.10). The value of the minimum is 1/Ck+1 .

Using this property we get

e(k+1) A 1

¤

e(0) A »1 + »n

Tk+1

»1 ’ »n

from which the thesis follows since in the case of a symmetric positive de¬nite

matrix

2ck+1

1

= .

1 + c2(k+1)

Ck+1

3

The generic k-th iteration of the conjugate gradient method is well de¬ned

only if the descent direction p(k) is non null. Besides, if p(k) = 0, then the

iterate x(k) must necessarily coincide with the solution x of the system.

Moreover, irrespectively of the choice of the parameters βk , one can show

(see [Axe94], p. 463) that the sequence x(k) generated by the CG method

is such that either x(k) = x, p(k) = 0, ±k = 0 for any k, or there must exist

an integer m such that x(m) = x, where x(k) = x, p(k) = 0 and ±k = 0 for

k = 0, 1, . . . , m ’ 1.

The particular choice made for βk in (4.45) ensures that m ¤ n. In ab-

sence of rounding errors, the CG method can thus be regarded as being a

direct method, since it terminates after a ¬nite number of steps. However,

for matrices of large size, it is usually employed as an iterative scheme,

156 4. Iterative Methods for Solving Linear Systems

where the iterations are stopped when the error gets below a ¬xed toler-

ance. In this respect, the dependence of the error reduction factor on the

condition number of the matrix is more favorable than for the gradient

method. We also notice that estimate (4.47) is often overly pessimistic and

does not account for the fact that in this method, unlike what happens for

the gradient method, the convergence is in¬‚uenced by the whole spectrum

of A, and not only by its extremal eigenvalues.

Remark 4.3 (E¬ect of rounding errors) The termination property of

the CG method is rigorously valid only in exact arithmetic. The cumulating

rounding errors prevent the descent directions from being A-conjugate and

can even generate null denominators in the computation of coe¬cients ±k

and βk . This latter phenomenon, known as breakdown, can be avoided by

introducing suitable stabilization procedures; in such an event, we speak

about stabilized gradient methods.

Despite the use of these strategies, it may happen that the CG method

fails to converge (in ¬nite arithmetic) after n iterations. In such a case,

the only reasonable possibility is to restart the iterative process, taking

as residual the last computed one. By so doing, the cyclic CG method or

CG method with restart is obtained, for which, however, the convergence

properties of the original CG method are no longer valid.

4.3.5 The Preconditioned Conjugate Gradient Method

If P is a symmetric and positive de¬nite preconditioning matrix, the pre-

conditioned conjugate gradient method (PCG) consists of applying the CG

method to the preconditioned system

P’1/2 AP’1/2 y = P’1/2 b, with y = P1/2 x.

In practice, the method is implemented without explicitly requiring the

computation of P1/2 or P’1/2 . After some algebra, the following scheme is

obtained:

given x(0) and setting r(0) = b ’ Ax(0) , z(0) = P’1 r(0) e p(0) = z(0) , the

k-th iteration reads

4.3 Stationary and Nonstationary Iterative Methods 157

T

p(k) r(k)

±k = T

p(k) Ap(k)

x(k+1) = x(k) + ±k p(k)

r(k+1) = r(k) ’ ±k Ap(k)

Pz(k+1) = r(k+1)

(Ap(k) )T z(k+1)

βk =

(Ap(k) )T p(k)

p(k+1) = z(k+1) ’ βk p(k) .

The computational cost is increased with respect to the CG method, as

one needs to solve at each step the linear system Pz(k+1) = r(k+1) . For this

system the symmetric preconditioners examined in Section 4.3.2 can be

used. The error estimate is the same as for the nonpreconditioned method,

provided to replace the matrix A by P’1 A.

In Program 20 an implementation of the PCG method is reported. For

a description of the input/output parameters, see Program 19.

Program 20 - conjgrad : Preconditioned conjugate gradient method

function [x, error, niter, ¬‚ag] = conjgrad(A, x, b, P, maxit, tol)

¬‚ag = 0; niter = 0; bnrm2 = norm( b );

if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

r = b - A*x; error = norm( r ) / bnrm2;

if ( error < tol ) return, end

for niter = 1:maxit

z = P \ r; rho = (r™*z);

if niter > 1

beta = rho / rho1; p = z + beta*p;

else

p = z;

end

q = A*p; alpha = rho / (p™*q );

x = x + alpha * p; r = r - alpha*q;

error = norm( r ) / bnrm2;

if ( error <= tol ), break, end

rho1 = rho;

end

if ( error > tol ) ¬‚ag = 1; end

158 4. Iterative Methods for Solving Linear Systems

Example 4.7 Let us consider again the linear system of Example 4.6. The CG

method has been run with the same input data as in the previous example. It

converges in 3 iterations for m = 16 and in 45 iterations for m = 400. Using the

same preconditioner as in Example 4.6, the number of iterations decreases from

•

45 to 26, in the case m = 400.

0

10

’2

10

’4

10

’6

10

’8

10

’10

10

’12

10

’14

10

0 5 10 15 20 25 30 35 40 45

FIGURE 4.8. Behavior of the residual, normalized to the right-hand side, as a

function of the number of iterations for the conjugate gradient method applied

to the systems of Example 4.6 in the case m = 400. The curve in dashed line

refers to the non preconditioned method, while the curve in solid line refers to

the preconditioned one

4.3.6 The Alternating-Direction Method

Assume that A = A1 +A2 , with A1 and A2 symmetric and positive de¬nite.

The alternating direction method (ADI), as introduced by Peaceman and

Rachford [PJ55], is an iterative scheme for (3.2) which consists of solving

the following systems ∀k ≥ 0

(I + ±1 A1 )x(k+1/2) = (I ’ ±1 A2 )x(k) + ±1 b,

(4.50)

= (I ’ ±2 A1 )x

(k+1) (k+1/2)

(I + ±2 A2 )x + ±2 b

where ±1 and ±2 are two real parameters. The ADI method can be cast in

the form (4.2) setting

B = (I + ±2 A2 )’1 (I ’ ±2 A1 )(I + ±1 A1 )’1 (I ’ ±1 A2 ),

f = ±1 (I ’ ±2 A1 )(I + ±1 A1 )’1 + ±2 I b.

Both B and f depend on ±1 and ±2 . The following estimate holds

(1) (2)

1 ’ ±2 »i 1 ’ ±1 »i

ρ(B) ¤ max max ,

(1) (2)

i=1,... ,n ±1 »i i=1,... ,n

1+ 1+ ±2 »i

4.4 Methods Based on Krylov Subspace Iterations 159

(i) (i)

where »1 and »2 , for i = 1, . . . , n, are the eigenvalues of A1 and A2 ,

respectively. The method converges if ρ(B) < 1, which is always veri¬ed if

(j)

±1 = ±2 = ± > 0. Moreover (see [Axe94]) if γ ¤ »i ¤ δ ∀i = 1, . . . , n,

∀j = 1, 2, for suitable γ and δ then the ADI method converges with the

√

choice ±1 = ±2 = 1/ δγ, provided that γ/δ tends to 0 as the size of A

grows. In such an event the corresponding spectral radius satis¬es

2

1’ γ/δ

ρ(B) ¤ .

1+ γ/δ

4.4 Methods Based on Krylov Subspace Iterations

In this section we introduce iterative methods based on Krylov subspace

iterations. For the proofs and further analysis, we refer to [Saa96], [Axe94]

and [Hac94].

Consider the Richardson method (4.24) with P=I; the residual at the

k-th step can be related to the initial residual as

k’1

(I ’ ±j A)r(0)

(k)

r = (4.51)

j=0

so that r(k) = pk (A)r(0) , where pk (A) is a polynomial in A of degree k. If

we introduce the space

Km (A; v) = span v, Av, . . . , Am’1 v , (4.52)

it immediately appears from (4.51) that r(k) ∈ Kk+1 (A; r(0) ). The space

de¬ned in (4.52) is called the Krylov subspace of order m. It is a subspace

of the set spanned by all the vectors u ∈ Rn that can be written as u =

pm’1 (A)v, where pm’1 is a polynomial in A of degree ¤ m ’ 1.

In an analogous manner as for (4.51), it is seen that the iterate x(k) of

the Richardson method is given by

k’1

(k) (0)

±j r(j)

x =x +

j=0

so that x(k) belongs to the following space

Wk = v = x(0) + y, y ∈ Kk (A; r(0) ) . (4.53)

k’1

Notice also that j=0 ±j r(j) is a polynomial in A of degree less than k ’ 1.

In the non preconditioned Richardson method we are thus looking for an

160 4. Iterative Methods for Solving Linear Systems

approximate solution to x in the space Wk . More generally, we can think

of devising methods that search for approximate solutions of the form

x(k) = x(0) + qk’1 (A)r(0) , (4.54)

where qk’1 is a polynomial selected in such a way that x(k) be, in a sense

that must be made precise, the best approximation of x in Wk . A method

that looks for a solution of the form (4.54) with Wk de¬ned as in (4.53) is

called a Krylov method.

A ¬rst question concerning Krylov subspace iterations is whether the

dimension of Km (A; v) increases as the order m grows. A partial answer is

provided by the following result.

Property 4.7 Let A ∈ Rn—n and v ∈ Rn . The Krylov subspace Km (A; v)

has dimension equal to m i¬ the degree of v with respect to A, denoted by

degA (v), is not less than m, where the degree of v is de¬ned as the minimum

degree of a monic non null polynomial p in A, for which p(A)v = 0.

The dimension of Km (A; v) is thus equal to the minimum between m and

the degree of v with respect to A and, as a consequence, the dimension

of the Krylov subspaces is certainly a nondecreasing function of m. Notice

that the degree of v cannot be greater than n due to the Cayley-Hamilton

Theorem (see Section 1.7).

Example 4.8 Consider the matrix A = tridiag4 (’1, 2, ’1). The vector v =

(1, 1, 1, 1)T has degree 2 with respect to A since p2 (A)v = 0 with p2 (A) = I4 ’

3A + A2 , while there is no monic polynomial p1 of degree 1 for which p1 (A)v = 0.

As a consequence, all Krylov subspaces from K2 (A; v) on, have dimension equal

to 2. The vector w = (1, 1, ’1, 1)T has, instead, degree 4 with respect to A. •

For a ¬xed m, it is possible to compute an orthonormal basis for Km (A; v)

using the so-called Arnoldi algorithm.

Setting v1 = v/ v 2 , this method generates an orthonormal basis {vi }

for Km (A; v1 ) using the Gram-Schmidt procedure (see Section 3.4.3). For

k = 1, . . . , m, the Arnoldi algorithm computes

T

hik = vi Avk , i = 1, 2, . . . , k,

(4.55)

k

wk = Avk ’ hik vi , hk+1,k = wk 2.

i=1

If wk = 0 the process terminates and in such a case we say that a breakdown

of the algorithm has occurred; otherwise, we set vk+1 = wk / wk 2 and the

algorithm restarts, incrementing k by 1.

4.4 Methods Based on Krylov Subspace Iterations 161

It can be shown that if the method terminates at the step m then the

vectors v1 , . . . , vm form a basis for Km (A; v). In such a case, if we denote

by Vm ∈ Rn—m the matrix whose columns are the vectors vi , we have

T T (4.56)

Vm AVm = Hm , Vm+1 AVm = Hm ,

where Hm ∈ R(m+1)—m is the upper Hessenberg matrix whose entries hij

are given by (4.55) and Hm ∈ Rm—m is the restriction of Hm to the ¬rst m

rows and m columns.

The algorithm terminates at an intermediate step k < m i¬ degA (v1 ) =

k. As for the stability of the procedure, all the considerations valid for the

Gram-Schmidt method hold. For more e¬cient and stable computational

variants of (4.55), we refer to [Saa96].

The functions arnoldi alg and GSarnoldi, invoked by Program 21, pro-

vide an implementation of the Arnoldi algorithm. In output, the columns

of V contain the vectors of the generated basis, while the matrix H stores

the coe¬cients hik computed by the algorithm. If m steps are carried out,

V = Vm and H(1 : m, 1 : m) = Hm .

Program 21 - arnoldi alg : The Arnoldi algorithm

function [V,H]=arnoldi alg(A,v,m)

v=v/norm(v,2); V=[v1]; H=[]; k=0;

while k <= m-1

[k,V,H] = GSarnoldi(A,m,k,V,H);

end

function [k,V,H]=GSarnoldi(A,m,k,V,H)

k=k+1; H=[H,V(:,1:k)™*A*V(:,k)];

s=0; for i=1:k, s=s+H(i,k)*V(:,i); end

w=A*V(:,k)-s; H(k+1,k)=norm(w,2);

if ( H(k+1,k) <= eps ) & ( k < m )

V=[V,w/H(k+1,k)];

else

k=m+1;

end

Having introduced an algorithm for generating the basis for a Krylov sub-

space of any order, we can now solve the linear system (3.2) by a Krylov

method. As already noticed, for all of these methods the iterate x(k) is

always of the form (4.54) and, for a given r(0) , the vector x(k) is selected as

being the unique element in Wk which satis¬es a criterion of minimal dis-

tance from x. Thus, the feature distinguishing two di¬erent Krylov methods

is the criterion for selecting x(k) .

The most natural idea consists of searching for x(k) ∈ Wk as the vector

which minimizes the Euclidean norm of the error. This approach, how-

162 4. Iterative Methods for Solving Linear Systems

ever, does not work in practice since x(k) would depend on the (unknown)

solution x.

Two alternative strategies can be pursued:

1. compute x(k) ∈ Wk enforcing that the residual r(k) is orthogonal to

any vector in Kk (A; r(0) ), i.e., we look for x(k) ∈ Wk such that

vT (b ’ Ax(k) ) = 0 ∀v ∈ Kk (A; r(0) ); (4.57)

2. compute x(k) ∈ Wk minimizing the Euclidean norm of the residual

r(k) 2 , i.e.

b ’ Ax(k) = min b ’ Av 2 . (4.58)

2

v∈Wk

Satisfying (4.57) leads to the Arnoldi method for linear systems (more

commonly known as FOM, full orthogonalization method), while satisfying

(4.58) yields the GMRES (generalized minimum residual) method.

In the two forthcoming sections we shall assume that k steps of the

Arnoldi algorithm have been carried out, in such a way that an orthonormal

basis for Kk (A; r(0) ) has been generated and stored into the column vectors

of the matrix Vk with v1 = r(0) / r(0) 2 . In such a case the new iterate x(k)

can always be written as

x(k) = x(0) + Vk z(k) , (4.59)

where z(k) must be selected according to a ¬xed criterion.

4.4.1 The Arnoldi Method for Linear Systems

Let us enforce that r(k) be orthogonal to Kk (A; r(0) ) by requiring that

(4.57) holds for all the basis vectors vi , i.e.

Vk r(k) = 0.

T

(4.60)

Since r(k) = b ’ Ax(k) with x(k) of the form (4.59), relation (4.60) becomes

Vk (b ’ Ax(0) ) ’ Vk AVk z(k) = Vk r(0) ’ Vk AVk z(k) = 0.

T T T T

(4.61)

Due to the orthonormality of the basis and the choice of v1 , Vk r(0) = T

r(0) 2 e1 , e1 being the ¬rst unit vector of Rm . Recalling (4.56), from (4.61)

it turns out that z(k) is the solution to the linear system

Hk z(k) = r(0) 2 e1 . (4.62)

Once z(k) is known, we can compute x(k) from (4.59). Since Hk is an upper

Hessenberg matrix, the linear system in (4.62) can be easily solved, for

instance, resorting to the LU factorization of Hk .

4.4 Methods Based on Krylov Subspace Iterations 163

We notice that the method, if working in exact arithmetic, cannot execute

more than n steps and that it terminates after m < n steps only if a

breakdown in the Arnoldi algorithm occurs. As for the convergence of the

method, the following result holds.

Theorem 4.13 In exact arithmetic the Arnoldi method yields the solution

of (3.2) after at most n iterations.

Proof. If the method terminates at the n-th iteration, then it must necessarily

be x(n) = x since Kn (A; r(0) ) = Rn . Conversely, if a breakdown occurs after m

iterations, for a suitable m < n, then x(m) = x. Indeed, inverting the ¬rst relation

in (4.56), we get

x(m) = x(0) + Vm z(m) = x(0) + Vm H’1 Vm r(0) = A’1 b.

T

m

3

In its naive form, FOM does not require an explicit computation of the

solution or the residual, unless a breakdown occurs. Therefore, monitoring

its convergence (by computing, for instance, the residual at each step) might

be computationally expensive. The residual, however, is available without

explicitly requiring to compute the solution since at the k-th step we have

b ’ Ax(k) = hk+1,k |eT zk |

2 k

and, as a consequence, one can decide to stop the method if

hk+1,k |eT zk |/ r(0) ¤µ (4.63)

2

k

µ > 0 being a ¬xed tolerance.

The most relevant consequence of Theorem 4.13 is that FOM can be

regarded as a direct method, since it yields the exact solution after a ¬nite

number of steps. However, this fails to hold when working in ¬‚oating point

arithmetic due to the cumulating rounding errors. Moreover, if we also

account for the high computational e¬ort, which, for a number of m steps

and a sparse matrix of order n with nz nonzero entries, is of the order of

2(nz + mn) ¬‚ops, and the large memory occupation needed to store the

matrix Vm , we conclude that the Arnoldi method cannot be used in the

practice, except for small values of m.

Several remedies to this drawback are available, one of which consisting

of preconditioning the system (using, for instance, one of the precondition-

ers proposed in Section 4.3.2). Alternatively, we can also introduce some

modi¬ed versions of the Arnoldi method following two approaches:

1. no more than m consecutive steps of FOM are taken, m being a small

¬xed number (usually, m 10). If the method fails to converge, we set

164 4. Iterative Methods for Solving Linear Systems

x(0) = x(m) and FOM is repeated for other m steps. This procedure

is carried out until convergence is achieved. This method, known as

FOM(m) or FOM with restart, reduces the memory occupation, only

requiring to store matrices with m columns at most;

2. a limitation is set on the number of directions involved in the orthog-

onalization procedure in the Arnoldi algorithm, yielding the incom-

plete orthogonalization method or IOM. In the practice, the k-th step

of the Arnoldi algorithm generates a vector vk+1 which is orthonor-

mal, at most, to the q preceding vectors, where q is ¬xed according

to the amount of available memory.

It is worth noticing that Theorem 4.13 does no longer hold for the methods

stemming from the two strategies above.

Program 22 provides an implementation of the FOM algorithm with a

stopping criterion based on the residual (4.63). The input parameter m is

the maximum admissible size of the Krylov subspace that is being gener-

ated and represents, as a consequence, the maximum admissible number of

iterations.

Program 22 - arnoldi met : The Arnoldi method for linear systems

function [x,k]=arnoldi met(A,b,m,x0,toll)

r0=b-A*x0; nr0=norm(r0,2);

if nr0 ˜= 0

v1=r0/nr0; V=[v1]; H=[]; k=0; istop=0;

while (k <= m-1) & (istop == 0)

[k,V,H] = GSarnoldi(A,m,k,V,H);

[nr,nc]=size(H); e1=eye(nc);

y=(e1(:,1)™*nr0)/H(1:nc,:);

residual = H(nr,nc)*abs(y*e1(:,nc));

if residual <= toll

istop = 1; y=y™;

end

end

if istop==0

[nr,nc]=size(H); e1=eye(nc);

y=(e1(:,1)™*nr0)/H(1:nc,:); y=y™;

end

x=x0+V(:,1:nc)*y;

else

x=x0;

end

Example 4.9 Let us solve the linear system Ax = b with A = tridiag100 (’1, 2,

’1) and b such that the solution is x = 1T . The initial vector is x(0) = 0T

and toll=10’10 . The method converges in 50 iterations and Figure 4.9 reports

its convergence history. Notice the sudden, dramatic, reduction of the residual,

4.4 Methods Based on Krylov Subspace Iterations 165

which is a typical warning that the last generated subspace Wk is su¬ciently rich

•

to contain the exact solution of the system.

2

10

0

10

’2

10

’4

10

’6

10

’8

10

’10

10

’12

10

’14

10

’16

10

0 10 20 30 40 50 60

FIGURE 4.9. The behavior of the residual as a function of the number of itera-

tions for the Arnoldi method applied to the linear system in Example 4.9

4.4.2 The GMRES Method

This method is characterized by selecting x(k) in such a way to minimize

the Euclidean norm of the residual at each k-th step. Recalling (4.59) we

have

r(k) = r(0) ’ AVk z(k) , (4.64)

but, since r(0) = v1 r(0) and (4.56) holds, relation (4.64) becomes

2

r(k) = Vk+1 ( r(0) 2 e1 ’ Hk z(k) ), (4.65)

where e1 is the ¬rst unit vector of Rk+1 . Therefore, in the GMRES method

the solution at step k can be computed through (4.59) as

r(0) 2 e1 ’ Hk z(k)

z(k) chosen in such a way to minimize (4.66)

2

(the matrix Vk+1 appearing in (4.65) does not change the value of · 2

since it is orthogonal). Having to solve at each step a least-squares problem

of size k, the GMRES method will be the more e¬ective the smaller is

the number of iterations. Exactly as for the Arnoldi method, the GMRES

method terminates at most after n iterations, yielding the exact solution.

Premature stops are due to a breakdown in the orthonormalization Arnoldi

algorithm. More precisely, we have the following result.

Property 4.8 A breakdown occurs for the GMRES method at a step m

(with m < n) i¬ the computed solution x(m) coincides with the exact solu-

tion to the system.

166 4. Iterative Methods for Solving Linear Systems

A basic implementation of the GMRES method is provided in Program 23.

This latter requires in input the maximum admissible size m for the Krylov

subspace and the tolerance toll on the Euclidean norm of the residual

normalized to the initial residual. This implementation of the method com-

putes the solution x(k) at each step in order to evaluate the residual, with

a consequent increase of the computational e¬ort.

Program 23 - GMRES : The GMRES method for linear systems

function [x,k]=gmres(A,b,m,toll,x0)

r0=b-A*x0; nr0=norm(r0,2);

if nr0 ˜= 0

v1=r0/nr0; V=[v1]; H=[]; k=0; residual=1;

while k <= m-1 & residual > toll,

[k,V,H] = GSarnoldi(A,m,k,V,H);

[nr,nc]=size(H); y=(H™*H) \ (H™*nr0*[1;zeros(nr-1,1)]);

x=x0+V(:,1:nc)*y; residual = norm(b-A*x,2)/nr0;

end

else

x=x0;

end

To improve the e¬ciency of the GMRES algorithm it is necessary to devise

a stopping criterion which does not require the explicit evaluation of the

residual at each step. This is possible, provided that the linear system with

upper Hessenberg matrix Hk is appropriately solved.

In practice, Hk is transformed into an upper triangular matrix Rk ∈

R (k+1)—k

with rk+1,k = 0 such that QT Rk = Hk , where Qk is a matrix

k

obtained as the product of k Givens rotations (see Section 5.6.3). Then,

since Qk is orthogonal, it can be seen that minimizing r(0) 2 e1 ’ Hk z(k) 2

is equivalent to minimize fk ’ Rk z(k) 2 , with fk = Qk r(0) 2 e1 . It can

also be shown that the k + 1-th component of fk is, in absolute value, the

Euclidean norm of the residual at the k-th step.

As FOM, the GMRES method entails a high computational e¬ort and

a large amount of memory, unless convergence occurs after few iterations.

For this reason, two variants of the algorithm are available, one named

GMRES(m) and based on the restart after m steps, the other named Quasi-

GMRES or QGMRES and based on stopping the Arnoldi orthogonalization

process. It is worth noting that these two methods do not enjoy Property

4.8.

Remark 4.4 (Projection methods) Denoting by Yk and Lk two generic

m-dimensional subspaces of Rn , we call projection method a process which

generates an approximate solution x(k) at step k, enforcing that x(k) ∈ Yk

4.4 Methods Based on Krylov Subspace Iterations 167

and that the residual r(k) = b ’ Ax(k) be orthogonal to Lk . If Yk = Lk , the

projection process is said to be orthogonal, oblique otherwise (see [Saa96]).

The Krylov subspace iterations can be regarded as being projection

methods. For instance, the Arnoldi method is an orthogonal projection

method where Lk = Yk = Kk (A; r(0) ), while the GMRES method is an

oblique projection method with Yk = Kk (A; r(0) ) and Lk = AYk . It is

worth noticing that some classical methods introduced in previous sections

fall into this category. For example, the Gauss-Seidel method is an orthogo-

nal projection method where at the k-th step Kk (A; r(0) ) = span(ek ), with

k = 1, . . . , n. The projection steps are carried out cyclically from 1 to n

until convergence.

4.4.3 The Lanczos Method for Symmetric Systems

The Arnoldi algorithm simpli¬es considerably if A is symmetric since the

matrix Hm is tridiagonal and symmetric (indeed, from (4.56) it turns out

that Hm must be symmetric, so that, being upper Hessenberg by construc-

tion, it must necessarily be tridiagonal). In such an event the method is

more commonly known as the Lanczos algorithm. For ease of notation, we

henceforth let ±i = hii and βi = hi’1,i .

An implementation of the Lanczos algorithm is provided in Program 24.

Vectors alpha and beta contain the coe¬cients ±i and βi computed by the

scheme.

Program 24 - Lanczos : The Lanczos method for linear systems

function [V,alpha,beta]=lanczos(A,m)

n=size(A); V=[0*[1:n]™,[1,0*[1:n-1]]™];

beta(1)=0; normb=1; k=1;

while k <= m & normb >= eps

vk = V(:,k+1); w = A*vk-beta(k)*V(:,k);

alpha(k)= w™*vk; w = w - alpha(k)*vk

normb = norm(w,2);

if normb ˜= 0

beta(k+1)=normb; V=[V,w/normb]; k=k+1;

end

end

[n,m]=size(V); V=V(:,2:m-1);

alpha=alpha(1:n); beta=beta(2:n);

The algorithm, which is far superior to Arnoldi™s one as far as memory

saving is concerned, is not numerically stable since only the ¬rst generated

vectors are actually orthogonal. For this reason, several stable variants have

been devised.

As in previous cases, also the Lanczos algorithm can be employed as a

solver for linear systems, yielding a symmetric form of the FOM method. It

168 4. Iterative Methods for Solving Linear Systems

can be shown that r(k) = γk vk+1 , for a suitable γk (analogously to (4.63))

so that the residuals are all mutually orthogonal.

Remark 4.5 (The conjugate gradient method) If A is symmetric and

positive de¬nite, starting from the Lanczos method for linear systems it is

possible to derive the conjugate gradient method already introduced in Sec-

tion 4.3.4 (see [Saa96]). The conjugate gradient method is a variant of the

Lanczos method where the orthonormalization process remains incomplete.

As a matter of fact, the A-conjugate directions of the CG method can

be characterized as follows. If we carry out at the generic k-th step the

LU factorization Hk = Lk Uk , with Lk (Uk ) lower (upper) bidiagonal, the

iterate x(k) of the Lanczos method for systems reads

x(k) = x(0) + Pk L’1 r(0) 2 e1 ,

k

with Pk = Vk U’1 . The column vectors of Pk are mutually A-conjugate.

k

Indeed, PT APk is symmetric and bidiagonal since

k

PT APk = U’T Hk U’1 = U’T Lk ,

k k k k

so that it must necessarily be diagonal. As a result, pT Api = 0 if i = j,

j

having denoted by pi the i-th column vector of matrix Pk .

As happens for the FOM method, also the GMRES method simpli¬es

if A is symmetric. The resulting scheme is called conjugate residuals or

CR method since it enjoys the property that the residuals are mutually

A-conjugate. Variants of this method are the generalized conjugate resid-

uals method (GCR) and the method commonly known as ORTHOMIN

(obtained by truncation of the orthonormalization process as done for the

IOM method).

4.5 The Lanczos Method for Unsymmetric Systems

The Lanczos orthogonalization process can be extended to deal with un-

symmetric matrices through a bi-orthogonalization procedure as follows.

m m

Two bases, {vi }i=1 and {zi }i=1 , are generated for the subspaces Km (A; v1 )

and Km (AT ; z1 ), respectively, with zT v1 = 1, such that

1

zT vj = δij , i, j = 1, . . . , m. (4.67)

i

Two sets of vectors satisfying (4.67) are said to be bi-orthogonal and can

be obtained through the following algorithm: setting β1 = γ1 = 0 and z0 =

v0 = 0T , at the generic k-th step, with k = 1, . . . , m, we set ±k = zT Avk ,

k

then we compute

vk+1 = Avk ’ ±k vk ’ βk vk’1 , zk+1 = AT zk ’ ±k zk ’ γk zk’1 .

˜ ˜

4.5 The Lanczos Method for Unsymmetric Systems 169

|˜T vk+1 | = 0 the algorithm is stopped, otherwise we set

zk+1 ˜

If γk+1 =

βk+1 = zT vk+1 /γk+1 and generate two new vectors in the basis as

˜k+1 ˜

˜ ˜

vk+1 = vk+1 /γk+1 , zk+1 = zk+1 /βk+1 .

If the process terminates after m steps, denoting by Vm and Zm the ma-

trices whose columns are the vectors of the basis that has been generated,

we have

ZT AVm = Tm ,

m

Tm being the following tridiagonal matrix

®

0

±1 β2

γ2 ±2 . . .

Tm = .

.. ..

. . βm

° »

0 γm ±m

As in the symmetric case, the bi-orthogonalization Lanczos algorithm can

be utilized to solve the linear system (3.2). For this purpose, for m ¬xed,

m m

once the bases {vi }i=1 and {zi }i=1 have been constructed, it su¬ces to set

x(m) = x(0) + Vm y(m) ,

where y(m) is the solution to the linear system Tm y(m) = r(0) 2 e1 . It

is also possible to introduce a stopping criterion based on the residual,

without computing it explicitly, since

= |γm+1 eT ym | vm+1 2 .

r(m) 2 m

An implementation of the Lanczos method for unsymmetric systems is

given in Program 25. If a breakdown of the algorithm occurs, i.e., if γk+1 =

0, the method stops returning in output a negative value of the variable

niter which denotes the number of iterations necessary to reduce the initial

residual by a factor toll.

170 4. Iterative Methods for Solving Linear Systems

Program 25 - Lanczosnosym : The Lanczos method for unsymmetric

systems

function [xk,nres,niter]=lanczosnosym(A,b,x0,m,toll)

r0=b-A*x0; nres0=norm(r0,2);

if nres0 ˜= 0

V=r0/nres0; Z=V; gamma(1)=0; beta(1)=0; k=1; nres=1;

while k <= m & nres > toll

vk=V(:,k); zk=Z(:,k);

if k==1, vk1=0*vk; zk1=0*zk;

else, vk1=V(:,k-1); zk1=Z(:,k-1); end

alpha(k)=zk™*A*vk;

tildev=A*vk-alpha(k)*vk-beta(k)*vk1;

tildez=A™*zk-alpha(k)*zk-gamma(k)*zk1;

gamma(k+1)=sqrt(abs(tildez™*tildev));

if gamma(k+1) == 0, k=m+2;

else

beta(k+1)=tildez™*tildev/gamma(k+1);

Z=[Z,tildez/beta(k+1)];

V=[V,tildev/gamma(k+1)];

end

if k˜=m+2

if k==1

Tk = alpha;

else

Tk=diag(alpha)+diag(beta(2:k),1)+diag(gamma(2:k),-1);

end

yk=Tk \ (nres0*[1,0*[1:k-1]]™);

xk=x0+V(:,1:k)*yk;

nres=abs(gamma(k+1)*[0*[1:k-1],1]*yk)*norm(V(:,k+1),2)/nres0;

k=k+1;

end

end

else

x=x0;

end

if k==m+2, niter=-k; else, niter=k-1; end

Example 4.10 Let us solve the linear system with matrix A = tridiag100 (’0.5, 2,

’1) and right-side b selected in such a way that the exact solution is x = 1T .

Using Program 25 with toll= 10’13 and a randomly generated x0, the algorithm

converges in 59 iterations. Figure 4.10 shows the convergence history reporting

•

the graph of r(k) 2 / r(0) 2 as a function of the number of iterations.

We conclude recalling that some variants of the unsymmetric Lanczos

method have been devised, that are characterized by a reduced compu-

tational cost. We refer the interested reader to the bibliography below for a

4.6 Stopping Criteria 171

0

10

’2

10

’4

10

’6

10

’8

10

’10

10

’12

10

’14

10

0 10 20 30 40 50 60

FIGURE 4.10. Graph of the residual normalized to the initial residual as a func-

tion of the number of iterations for the Lanczos method applied to the system in

Example 4.10

complete description of the algorithms and to the programs included in the

MATLAB version of the public domain library templates for their e¬cient

implementation [BBC+ 94].

1. The bi-conjugate gradient method (BiCG): it can be derived by the

unsymmetric Lanczos method in the same way as the conjugate gra-

dient method is obtained from the FOM method [Fle75];

2. the Quasi-Minimal Residual method (QMR): it is analogous to the

GMRES method, the only di¬erence being the fact that the Arnoldi

orthonormalization process is replaced by the Lanczos bi-orthogona-

lization;

3. the conjugate gradient squared method (CGS): the matrix-vector prod-

ucts involving the transposed matrix AT are removed. A variant of

this method, known as BiCGStab, is characterized by a more reg-

ular convergence than provided by the CGS method (see [Son89],

[vdV92]).

4.6 Stopping Criteria

In this section we address the problem of how to estimate the error intro-

duced by an iterative method and the number kmin of iterations needed to

reduce the initial error by a factor µ.

In practice, kmin can be obtained by estimating the convergence rate of

(4.2), i.e. the rate at which e(k) ’ 0 as k tends to in¬nity. From (4.4),

we get

e(k)

¤ Bk ,

(0)

e

172 4. Iterative Methods for Solving Linear Systems

so that Bk is an estimate of the reducing factor of the norm of the error

after k steps. Typically, the iterative process is continued until e(k) has

reduced with respect to e(0) by a certain factor µ < 1, that is

e(k) ¤ µ e(0) . (4.68)

If we assume that ρ(B) < 1, then Property 1.13 implies that there exists

a suitable matrix norm · such that B < 1. As a consequence, Bk

tends to zero as k tends to in¬nity, so that (4.68) can be satis¬ed for a

su¬ciently large k such that Bk ¤ µ holds. However, since Bk < 1, the

previous inequality amounts to requiring that

1

k ≥ log(µ)/ = ’ log(µ)/Rk (B),

log Bk (4.69)

k

where Rk (B) is the average convergence rate introduced in De¬nition 4.2.

From a practical standpoint, (4.69) is useless, being nonlinear in k; if, how-

ever, the asymptotic convergence rate is adopted, instead of the average

one, the following estimate for kmin is obtained

’ log(µ)/R(B).

kmin (4.70)

This latter estimate is usually rather optimistic, as con¬rmed by Example

4.11.

Example 4.11 For the matrix A3 of Example 4.2, in the case of Jacobi method,

letting µ = 10’5 , condition (4.69) is satis¬ed with kmin = 16, while (4.70) yields

kmin = 15, with a good agreement between the two estimates. Instead, on the

matrix A4 of Example 4.2, we ¬nd that (4.69) is satis¬ed with kmin = 30, while

•

(4.70) yields kmin = 26.

4.6.1 A Stopping Test Based on the Increment

From the recursive error relation e(k+1) = Be(k) , we get

e(k+1) ¤ B e(k) . (4.71)

Using the triangular inequality we get

e(k+1) ¤ B ( e(k+1) + x(k+1) ’ x(k) ),

from which it follows that

B

x ’ x(k+1) ¤ x(k+1) ’ x(k) . (4.72)

1’ B

In particular, taking k = 0 in (4.72) and applying recursively (4.71) we also

get

B k+1 (1)

x ’ x(k+1) ¤ x ’ x(0) ,

1’ B

4.6 Stopping Criteria 173

which can be used to estimate the number of iterations necessary to ful¬ll

the condition e(k+1) ¤ µ, for a given tolerance µ.

In the practice, B can be estimated as follows: since

x(k+1) ’ x(k) = ’(x ’ x(k+1) ) + (x ’ x(k) ) = B(x(k) ’ x(k’1) ),

a lower bound of B is provided by c = δk+1 /δk , where δj+1 = x(j+1) ’

x(j) , with j = k ’ 1, k. Replacing B by c, the right-hand side of (4.72)

suggests using the following indicator for e(k+1)

2

δk+1

(k+1)

= . (4.73)

δk ’ δk+1

Due to the kind of approximation of B that has been used, the reader is

warned that (k+1) should not be regarded as an upper bound for e(k+1) .

However, often (k+1) provides a reasonable indication about the true error

behavior, as we can see in the following example.

Example 4.12 Consider the linear system Ax=b with

® ®

4 1 1 6

°2 0 », b = ° ’7 » ,

’9

A=

’8 ’6 ’14

0

which admits the unit vector as exact solution. Let us apply the Jacobi method

and estimate the error at each step by using (4.73). Figure 4.11 shows an ac-

ceptable agreement between the behavior of the error e(k+1) ∞ and that of its

•

estimate (k+1) .

1

10

0

10

’1

10

’2

10

’3

10

’4

10