<< . .

. 18
( : 95)

. . >>

mixed-precision iterative re¬nement (shortly, MPR), as compared to ¬xed-
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),

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].

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
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

[(QT b)i ]2 .
¦(x ) =

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

Recalling that R is upper trapezoidal, we have
˜ ˜
Rx ’ Q b = Rx ’ QT b
T 2 2
[(QT b)i ]2 ,
2 2

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
Ax— ’ b ¤ min Ax ’ b 2 .
2 n
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
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

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
= + ,
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— =
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 ’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
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

<< . .

. 18
( : 95)

. . >>