<< . .

. 11
( : 95)



. . >>

for j = n:-1:2,
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,

<< . .

. 11
( : 95)



. . >>