<< . .

. 21
( : 95)

. . >>

e =B e +

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

(k+1) (k)
= bi ’
Aii xi Aij xj , i = 1, . . . , p,

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

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


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’E D’F .
PSSOR (ω) = (4.22)
ω ω

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



0.8 SSOR


ρ 0.5




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)
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);
r = b - a * x; err = norm (r) / r0;

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)
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);
x (i) = omega * ( b(i) - s) / a(i,i) + (1 - omega) * xold (i);
x = x™; xold = x; r = b - a * x; err = norm (r) / r0;

4.3 Stationary and Nonstationary Iterative
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)

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,
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 ;
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¬
> 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
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
±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
desired value of ρopt is obtained.

|1 ’ ±»1 |

|1 ’ ±»k |


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

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

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,


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

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-
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
pii =  a2  .

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.

<< . .

. 21
( : 95)

. . >>