b(j)=b(j)/U(j,j); b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);

end; b(1) = b(1)/U(1,1);

When large triangular systems must be solved, only the triangular portion

of the matrix should be stored leading to considerable saving of memory

resources.

3.2 Solution of Triangular Systems 67

3.2.2 Rounding Error Analysis

The analysis developed so far has not accounted for the presence of round-

ing errors. When including these, the forward and backward substitution

algorithms no longer yield the exact solutions to the systems Lx=b and

Uy=b, but rather provide approximate solutions x that can be regarded

as being exact solutions to the perturbed systems

(L + δL)x = b, (U + δU)x = b,

where δL = (δlij ) and δU = (δuij ) are perturbation matrices. In order

to apply the estimates (3.9) carried out in Section 3.1.2, we must provide

estimates of the perturbation matrices, δL and δU, as a function of the

entries of L and U, of their size and of the characteristics of the ¬‚oating-

point arithmetic. For this purpose, it can be shown that

nu

|δT| ¤ |T|, (3.24)

1 ’ nu

where T is equal to L or U, u = 1 β 1’t is the roundo¬ unit de¬ned in (2.34).

2

Clearly, if nu < 1 from (3.24) it turns out that, using a Taylor expansion,

|δT| ¤ nu|T| + O(u2 ). Moreover, from (3.24) and (3.9) it follows that, if

nuK(T) < 1, then

x’x nuK(T)

¤ = nuK(T) + O(u2 ) (3.25)

1 ’ nuK(T)

x

for the norms · 1 , · ∞ and the Frobenius norm. If u is su¬ciently

small (as typically happens), the perturbations introduced by the rounding

errors in the solution of a triangular system can thus be neglected. As

a consequence, the accuracy of the solution computed by the forward or

backward substitution algorithm is generally very high.

These results can be improved by introducing some additional assump-

tions on the entries of L or U. In particular, if the entries of U are such

that |uii | ≥ |uij | for any j > i, then

nu

|xi ’ xi | ¤ 2n’i+1 max|xj |, 1 ¤ i ¤ n.

1 ’ nu j≥i

The same result holds if T=L, provided that |lii | ≥ |lij | for any j < i, or if

L and U are diagonally dominant. The previous estimates will be employed

in Sections 3.3.1 and 3.4.2.

For the proofs of the results reported so far, see [FM67], [Hig89] and

[Hig88].

3.2.3 Inverse of a Triangular Matrix

The algorithm (3.23) can be employed to explicitly compute the inverse

of an upper triangular matrix. Indeed, given an upper triangular matrix

68 3. Direct Methods for the Solution of Linear Systems

U, the column vectors vi of the inverse V=(v1 , . . . , vn ) of U satisfy the

following linear systems

Uvi = ei , i = 1, . . . , n (3.26)

where {ei } is the canonical basis of Rn (de¬ned in Example 1.3). Solving

for vi thus requires the application of algorithm (3.23) n times to (3.26).

This procedure is quite ine¬cient since at least half the entries of the

inverse of U are null. Let us take advantage of this as follows. Denote by

vk = (v1k , . . . , vkk )T the vector of size k such that

U(k) vk = lk (3.27)

k = 1, . . . , n

where U(k) is the principal submatrix of U of order k and lk the vector of

Rk having null entries, except the ¬rst one which is equal to 1. Systems

(3.27) are upper triangular, but have order k and can be again solved using

the method (3.23). We end up with the following inversion algorithm for

upper triangular matrices: for k = n, n ’ 1, . . . , 1 compute

vkk = u’1 ,

kk

k

(3.28)

’u’1 uij vjk , for i = k ’ 1, k ’ 2, . . . , 1.

vik = ii

j=i+1

At the end of this procedure the vectors vk furnish the non vanishing entries

of the columns of U’1 . The algorithm requires about n3 /3 + (3/4)n2 ¬‚ops.

Once again, due to rounding errors, the algorithm (3.28) no longer yields

the exact solution, but an approximation of it. The error that is introduced

can be estimated using the backward a priori analysis carried out in Section

3.1.3.

A similar procedure can be constructed from (3.22) to compute the in-

verse of a lower triangular system.

3.3 The Gaussian Elimination Method (GEM) and

LU Factorization

The Gaussian elimination method aims at reducing the system Ax=b to an

equivalent system (that is, having the same solution) of the form Ux=b,

where U is an upper triangular matrix and b is an updated right side

vector. This latter system can then be solved by the backward substitution

method. Let us denote the original system by A(1) x = b(1) . During the

reduction procedure we basically employ the property which states that

replacing one of the equations by the di¬erence between this equation and

another one multiplied by a non null constant yields an equivalent system

(i.e., one with the same solution).

3.3 The Gaussian Elimination Method (GEM) and LU Factorization 69

Thus, consider a nonsingular matrix A ∈ Rn—n , and suppose that the

diagonal entry a11 is non vanishing. Introducing the multipliers

(1)

ai1

mi1 = , i = 2, 3, . . . , n,

(1)

a11

(1)

where aij denote the elements of A(1) , it is possible to eliminate the un-

known x1 from the rows other than the ¬rst one by simply subtracting

from row i, with i = 2, . . . , n, the ¬rst row multiplied by mi1 and doing

the same on the right side. If we now de¬ne

(2) (1) (1)

aij = aij ’ mi1 a1j , i, j = 2, . . . , n,

(2) (1) (1)

’ mi1 b1 ,

bi = bi i = 2, . . . , n,

(1)

denote the components of b(1) , we get a new system of the form

where bi

® ® ®

(1) (1) (1) (1)

a11 a12 ... a1n b1

x1

(2) (2) (2)

x2

0 a22 ... a2n b2

= ,

°

.

. . . .

»

.

. . . .

° » ° »

.

. . . .

(2) (2) (2)

xn

0 an2 ... ann bn

which we denote by A(2) x = b(2) , that is equivalent to the starting one.

Similarly, we can transform the system in such a way that the unknown

x2 is eliminated from rows 3, . . . , n. In general, we end up with the ¬nite

sequence of systems

A(k) x = b(k) , 1 ¤ k ¤ n, (3.29)

where, for k ≥ 2, matrix A(k) takes the following form

®

(1) (1) (1)

a11 a12 ... ... ... a1n

(2) (2)

0 a22 a2n

. .

..

. .

.

. .

(k)

= ,

A

(k) (k)

0 ... 0 akk ... akn

. . . .

° »

. . . .

. . . .

(k) (k)

0 ... 0 ank ... ann

70 3. Direct Methods for the Solution of Linear Systems

(i)

having assumed that aii = 0 for i = 1, . . . , k ’ 1. It is clear that for k = n

we obtain the upper triangular system A(n) x = b(n)

® (1) ® ® (1)

(1) (1)

a11 a12 . . . . . . a1n b1

x1

(2) (2)

a2n x2

(2)

0 b2

a22

.

. .

. . . .

.. .

. . = . .

. .

. .

.

.. . »°. .» °.»

°0 . . .

(n) (n)

xn

0 ann bn

Consistently with the notations that have been previously introduced, we

(k)

denote by U the upper triangular matrix A(n) . The entries akk are called

pivots and must obviously be non null for k = 1, . . . , n ’ 1.

In order to highlight the formulae which transform the k-th system into

(k)

the k + 1-th one, for k = 1, . . . , n ’ 1 we assume that akk = 0 and de¬ne

the multiplier

(k)

aik

mik = , i = k + 1, . . . , n. (3.30)

(k)

akk

Then we let

(k+1) (k) (k)

= aij ’ mik akj , i, j = k + 1, . . . , n

aij

(3.31)

(k+1) (k) (k)

’

bi = bi mik bk , i = k + 1, . . . , n.

Example 3.2 Let us use GEM to solve the following system

± 1 1 11

x 1 + 2 x2 + x =

33

6

(1) (1) 1

x1 + 1 x2 1 13

+ x =

(A x = b ) ,

43

2 3 12

1

+ 1 x2 1 47

x + x =

31 53

4 60

which admits the solution x=(1, 1, 1)T . At the ¬rst step we compute the mul-

tipliers m21 = 1/2 and m31 = 1/3, and subtract from the second and third

equation of the system the ¬rst row multiplied by m21 and m31 , respectively. We

obtain the equivalent system

± 1 1 11

x1 + x + x =

22 33

6

(A(2) x = b(2) ) 1 1 1

0 + 12 x2 + 12 x3 = .

6

1 4 31

0 + 12 x2 + 45 x3 = 180

If we now subtract the second row multiplied by m32 = 1 from the third one, we

end up with the upper triangular system

± 1 1 11

x1 + x + x =

22 33

6

(A(3) x = b(3) ) 1 1 1

0 + 12 x2 + x = ,

12 3

6

1 1

0+ 0 + 180 x3 = 180

3.3 The Gaussian Elimination Method (GEM) and LU Factorization 71

from which we immediately compute x3 = 1 and then, by back substitution, the

•

remaining unknowns x1 = x2 = 1.

Remark 3.2 The matrix in Example 3.2 is called the Hilbert matrix of

order 3. In the general n — n case, its entries are

hij = 1/(i + j ’ 1), i, j = 1, . . . , n. (3.32)

As we shall see later on, this matrix provides the paradigm of an ill-

conditioned matrix.

To complete Gaussian elimination 2(n ’ 1)n(n + 1)/3 + n(n ’ 1) ¬‚ops are

required, plus n2 ¬‚ops to backsolve the triangular system U x = b(n) .

Therefore, about (2n3 /3 + 2n2 ) ¬‚ops are needed to solve the linear sys-

tem using GEM. Neglecting the lower order terms, we can state that the

Gaussian elimination process has a cost of 2n3 /3 ¬‚ops.

(k)

As previously noticed, GEM terminates safely i¬ the pivotal elements akk ,

for k = 1, . . . , n ’ 1, are non vanishing. Unfortunately, having non null

diagonal entries in A is not enough to prevent zero pivots to arise during

the elimination process. For example, matrix A in (3.33) is nonsingular and

has nonzero diagonal entries

® ®

12 3

123

A = ° 2 4 5 » , A(2) = ° 0 0 ’1 » . (3.33)

0 ’6 ’12

789

Nevertheless, when GEM is applied, it is interrupted at the second step

(2)

since a22 = 0.

More restrictive conditions on A are thus needed to ensure the appli-

cability of the method. We shall see in Section 3.3.1 that if the leading

dominating minors di of A are nonzero for i = 1, . . . , n ’ 1, then the corre-

(i)

sponding pivotal entries aii must necessarily be non vanishing. We recall

that di is the determinant of Ai , the i-th principal submatrix made by the

¬rst i rows and columns of A. The matrix in the previous example does

not satisfy this condition, having d1 = 1 and d2 = 0.

Classes of matrices exist such that GEM can be always safely employed in

its basic form (3.31). Among them, we recall the following ones:

1. matrices diagonally dominant by rows;

2. matrices diagonally dominant by columns. In such a case one can even

show that the multipliers are in module less than or equal to 1 (see

Property 3.2);

3. matrices symmetric and positive de¬nite (see Theorem 3.6).

For a rigorous derivation of these results, we refer to the forthcoming sec-

tions.

72 3. Direct Methods for the Solution of Linear Systems

3.3.1 GEM as a Factorization Method

In this section we show how GEM is equivalent to performing a factorization

of the matrix A into the product of two matrices, A=LU, with U=A(n) .

Since L and U depend only on A and not on the right hand side, the same

factorization can be reused when solving several linear systems having the

same matrix A but di¬erent right hand side b, with a considerable reduction

of the operation count (indeed, the main computational e¬ort, about 2n3 /3

¬‚ops, is spent in the elimination procedure).

Let us go back to Example 3.2 concerning the Hilbert matrix H3 . In

practice, to pass from A(1) =H3 to the matrix A(2) at the second step, we

have multiplied the system by the matrix

® ®

100 10 0

M1 = ’ 1 1 0 = ’m21 1 0 .

° »° »

2

’m31 0

’3 0 1

1

1

Indeed,

®

1 1

1 2 3

(1)

= A(2) .

= 0 1 1

M1 A = M1 A

° »

12 12

1 4

0 12 45

Similarly, to perform the second (and last) step of GEM, we must multiply

A(2) by the matrix

® ®

1 00 1 00

M2 = 0 1 0 = 0 1 0 ,

° »° »

0 ’1 1 0 ’m32 1

where A(3) = M2 A(2) . Therefore

M2 M1 A = A(3) = U. (3.34)

On the other hand, matrices M1 and M2 are lower triangular, their product

is still lower triangular, as is their inverse; thus, from (3.34) one gets

A = (M2 M1 )’1 U = LU,