j=0

The ¬rst term represents the error that is made by the iterative method

in exact arithmetic; if the method is convergent, this error is negligible for

su¬ciently large values of k. The second term refers instead to rounding

error propagation; its analysis is quite technical and is carried out, for

instance, in [Hig88] in the case of Jacobi, Gauss-Seidel and SOR methods.

4.2 Linear Iterative Methods 133

4.2.5 Block Matrices

The methods of the previous sections are also referred to as point (or line)

iterative methods, since they act on single entries of matrix A. It is possible

to devise block versions of the algorithms, provided that D denotes the block

diagonal matrix whose entries are the m — m diagonal blocks of matrix A

(see Section 1.6).

The block Jacobi method is obtained taking again P=D and N=D-A. The

method is well-de¬ned only if the diagonal blocks of D are nonsingular. If

A is decomposed in p — p square blocks, the block Jacobi method is

p

(k+1) (k)

= bi ’

Aii xi Aij xj , i = 1, . . . , p,

j=1

j=i

having also decomposed the solution vector and the right side in blocks of

size p, denoted by xi and bi , respectively. As a result, at each step, the block

Jacobi method requires solving p linear systems of matrices Aii . Theorem

4.3 is still valid, provided that D is substituted by the corresponding block

diagonal matrix.

In a similar manner, the block Gauss-Seidel and block SOR methods can

be introduced.

4.2.6 Symmetric Form of the Gauss-Seidel and SOR Methods

Even if A is a symmetric matrix, the Gauss-Seidel and SOR methods gen-

erate iteration matrices that are not necessarily symmetric. For that, we

introduce in this section a technique that allows for symmetrizing these

schemes. The ¬nal aim is to provide an approach for generating symmetric

preconditioners (see Section 4.3.2).

Firstly, let us remark that an analogue of the Gauss-Seidel method can

be constructed, by simply exchanging E with F. The following iteration

can thus be de¬ned, called the backward Gauss-Seidel method

(D ’ F)x(k+1) = Ex(k) + b

with iteration matrix given by BGSb = (D ’ F)’1 E.

The symmetric Gauss-Seidel method is obtained by combining an itera-

tion of Gauss-Seidel method with an iteration of backward Gauss-Seidel

method. Precisely, the k-th iteration of the symmetric Gauss-Seidel method

is

(D ’ E)x(k+1/2) = Fx(k) + b, (D ’ F)x(k+1) = Ex(k+1/2) + b.

134 4. Iterative Methods for Solving Linear Systems

Eliminating x(k+1/2) , the following scheme is obtained

x(k+1) = BSGS x(k) + bSGS ,

BSGS = (D ’ F)’1 E(D ’ E)’1 F, (4.21)

bSGS = (D ’ F)’1 [E(D ’ E)’1 + I]b.

The preconditioning matrix associated with (4.21) is

PSGS = (D ’ E)D’1 (D ’ F).

The following result can be proved (see [Hac94]).

Property 4.5 If A is a symmetric positive de¬nite matrix, the symmet-

ric Gauss-Seidel method is convergent, and, moreover, BSGS is symmetric

positive de¬nite.

In a similar manner, de¬ning the backward SOR method

(D ’ ωF)x(k+1) = [ωE + (1 ’ ω)D] x(k) + ωb,

and combining it with a step of SOR method, the following symmetric SOR

method or SSOR, is obtained

x(k+1) = Bs (ω)x(k) + bω

where

Bs (ω) = (D ’ ωF)’1 (ωE + (1 ’ ω)D)(D ’ ωE)’1 (ωF + (1 ’ ω)D),

bω = ω(2 ’ ω)(D ’ ωF)’1 D(D ’ ωE)’1 b.

The preconditioning matrix of this scheme is

1 ω 1

D’1

D’E D’F .

PSSOR (ω) = (4.22)

2’ω

ω ω

If A is symmetric and positive de¬nite, the SSOR method is convergent if

0 < ω < 2 (see [Hac94] for the proof). Typically, the SSOR method with an

optimal choice of the relaxation parameter converges more slowly than the

corresponding SOR method. However, the value of ρ(Bs (ω)) is less sensitive

to a choice of ω around the optimal value (in this respect, see the behavior

of the spectral radii of the two iteration matrices in Figure 4.1). For this

reason, the optimal value of ω that is chosen in the case of SSOR method

is usually the same used for the SOR method (for further details, we refer

to [You71]).

4.2 Linear Iterative Methods 135

1

0.9

0.8 SSOR

0.7

0.6

ρ 0.5

SOR

0.4

0.3

0.2

0.1

0

0 0.5 1 1.5 2

ω

FIGURE 4.1. Spectral radius of the iteration matrix of SOR and SSOR methods,

as a function of the relaxation parameter ω for the matrix tridiag10 (’1, 2, ’1)

4.2.7 Implementation Issues

We provide the programs implementing the Jacobi and Gauss-Seidel meth-

ods in their point form and with relaxation.

In Program 15 the JOR method is implemented (the Jacobi method is

obtained as a special case setting omega = 1). The stopping test monitors

the Euclidean norm of the residual at each iteration, normalized to the

value of the initial residual.

Notice that each component x(i) of the solution vector can be computed

independently; this method can thus be easily parallelized.

Program 15 - JOR : JOR method

function [x, iter]= jor ( a, b, x0, nmax, toll, omega)

[n,n]=size(a);

iter = 0; r = b - a * x0; r0 = norm(r); err = norm (r); x = x0;

while err > toll & iter < nmax

iter = iter + 1;

for i=1:n

s = 0;

for j = 1:i-1, s = s + a (i,j) * x (j); end

for j = i+1:n, s = s + a (i,j) * x (j); end

x (i) = omega * ( b(i) - s) / a(i,i) + (1 - omega) * x(i);

end

r = b - a * x; err = norm (r) / r0;

end

Program 16 implements the SOR method. Taking omega=1 yields the

Gauss-Seidel method.

136 4. Iterative Methods for Solving Linear Systems

Unlike the Jacobi method, this scheme is fully sequential. However, it can

be e¬ciently implemented without storing the solution of the previous step,

with a saving of memory storage.

Program 16 - SOR : SOR method

function [x, iter]= sor ( a, b, x0, nmax, toll, omega)

[n,n]=size(a);

iter = 0; r = b - a * x0; r0 = norm (r); err = norm (r); xold = x0;

while err > toll & iter < nmax

iter = iter + 1;

for i=1:n

s = 0;

for j = 1:i-1, s = s + a (i,j) * x (j); end

for j = i+1:n

s = s + a (i,j) * xold (j);

end

x (i) = omega * ( b(i) - s) / a(i,i) + (1 - omega) * xold (i);

end

x = x™; xold = x; r = b - a * x; err = norm (r) / r0;

end

4.3 Stationary and Nonstationary Iterative

Methods

Denote by

RP = I ’ P’1 A

the iteration matrix associated with (4.7). Proceeding as in the case of

relaxation methods, (4.7) can be generalized introducing a relaxation (or

acceleration) parameter ±. This leads to the following stationary Richard-

son method

x(k+1) = x(k) + ±P’1 r(k) , k ≥ 0. (4.23)

More generally, allowing ± to depend on the iteration index, the nonsta-

tionary Richardson method or semi-iterative method given by

x(k+1) = x(k) + ±k P’1 r(k) , k ≥ 0. (4.24)

The iteration matrix at the k-th step for these methods (depending on k)

is

R(±k ) = I ’ ±k P’1 A,

4.3 Stationary and Nonstationary Iterative Methods 137

with ±k = ± in the stationary case. If P=I, the methods will be called

nonpreconditioned. The Jacobi and Gauss-Seidel methods can be regarded

as stationary Richardson methods with ± = 1, P = D and P = D ’ E,

respectively.

We can rewrite (4.24) (and, thus, also (4.23)) in a form of greater inter-

est for computation. Letting z(k) = P’1 r(k) (the so-called preconditioned

residual), we get x(k+1) = x(k) + ±k z(k) and r(k+1) = b ’ Ax(k+1) =

r(k) ’±k Az(k) . To summarize, a nonstationary Richardson method requires

at each k + 1-th step the following operations:

solve the linear system Pz(k) = r(k) ;

compute the acceleration parameter ±k ;

(4.25)

update the solution x(k+1) = x(k) + ±k z(k) ;

update the residual r(k+1) = r(k) ’ ±k Az(k) .

4.3.1 Convergence Analysis of the Richardson Method

Let us ¬rst consider the stationary Richardson methods for which ±k = ±

for k ≥ 0. The following convergence result holds.

Theorem 4.8 For any nonsingular matrix P, the stationary Richardson

method (4.23) is convergent i¬

2Re»i

> 1 ∀i = 1, . . . , n, (4.26)

±|»i |2

where »i ∈ C are the eigenvalues of P’1 A.

Proof. Let us apply Theorem 4.1 to the iteration matrix R± = I ’ ±P’1 A. The

condition |1 ’ ±»i | < 1 for i = 1, . . . , n yields the inequality

(1 ’ ±Re»i )2 + ±2 (Im»i )2 < 1

3

from which (4.26) immediately follows.

Let us notice that, if the sign of the real parts of the eigenvalues of P’1 A

is not constant, the stationary Richardson method cannot converge.

More speci¬c results can be obtained provided that suitable assumptions

are made on the spectrum of P’1 A.

Theorem 4.9 Assume that P is a nonsingular matrix and that P’1 A has

positive real eigenvalues, ordered in such a way that »1 ≥ »2 ≥ . . . ≥

»n > 0. Then, the stationary Richardson method (4.23) is convergent i¬

0 < ± < 2/»1 . Moreover, letting

2

±opt = (4.27)

»1 + »n

138 4. Iterative Methods for Solving Linear Systems

the spectral radius of the iteration matrix R± is minimum if ± = ±opt , with

»1 ’ » n

ρopt = min [ρ(R± )] = . (4.28)

»1 + »n

±

Proof. The eigenvalues of R± are given by »i (R± ) = 1 ’ ±»i , so that (4.23) is

convergent i¬ |»i (R± )| < 1 for i = 1, . . . , n, that is, if 0 < ± < 2/»1 . It follows

(see Figure 4.2) that ρ(R± ) is minimum when 1 ’ ±»n = ±»1 ’ 1, that is, for

± = 2/(»1 + »n ), which furnishes the desired value for ±opt . By substitution, the

3

desired value of ρopt is obtained.

|1 ’ ±»1 |

ρ=1

|1 ’ ±»k |

ρopt

|1 ’ ±»n |

1 2 1

±opt ±

»1 »1 »n

FIGURE 4.2. Spectral radius of R± as a function of the eigenvalues of P’1 A

If P’1 A is symmetric positive de¬nite, it can be shown that the convergence

of the Richardson method is monotone with respect to either · 2 and · A .

In such a case, using (4.28), we can also relate ρopt to K2 (P’1 A) as follows

K2 (P’1 A) ’ 1 2 A’1 P 2

(4.29)

ρopt = , ±opt = .

K2 (P’1 A) + 1 K2 (P’1 A) + 1

The choice of a suitable preconditioner P is, therefore, of paramount im-

portance for improving the convergence of a Richardson method. Of course,

such a choice should also account for the need of keeping the computational

e¬ort as low as possible. In Section 4.3.2, some preconditioners of common

use in practice will be described.

Corollary 4.1 Let A be a symmetric positive de¬nite matrix. Then, the

non preconditioned stationary Richardson method is convergent and

¤ ρ(R± ) e(k) k ≥ 0.

e(k+1) A, (4.30)

A

4.3 Stationary and Nonstationary Iterative Methods 139

The same result holds for the preconditioned Richardson method, provided

that the matrices P, A and P’1 A are symmetric positive de¬nite.

Proof. The convergence is a consequence of Theorem 4.8. Moreover, we notice

that

¤ A1/2 R± A’1/2

e(k+1) = R± e(k) = A1/2 R± e(k) A1/2 e(k) 2.

A A 2 2

The matrix R± is symmetric positive de¬nite and is similar to A1/2 R± A’1/2 .

Therefore,

A1/2 R± A’1/2 2 = ρ(R± ).

The result (4.30) follows by noting that A1/2 e(k) 2 = e(k) A . A similar proof

can be carried out in the preconditioned case, provided we replace A with P’1 A.

3

Finally, the inequality (4.30) holds even if only P and A are symmetric

positive de¬nite (for the proof, see [QV94], Chapter 2).

4.3.2 Preconditioning Matrices

All the methods introduced in the previous sections can be cast in the form

(4.2), so that they can be regarded as being methods for solving the system

(I ’ B)x = f = P’1 b.

On the other hand, since B=P’1 N, system (3.2) can be equivalently refor-

mulated as

P’1 Ax = P’1 b. (4.31)

The latter is the preconditioned system, being P the preconditioning matrix

or left preconditioner. Right and centered preconditioners can be introduced

as well, if system (3.2) is transformed, respectively, as

AP’1 y = b, y = Px,

or

P’1 AP’1 y = P’1 b, y = PR x.

L R L

There are point preconditioners or block preconditioners, depending on

whether they are applied to the single entries of A or to the blocks of

a partition of A. The iterative methods considered so far correspond to

¬xed-point iterations on a left-preconditioned system. As stressed by (4.25),

computing the inverse of P is not mandatory; actually, the role of P is to

“preconditioning” the residual r(k) through the solution of the additional

system Pz(k) = r(k) .

140 4. Iterative Methods for Solving Linear Systems

Since the preconditioner acts on the spectral radius of the iteration ma-

trix, it would be useful to pick up, for a given linear system, an optimal

preconditioner, i.e., a preconditioner which is able to make the number of

iterations required for convergence independent of the size of the system.

Notice that the choice P=A is optimal but, trivially, “ine¬cient”; some

alternatives of greater computational interest will be examined below.

There is a lack of general theoretical results that allow to devise optimal

preconditioners. However, an established “rule of thumb” is that P is a

good preconditioner for A if P’1 A is near to being a normal matrix and if

its eigenvalues are clustered within a su¬ciently small region of the com-

plex ¬eld. The choice of a preconditioner must also be guided by practical

considerations, noticeably, its computational cost and its memory require-

ments.

Preconditioners can be divided into two main categories: algebraic and

functional preconditioners, the di¬erence being that the algebraic precon-

ditioners are independent of the problem that originated the system to

be solved, and are actually constructed via algebraic procedure, while the

functional preconditioners take advantage of the knowledge of the problem

and are constructed as a function of it. In addition to the preconditioners

already introduced in Section 4.2.6, we give a description of other algebraic

preconditioners of common use.

1. Diagonal preconditioners: choosing P as the diagonal of A is generally

e¬ective if A is symmetric positive de¬nite. A usual choice in the non

symmetric case is to set

« 1/2

n

pii = a2 .

ij

j=1

Block diagonal preconditioners can be constructed in a similar man-

ner. We remark that devising an optimal diagonal preconditioner is

far from being trivial, as previously noticed in Section 3.12.1 when

dealing with the scaling of a matrix.

2. Incomplete LU factorization (shortly ILU) and Incomplete Cholesky

factorization (shortly IC).

An incomplete factorization of A is a process that computes P =

Lin Uin , where Lin is a lower triangular matrix and Uin is an upper

triangular matrix. These matrices are approximations of the exact

matrices L, U of the LU factorization of A and are chosen in such a

way that the residual matrix R = A’Lin Uin satis¬es some prescribed

requirements, such as having zero entries in speci¬ed locations.