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