<<

. 7
( 26)



>>

nonstationary Richardson method employing (4.37) to evaluate the acceler-
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

<<

. 7
( 26)



>>