<< . .

. 22
( : 95)

. . >>

For a given matrix M, the L-part (U-part) of M will mean henceforth
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
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)
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

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

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
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
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
entry mik of Lin and the entries aij of Uin , j = i + 1, . . . , n, are
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 ’
δ 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

Program 18 - ilup : ILU(p) factorization

function [a,lev] = ilup (a,p)
for i=1:n, for j=1:n
if (a(i,j) ˜= 0) | (i==j)
end, end
for i=2:n,
for k=1:i-1
if lev(i,k) <= p
for j=k+1:n
if a(i,j) ˜= 0
for j=1:n, if lev(i,j) > p, a(i,j) = 0; 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

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.
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-
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
We ¬rst notice that, for such matrices, solving system (3.2) is equivalent
to ¬nding the minimizer x ∈ Rn of the quadratic form
¦(y) = yT Ay ’ yT b,
which is called the energy of system (3.2). Indeed, the gradient of ¦ is given
∇¦(y) = (AT + A)y ’ b = Ay ’ b. (4.34)
As a consequence, if ∇¦(x) = 0 then x is a solution of the original system.
Conversely, if x is a solution, then
¦(y) = ¦(x + (y ’ x)) = ¦(x) + (y ’ x)T A(y ’ x), ∀y ∈ Rn
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

y’x = ¦(y) ’ ¦(x)
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) ) =
Di¬erentiating with respect to ± and setting it equal to zero, yields the
desired value of ±k
r(k) r(k)
±k = (4.37)
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)
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)

<< . .

. 22
( : 95)

. . >>