<< . .

. 25
( : 95)

. . >>

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)

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.

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.

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.


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)

µ > 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

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);
residual = H(nr,nc)*abs(y*e1(:,nc));
if residual <= toll
istop = 1; y=y™;
if istop==0
[nr,nc]=size(H); e1=eye(nc);
y=(e1(:,1)™*nr0)/H(1:nc,:); y=y™;

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.










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
r(k) = r(0) ’ AVk z(k) , (4.64)
but, since r(0) = v1 r(0) and (4.56) holds, relation (4.64) becomes

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)

(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;

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

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

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;
[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 ,

with Pk = Vk U’1 . The column vectors of Pk are mutually A-conjugate.
Indeed, PT APk is symmetric and bidiagonal since

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,
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

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

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 ,
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 ,

Tm being the following tridiagonal matrix
® 
±1 β2
 

<< . .

. 25
( : 95)

. . >>