precision iterative re¬nement (FPR).

It can be shown that, if |A’1 | |L| |U| ∞ is su¬ciently small, then at

each step i of the algorithm, the relative error x’x(i) ∞ / x ∞ is reduced

by a factor ρ, which is given by

ρ 2 n cond(A, x)u (FPR),

ρ (MPR),

u

where ρ is independent of the condition number of A in the case of MPR.

Slow convergence of FPR is a clear indication of the ill-conditioning of the

112 3. Direct Methods for the Solution of Linear Systems

matrix, as it can be shown that, if p is the number of iterations for the

method to converge, then K∞ (A) β t(1’1/p) .

Even if performed in ¬xed precision, iterative re¬nement is worth using

since it improves the overall stability of any direct method for solving the

system. We refer to [Ric81], [Ske80], [JW77] [Ste73], [Wil63] and [CMSW79]

for an overview of this subject.

3.13 Undetermined Systems

We have seen that the solution of the linear system Ax=b exists and is

unique if n = m and A is nonsingular. In this section we give a meaning

to the solution of a linear system both in the overdetermined case, where

m > n, and in the underdetermined case, corresponding to m < n. We

notice that an underdetermined system generally has no solution unless

the right side b is an element of range(A).

For a detailed presentation, we refer to [LH74], [GL89] and [Bj¨88].

o

—

Given A∈ R with m ≥ n, b∈ R , we say that x ∈ R is a solution

m—n m n

of the linear system Ax=b in the least-squares sense if

¦(x— ) = Ax— ’ b ¤ min Ax ’ b

2 2

= min ¦(x). (3.73)

2 2

n n

x∈R x∈R

The problem thus consists of minimizing the Euclidean norm of the resid-

ual. The solution of (3.73) can be found by imposing the condition that the

gradient of the function ¦ in (3.73) must be equal to zero at x— . From

¦(x) = (Ax ’ b)T (Ax ’ b) = xT AT Ax ’ 2xT AT b + bT b,

we ¬nd that

∇¦(x— ) = 2AT Ax— ’ 2AT b = 0,

from which it follows that x— must be the solution of the square system

AT Ax— = AT b (3.74)

known as the system of normal equations. The system is nonsingular if

A has full rank and in such a case the least-squares solution exists and

is unique. We notice that B = AT A is a symmetric and positive de¬nite

matrix. Thus, in order to solve the normal equations, one could ¬rst com-

pute the Cholesky factorization B = HT H and then solve the two systems

HT y = AT b and Hx— = y. However, due to roundo¬ errors, the com-

putation of AT A may be a¬ected by a loss of signi¬cant digits, with a

consequent loss of positive de¬niteness or nonsingularity of the matrix, as

happens in the following example (implemented in MATLAB) where for a

3.13 Undetermined Systems 113

the corresponding matrix f l(AT A) turns out to

matrix A with full rank,

be singular

®

1 1

1 1

A = ° 2’27 » , f l(AT A) =

0 .

1 1

2’27

0

Therefore, in the case of ill-conditioned matrices it is more convenient to

utilize the QR factorization introduced in Section 3.4.3. Indeed, the follow-

ing result holds.

Theorem 3.8 Let A ∈ Rm—n , with m ≥ n, be a full rank matrix. Then

the unique solution of (3.73) is given by

x— = R’1 QT b

˜˜ (3.75)

˜ ˜

where R ∈ Rn—n and Q ∈ Rm—n are the matrices de¬ned in (3.48) starting

from the QR factorization of A. Moreover, the minimum of ¦ is given by

m

—

[(QT b)i ]2 .

¦(x ) =

i=n+1

Proof. The QR factorization of A exists and is unique since A has full rank.

Thus, there exist two matrices, Q∈ Rm—m and R∈ Rm—n such that A=QR, where

Q is orthogonal. Since orthogonal matrices preserve the Euclidean scalar product

(see Property 1.8), it follows that

Ax ’ b = Rx ’ QT b 2 .

2

2 2

Recalling that R is upper trapezoidal, we have

m

˜ ˜

Rx ’ Q b = Rx ’ QT b

T 2 2

[(QT b)i ]2 ,

+

2 2

i=n+1

so that the minimum is achieved when x = x— . 3

For more details about the analysis of the computational cost the algo-

rithm (which depends on the actual implementation of the QR factoriza-

tion), as well as for results about its stability, we refer the reader to the

texts quoted at the beginning of the section.

If A does not have full rank, the solution techniques above fail, since in

this case if x— is a solution to (3.73), the vector x— + z, with z ∈ ker(A), is

a solution too. We must therefore introduce a further constraint to enforce

the uniqueness of the solution. Typically, one requires that x— has minimal

Euclidean norm, so that the least-squares problem can be formulated as

¬nd x— ∈ Rn with minimal Euclidean norm such that

(3.76)

Ax— ’ b ¤ min Ax ’ b 2 .

2

2

2 n

x∈R

114 3. Direct Methods for the Solution of Linear Systems

This problem is consistent with (3.73) if A has full rank, since in this case

(3.73) has a unique solution which necessarily must have minimal Euclidean

norm.

The tool for solving (3.76) is the singular value decomposition (or SVD,

see Section 1.9), for which the following theorem holds.

Theorem 3.9 Let A ∈ Rm—n with SVD given by A = UΣVT . Then the

unique solution to (3.76) is

x— = A† b (3.77)

where A† is the pseudo-inverse of A introduced in De¬nition 1.15.

Proof. Using the SVD of A, problem (3.76) is equivalent to ¬nding w = VT x

such that w has minimal Euclidean norm and

Σw ’ UT b ¤ Σy ’ UT b 2 , ∀y ∈ Rn .

2

2 2

If r is the number of nonzero singular values σi of A, then

r m

2 2

Σw ’ U b σi wi ’ (U b)i

T 2 T

(UT b)i

= + ,

2

i=1 i=r+1

which is minimum if wi = (UT b)i /σi for i = 1, . . . , r. Moreover, it is clear that

among the vectors w of Rn having the ¬rst r components ¬xed, the one with

minimal Euclidean norm has the remaining n ’ r components equal to zero.

Thus the solution vector is w— = Σ† UT b, that is, x— = VΣ† UT b = A† b, where

Σ† is the diagonal matrix de¬ned in (1.11). 3

As for the stability of problem (3.76), we point out that if the matrix

A does not have full rank, the solution x— is not necessarily a continuous

function of the data, so that small changes on these latter might produce

large variations in x— . An example of this is shown below.

Example 3.11 Consider the system Ax = b with

® ®

10 1

A = ° 0 0 » , b = ° 2 » , rank(A) = 1.

00 3

Using the MATLAB function svd we can compute the SVD of A. Then computing

the pseudo-inverse, one ¬nds the solution vector x— = (1, 0)T . If we perturb the

null entry a22 , with the value 10’12 , the perturbed matrix has (full) rank 2

and the solution (which is unique in the sense of (3.73)) is now given by x— =

T

1, 2 · 1012 . •

We refer the reader to Section 5.8.3 for the approximate computation of

the SVD of a matrix.

3.14 Applications 115

In the case of underdetermined systems, for which m < n, if A has full

rank the QR factorization can still be used. In particular, when applied

to the transpose matrix AT , the method yields the solution of minimal

euclidean norm. If, instead, the matrix has not full rank, one must resort

to SVD.

Remark 3.8 If m = n (square system), both SVD and QR factorization

can be used to solve the linear system Ax=b, as alternatives to GEM.

Even though these algorithms require a number of ¬‚ops far superior to

GEM (SVD, for instance, requires 12n3 ¬‚ops), they turn out to be more

accurate when the system is ill-conditioned and nearly singular.

Example 3.12 Compute the solution to the linear system H15 x=b, where H15

is the Hilbert matrix of order 15 (see (3.32)) and the right side is chosen in

such a way that the exact solution is the unit vector x = 1. Using GEM with

partial pivoting yields a solution a¬ected by a relative error larger than 100%. A

solution of much better quality is obtained by passing through the computation

of the pseudo-inverse, where the entries in Σ that are less than 10’13 are set

•

equal to zero.

3.14 Applications

In this section we present two problems, suggested by structural mechan-

ics and grid generation in ¬nite element analysis, whose solutions require

solving large linear systems.

3.14.1 Nodal Analysis of a Structured Frame

Let us consider a structured frame which is made by rectilinear beams con-

nected among them through hinges (referred to as the nodes) and suitably

constrained to the ground. External loads are assumed to be applied at

the nodes of the frame and for any beam in the frame the internal actions

amount to a unique force of constant strength and directed as the beam

itself. If the normal stress acting on the beam is a traction we assume

that it has positive sign, otherwise the action has negative sign. Structured

frames are frequently employed as covering structures for large size public

buildings like exhibition stands, railway stations or airport halls.

To determine the internal actions in the frame, that are the unknowns

of the mathematical problem, a nodal analysis is used (see [Zie77]): the

equilibrium with respect to translation is imposed at every node of the

frame yielding a sparse and large-size linear system. The resulting matrix

has a sparsity pattern which depends on the numbering of the unknowns

and that can strongly a¬ect the computational e¬ort of the LU factorization

116 3. Direct Methods for the Solution of Linear Systems

due to ¬ll-in. We will show that the ¬ll-in can be dramatically reduced by

a suitable reordering of the unknowns.

The structure shown in Figure 3.9 is arc-shaped and is symmetric with

respect to the origin. The radii r and R of the inner and outer circles are

equal to 1 and 2, respectively. An external vertical load of unit size directed

downwards is applied at (0, 1) while the frame is constrained to ground

through a hinge at (’(r + R), 0) and a bogie at (r + R, 0). To generate

the structure we have partitioned the half unit circle in nθ uniform slices,

resulting in a total number of n = 2(nθ + 1) nodes and a matrix size of

m = 2n. The structure in Figure 3.9 has nθ = 7 and the unknowns are

numbered following a counterclockwise labeling of the beams starting from

the node at (1, 0).

We have represented the structure along with the internal actions com-

puted by solving the nodal equilibrium equations where the width of the

beams is proportional to the strength of the computed action. Black is

used to identify tractions whereas gray is associated with compressions. As

expected the maximum traction stress is attained at the node where the

external load is applied.

2.5

2

1.5

1

0.5

0

’0.5

’2 ’1.5 ’1 ’0.5 0 0.5 1 1.5 2

FIGURE 3.9. A structured frame loaded at the point (0, 1)

We show in Figure 3.10 the sparsity pattern of matrix A (left) and that

of the L-factor of its LU factorization with partial pivoting (right) in the

case nθ = 40 which corresponds to a size of 164 — 164. Notice the large

¬ll-in e¬ect arising in the lower part of L which results in an increase of

the nonzero entries from 645 (before the factorization) to 1946 (after the

factorization).

3.14 Applications 117

0 0

20 20

40 40

60 60

80 80

100 100

120 120

140 140

160 160

0 20 40 60 80 100 120 140 160 0 20 40 60 80 100 120 140 160

nz = 645 nz = 1946

FIGURE 3.10. Sparsity pattern of matrix A (left) and of the L-factor of the LU

factorization with partial pivoting (right) in the case nθ = 40

In view of the solution of the linear system by a direct method, the

increase of the nonzero entries demands for a suitable reordering of the

unknowns. For this purpose we use the MATLAB function symrcm which

implements the symmetric reverse Cuthill-McKee algorithm described in

Section 3.9.1. The sparsity pattern, after reordering, is shown in Figure 3.11

(left) while the L-factor of the LU factorization of the reordered matrix

is shown in Figure 3.11 (right). The results indicate that the reordering

procedure has “scattered” the sparsity pattern throughout the matrix with

a relatively modest increase of the nonzero entries from 645 to 1040.

0 0

20 20

40 40

60 60

80 80

100 100

120 120

140 140

160 160

0 20 40 60 80 100 120 140 160 0 20 40 60 80 100 120 140 160

nz = 645 nz = 1040

FIGURE 3.11. Sparsity pattern of matrix A (left) after a reordering with the sym-

metric reverse Cuthill-McKee algorithm and the L-factor of the LU factorization

of the reordered matrix with partial pivoting (right) in the case nθ = 40

The e¬ectiveness of the symmetric reverse Cuthill-McKee reordering pro-

cedure is demonstrated in Figure 3.12 which shows the number of nonzero

entries nz in the L-factor of A as a function of the size m of the matrix