the lower (upper) triangular part of M. Moreover, we assume that the

factorization process can be carried out without resorting to pivoting.

4.3 Stationary and Nonstationary Iterative Methods 141

The basic approach to incomplete factorization, consists of requiring

the approximate factors Lin and Uin to have the same sparsity pat-

tern as the L-part and U-part of A, respectively. A general algorithm

for constructing an incomplete factorization is to perform Gauss elim-

(k) (k)

ination as follows: at each step k, compute mik = aik /akk only if

aik = 0 for i = k + 1, . . . , n. Then, compute for j = k + 1, . . . , n

(k+1)

aij only if aij = 0. This algorithm is implemented in Program 17

where the matrices Lin and Uin are progressively overwritten onto

the L-part and U-part of A.

Program 17 - basicILU : Incomplete LU factorization

function [a] = basicILU(a)

[n,n]=size(a);

for k=1:n-1, for i=k+1:n,

if a(i,k) ˜= 0

a(i,k) = a(i,k) / a(k,k);

for j=k+1:n

if a(i,j) ˜= 0

a(i,j) = a(i,j) -a(i,k)*a(k,j);

end

end

end

end, end

We notice that having Lin and Uin with the same patterns as the

L and U-parts of A, respectively, does not necessarily imply that R

has the same sparsity pattern as A, but guarantees that rij = 0 if

aij = 0, as is shown in Figure 4.3.

The resulting incomplete factorization is known as ILU(0), where “0”

means that no ¬ll-in has been introduced in the factorization process.

An alternative strategy might be to ¬x the structure of Lin and Uin

irrespectively of that of A, in such a way that some computational

criteria are satis¬ed (for example, that the incomplete factors have

the simplest possible structure).

The accuracy of the ILU(0) factorization can obviously be improved

by allowing some ¬ll-in to arise, and thus, by accepting nonzero entries

in the factorization whereas A has elements equal to zero. To this

purpose, it is convenient to introduce a function, which we call ¬ll-

in level, that is associated with each entry of A and that is being

modi¬ed during the factorization process. If the ¬ll-in level of an

142 4. Iterative Methods for Solving Linear Systems

0

1

2

3

4

5

6

7

8

9

10

11

0 1 2 3 4 5 6 7 8 9 10 11

FIGURE 4.3. The sparsity pattern of the original matrix A is represented by the

squares, while the pattern of R = A’Lin Uin , computed by Program 17, is drawn

by the bullets

element is greater than an admissible value p ∈ N, the corresponding

entry in Uin or Lin is set equal to zero.

Let us explain how this procedure works, assuming that the matri-

ces Lin and Uin are progressively overwritten to A (as happens in

(k)

Program 4). The ¬ll-in level of an entry aij is denoted by levij ,

where the dependence on k is understood, and it should provide a

reasonable estimate of the size of the entry during the factorization

process. Actually, we are assuming that if levij = q then |aij | δq

(k)

with δ ∈ (0, 1), so that q is greater when |aij | is smaller.

At the starting step of the procedure, the level of the nonzero entries

of A and of the diagonal entries is set equal to 0, while the level of

the null entries is set equal to in¬nity. For any row i = 2, . . . , n, the

following operations are performed: if levik ¤ p, k = 1, . . . , i ’ 1, the

(k+1)

entry mik of Lin and the entries aij of Uin , j = i + 1, . . . , n, are

(k+1)

updated. Moreover, if aij = 0 the value levij is updated as being

the minimum between the available value of levij and levik +levkj +1.

(k+1) (k) (k)

The reason of this choice is that |aij | = |aij ’ mik akj | |δ levij ’

(k+1)

δ levik +levkj +1 |, so that one can assume that the size of |aij | is the

maximum between δ levij and δ levik +levkj +1 .

The above factorization process is called ILU(p) and turns out to be

extremely e¬cient (with p small) provided that it is coupled with a

suitable matrix reordering (see Section 3.9).

Program 18 implements the ILU(p) factorization; it returns in out-

put the approximate matrices Lin and Uin (overwritten to the input

matrix a), with the diagonal entries of Lin equal to 1, and the ma-

4.3 Stationary and Nonstationary Iterative Methods 143

trix lev containing the ¬ll-in level of each entry at the end of the

factorization.

Program 18 - ilup : ILU(p) factorization

function [a,lev] = ilup (a,p)

[n,n]=size(a);

for i=1:n, for j=1:n

if (a(i,j) ˜= 0) | (i==j)

lev(i,j)=0;

else

lev(i,j)=Inf;

end

end, end

for i=2:n,

for k=1:i-1

if lev(i,k) <= p

a(i,k)=a(i,k)/a(k,k);

for j=k+1:n

a(i,j)=a(i,j)-a(i,k)*a(k,j);

if a(i,j) ˜= 0

lev(i,j)=min(lev(i,j),lev(i,k)+lev(k,j)+1);

end

end

end

end

for j=1:n, if lev(i,j) > p, a(i,j) = 0; end, end

end

Example 4.4 Consider the matrix A ∈ R46—46 associated with the ¬nite

‚2· ‚2·

di¬erence approximation of the Laplace operator ∆· = ‚x2 + ‚y2 (see

Section 12.6). This matrix can be generated with the following MATLAB

commands: G=numgrid(™B™,10); A=delsq(G) and corresponds to the dis-

cretization of the di¬erential operator on a domain having the shape of the

exterior of a butter¬‚y and included in the square [’1, 1]2 (see Section 12.6).

The number of nonzero entries of A is 174. Figure 4.4 shows the pattern of

matrix A (drawn by the bullets) and the entries in the pattern added by

the ILU(1) and ILU(2) factorizations due to ¬ll-in (denoted by the squares

and the triangles, respectively). Notice that these entries are all contained

•

within the envelope of A since no pivoting has been performed.

The ILU(p) process can be carried out without knowing the actual

values of the entries of A, but only working on their ¬ll-in levels.

Therefore, we can distinguish between a symbolic factorization (the

generation of the levels) and an actual factorization (the computation

of the entries of ILU(p) starting from the informations contained in

144 4. Iterative Methods for Solving Linear Systems

0

5

10

15

20

25

30

35

40

45

0 5 10 15 20 25 30 35 40 45

FIGURE 4.4. Pattern of the matrix A in Example 4.4 (bullets); entries added by

the ILU(1) and ILU(2) factorizations (squares and triangles, respectively)

the level function). The scheme is thus particularly e¬ective when

several linear systems must be solved, with matrices having the same

structure but di¬erent entries.

On the other hand, for certain classes of matrices, the ¬ll-in level

does not always provide a sound indication of the actual size attained

by the entries. In such cases, it is better to monitor the size of the

entries of R by neglecting each time the entries that are too small.

(k+1)

For instance, one can drop out the entries aij such that

(k+1) (k+1) (k+1) 1/2

|aij | ¤ c|aii ajj | , i, j = 1, . . . , n,

with 0 < c < 1 (see [Axe94]).

In the strategies considered so far, the entries of the matrix that are

dropped out can no longer be recovered in the incomplete factoriza-

tion process. Some remedies exist for this drawback: for instance, at

the end of each k-th step of the factorization, one can sum, row by

row, the discarded entries to the diagonal entries of Uin . By doing

so, an incomplete factorization known as MILU (Modi¬ed ILU) is

obtained, which enjoys the property of being exact with respect to

the constant vectors, i.e., such that R1T = 0T (see [Axe94] for other

formulations). In the practice, this simple trick provides, for a wide

class of matrices, a better preconditioner than obtained with the ILU

method. In the case of symmetric positive de¬nite matrices one can

resort to the Modi¬ed Incomplete Cholesky Factorization (MICh).

We conclude by mentioning the ILUT factorization, which collects the

features of ILU(p) and MILU. This factorization can also include par-

tial pivoting by columns with a slight increase of the computational

4.3 Stationary and Nonstationary Iterative Methods 145

cost. For an e¬cient implementation of incomplete factorizations, we

refer to the MATLAB function luinc in the toolbox sparfun.

The existence of the ILU factorization is not guaranteed for all non-

singular matrices (see for an example [Elm86]) and the process stops

if zero pivotal entries arise. Existence theorems can be proved if A is

an M-matrix [MdV77] or diagonally dominant [Man80]. It is worth

noting that sometimes the ILU factorization turns out to be more

stable than the complete LU factorization [GM83].

3. Polynomial preconditioners: the preconditioning matrix is de¬ned as

P’1 = p(A),

where p is a polynomial in A, usually of low degree.

A remarkable example is given by Neumann polynomial precondi-

tioners. Letting A = D ’ C, we have A = (I ’ CD’1 )D, from which

A’1 = D’1 (I ’ CD’1 )’1 = D’1 (I + CD’1 + (CD’1 )2 + . . . ).

A preconditioner can then be obtained by truncating the series above

at a certain power p. This method is actually e¬ective only if ρ(CD’1 )

< 1, which is the necessary condition in order the series to be con-

vergent.

4. Least-squares preconditioners: A’1 is approximated by a least-squares

polynomial ps (A) (see Section 3.13). Since the aim is to make ma-

trix I ’ P’1 A as close as possible to the null matrix, the least-

squares approximant ps (A) is chosen in such a way that the function

•(x) = 1’ps (x)x is minimized. This preconditioning technique works

e¬ectively only if A is symmetric and positive de¬nite.

For further results on preconditioners, see [dV89] and [Axe94].

Example 4.5 Consider the matrix A∈ R324—324 associated with the ¬nite di¬er-

ence approximation of the Laplace operator on the square [’1, 1]2 . This matrix

can be generated with the following MATLAB commands: G=numgrid(™N™,20);

A=delsq(G). The condition number of the matrix is K2 (A) = 211.3. In Table

4.1 we show the values of K2 (P’1 A) computed using the ILU(p) and Neumann

preconditioners, with p = 0, 1, 2, 3. In the last case D is the diagonal part of A. •

Remark 4.2 Let A and P be real symmetric matrices of order n, with P

positive de¬nite. The eigenvalues of the preconditioned matrix P’1 A are

solutions of the algebraic equation

Ax = »Px, (4.32)

146 4. Iterative Methods for Solving Linear Systems

p ILU(p) Neumann

0 22.3 211.3

1 12 36.91

2 8.6 48.55

3 5.6 18.7

TABLE 4.1. Spectral condition numbers of the preconditioned matrix A of Ex-

ample 4.5 as a function of p

where x is an eigenvector associated with the eigenvalue ». Equation (4.32)

is an example of generalized eigenvalue problem (see Section 5.9 for a thor-

ough discussion) and the eigenvalue » can be computed through the fol-

lowing generalized Rayleigh quotient

(Ax, x)

»= .

(Px, x)

Applying the Courant-Fisher Theorem (see Section 5.11) yields

»min (A) »max (A)

¤»¤ . (4.33)

»max (P) »min (P)

Relation (4.33) provides a lower and upper bound for the eigenvalues of the

preconditioned matrix as a function of the extremal eigenvalues of A and

P, and therefore it can be pro¬tably used to estimate the condition number

of P’1 A.

4.3.3 The Gradient Method

The expression of the optimal parameter that has been provided in Theo-

rem 4.9 is of limited usefulness in practical computations, since it requires

the knowledge of the extremal eigenvalues of the matrix P’1 A. In the spe-

cial case of symmetric and positive de¬nite matrices, however, the optimal

acceleration parameter can be dynamically computed at each step k as

follows.

We ¬rst notice that, for such matrices, solving system (3.2) is equivalent

to ¬nding the minimizer x ∈ Rn of the quadratic form

1

¦(y) = yT Ay ’ yT b,

2

which is called the energy of system (3.2). Indeed, the gradient of ¦ is given

by

1

∇¦(y) = (AT + A)y ’ b = Ay ’ b. (4.34)

2

As a consequence, if ∇¦(x) = 0 then x is a solution of the original system.

Conversely, if x is a solution, then

1

¦(y) = ¦(x + (y ’ x)) = ¦(x) + (y ’ x)T A(y ’ x), ∀y ∈ Rn

2

4.3 Stationary and Nonstationary Iterative Methods 147

and thus, ¦(y) > ¦(x) if y = x, i.e. x is a minimizer of the functional ¦.

Notice that the previous relation is equivalent to

1

y’x = ¦(y) ’ ¦(x)

2

(4.35)

A

2

where · A is the A-norm or energy norm, de¬ned in (1.28).

The problem is thus to determine the minimizer x of ¦ starting from a

point x(0) ∈ Rn and, consequently, to select suitable directions along which

moving to get as close as possible to the solution x. The optimal direction,

that joins the starting point x(0) to the solution point x, is obviously un-

known a priori. Therefore, we must take a step from x(0) along another

direction d(0) , and then ¬x along this latter a new point x(1) from which

to iterate the process until convergence.

Thus, at the generic step k, x(k+1) is computed as

x(k+1) = x(k) + ±k d(k) , (4.36)

where ±k is the value which ¬xes the length of the step along d(k) . The most

natural idea is to take the descent direction of maximum slope ∇¦(x(k) ),

which yields the gradient method or steepest descent method.

On the other hand, due to (4.34), ∇¦(x(k) ) = Ax(k) ’ b = ’r(k) , so that

the direction of the gradient of ¦ coincides with that of residual and can

be immediately computed using the current iterate. This shows that the

gradient method, as well as the Richardson method, moves at each step k

along the direction d(k) = r(k) .

To compute the parameter ±k let us write explicitly ¦(x(k+1) ) as a func-

tion of a parameter ±

1 (k)

(x + ±r(k) )T A(x(k) + ±r(k) ) ’ (x(k) + ±r(k) )T b.

¦(x(k+1) ) =

2

Di¬erentiating with respect to ± and setting it equal to zero, yields the

desired value of ±k

T

r(k) r(k)

±k = (4.37)

T

r(k) Ar(k)

which depends only on the residual at the k-th step. For this reason, the

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)