Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 35
Текст из файла (страница 35)
Show that a boundholds of the form in Theorem 8.9, namely, for(8.20)Give an example of a triangular system solver for which (8.20) is not satisfied.PreviousHomeChapter 9LU Factorization and Linear EquationsIt appears that Gauss and Doolittle appliedthe method only to symmetric equations.More recent authors, for example, Aitken, Banachiewicz, Dwyer, and Crout . . .have emphasized the use of the method, or variations of it,in connection with non-symmetric problems . .
.Banachiewicz . . . saw the point . . .that the basic problem is really one of matrix factorization,or “decomposition” as he called it.-PAUL S. DWYER, Linear Computations (1951)Intolerable pivot-growth [with partial pivoting] is a phenomenon that happensonly to numerical analysts who are looking for that phenomenon.-WILLIAM M.
KAHAN, Numerical Linear Algebra (1966)By 1949 the major components of thePilot ACE were complete and undergoing trials . . .During 1951 a programme for solving simultaneouslinear algebraic equations was used for the first time.26th June, 1951 was a landmark in the history of the machine,for on that day it first rivalled alternative computing methodsby yielding by 3 p.m. the solution toa set of 17 equations submitted the same morning.-MICHAEL WOODGER, The History and Present Use ofDigital Computers at the National Physical Laboratory (1958).The closer one looks,the more subtle and remarkable Gaussian elimination appears.-LLOYD N.
TREFETHEN, Three Mysteries of Gaussian Elimination (1985)169Next170LU FACTORIZATIONANDLINEAR E QUATIONS9.1. Gaussian EliminationWe begin by giving a traditional description of Gaussian elimination (GE) forsolving a linear system Ax = b, whereis nonsingular.The strategy of GE is to reduce a problem we can’t solve (a full linearsystem) to one that we can (a triangular system), using elementary row operations.
There are n - 1 stages, beginning with A(1) := A, b (1) := b, andfinishing with the upper triangular system A(n)x = b (n) .At the kth stage we have converted the original system to A(k)x = b(k),withupper triangular. The purpose of the kth stage ofthe elimination is to zero the elements below the diagonal in the kth columnof A ( k ). This is accomplished by the operationsAt the end of the (n - 1)stwhere the multipliersstage we have the upper triangular system A(n)x = b (n) , which is solved byback substitution. For an n × n matrix, GE requires 2n 3/3 flops.There are two problems with the method as described. First. there isa breakdown with division by zero ifSecond, if we are working infinite precision and some multiplier mik is large, then there is a possible loss ofsignificance: in the subtractionlow-order digits ofcould belost.
Losing these digits could correspond to making a relatively large changeto the original matrix A. The simplest example of this phenomenon is for thematrixhere.andifwhich wouldbe the exact answer if we changed a22 from 1 to 0.These observations motivate the strategy of partial pivoting.
At the startof the kth stage, the kth and rth rows are interchanged, wherePartial pivoting ensures that the multipliers are nicely bounded:A more expensive pivoting strategy, which interchanges both rows andcolumns, is complete pivoting.9.1 GAUSSIAN ELIMINATION171At the start of the kth stage, rows k and r and columns k and s areinterchanged, whereNote that this requires O(n3 ) comparisons in total, compared with O(n2 )for partial pivoting. Because of the searching overhead, and because partialpivoting works so well, complete pivoting is rarely used in practice.Much insight into GE is obtained by expressing it in matrix notation. Wecan writeThe matrix Mk can be expressed compactly as Mk =where ek isfor i < k.
To invert Mk, just flip the signsthe kth unit vector andof the multipliers:Overall,and soThe conclusion is that GE computes an LU factorization of A : A = LU,where L is unit lower triangular and U is upper triangular.We introduce the shorthand notation Ak := A(1:k, 1:k).Theorem 9.1. There exists a unique LU factorization ofif andonly if Ak is nonsingular for k = 1:n - 1. If Ak is singular for some 1 < k <n - 1 then the factorization may exist, but if so it is not unique.172LU F ACTORIZATIONANDLINEAR E QUATIONSProof. Suppose Ak is nonsingular for k = 1:n - 1.
The existence of anLU factorization can be proved by. examining the steps of GE. but a moreelegant proof, which also gives uniqueness, can be obtained by an inductive bordering construct ion. Suppose Ak-1 has the unique LU factorizationA k-1 = Lk- 1 U k-1 (this supposition clearly holds for k - 1 = 1). We look fora factorizationThe equations to be satisfied are Lk- 1u =and a kk = 1 Tu +u k k . The matrices Lk-1 and Uk-1 are nonsingular, since 0det(A k-1) =det(L k-1) det(U k-1). so the equations have a unique solution.
completing theinduction.We prove the converse, under the assumption that A is nonsingular; for thecase A singular see Problem 9.1. Suppose an LU factorization exists. ThenAk = LkUk for k = 1:n, which givesdet(A k) = det(U k) = u 11 . . . ukk.(9.1)Setting k = n we find that 0det(A) = u 11 . . . unn, and hence det(A k) =u 11 . . . ukk 0, k = 1:n - 1.Examples that illustrate the last part of the theorem arewhich does not have an LU factorizawhich holds for any l, andtion.Visually, the condition of Theorem 9.1 is (for n = 5) that the indicatedsubmatrices must be nonsingular:From (9.1) follows the expression u kk = det(A k)/det(A k-l). In fact, allthe elements of L and U can be expressed by determinant al formulae (see,e.g., Gantmacher [413, 1959, p. 35] or Householder [587, 1964.
p. 11]):(9.2a)(9.2b)9.1 GAUSSIAN ELIMINATION173The effect of partial pivoting is easily seen by considering the case n = 4.We havewhere, for example,For k =1,2,3,is the same as Mk except the multipliers are interchanged. Hence,for n = 4, GE with partial pivoting (GEPP) applied to A is equivalent to GEwithout pivoting applied to the row-permuted matrix PA. This conclusionis true for any n: GEPP computes a factorization PA = LU.
Similarly, GEwith complete pivoting computes a factorization PAQ = LU, where P and Qare permutation matrices.Exploitation of the LU factorization streamlines both the error analysisand the practical solution of linear systems. Solution of Ax = b breaks intoa factorization phase, PA = LU for partial pivoting (O(n 3 ) flops), and asubstitution phase, where the triangular systems Ly = Pb, Ux = y are solved(O(n2 ) flops). If more than one system is to be solved with the same coefficientmatrix but different right-hand sides, the factorization can be reused, with aconsequent saving in work.Computing an LU factorization A = LU is equivalent to solving the equationsIf these nonlinear equations are examined in the right order, they are easilysolved.
For added generality let(m > n) and consider an LUfactorization withand(L is lower trapezoidal: lij = 0for i < j). Suppose we know the first k - 1 columns of L and the first k - 1rows of U. Setting lkk = 1,(9.3)(9.4)We can solve for the boxed elements in the kth row of U and then the kthcolumn of L. This process is called Doolittle’s method.Algorithm 9.2 (Doolittle’s method). This algorithm computes an LU factorization A =where m > n (assuming the factorization exists),by Doolittle’s method.LU F ACTORIZATION174ANDLINEAR E QUATIONSfor k = 1: nfor j = k :n(*)endfor i = k + 1:m(**)endendCost: n 2 (m - n/3) flops.Doolittle’s met hod is mathematically equivalent to GE without pivoting,for we have, in (9.3).(9.5)and similarly for (9.3).
Had we chosen the normalization u ii1, we wouldhave obtained the Crout method. The Crout and Doolittle methods are wellsuited to calculations by hand or with a desk calculator, because they obviatethe need to store the intermediate quantitiesThey are also attractivewhen we can accumulate inner products in extended precision.It is straightforward to incorporate partial pivoting into Doolittle’s method(see. e.g.. Stewart [941, 1973, p. 138]). However, complete pivoting cannot beincorporated without changing the met hod.9.2.
Error AnalysisThe error analysis of GE is a combination of the error analyses of innerproducts and substitution. When this fact is realized. the analysis becomesstraight forward. The key observation leading to this viewpoint is that allmathematically equivalent variants of GE satisfy a common error bound. Tosee why, first note the connection between standard GE, as we have describedit, and the Doolittle met hod, as shown in (9.5). Whether the inner productin (9.5) is calculated as one operation. or whether its terms are calculatedmany operations apart, precisely the same rounding errors are sustained (assuming that extended precision accumulation of inner products is not used):all that changes is the moment when those rounding errors are committed.
Ifwe allow inner products to be reordered. so that, for example, the summation(*) in Algorithm 9.2 is calculated with the index i decreasing from k - 1 to1, instead of increasing from 1 to k - 1. then the actual rounding errors aredifferent but a common bound holds for all orderings.It suffices, then, to analyse the Doolittle method. It also suffices to analysethe met hod without pivoting, because GE with partial or complete pivotingis equivalent to GE without pivoting applied to a permuted matrix.9.2 ERROR ANALYSIS175The assignments (*) and (**) in Algorithm 9.2 are of the form y = ( c which is analysed in Lemma 8.4.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
