Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 40
Текст из файла (страница 40)
Even Skeel’s optimal scaling does not guarantee a smallcomponentwise relative backward error.Some programs for GEPP incorporate row scaling implicitly. They compute row scale factors d 1. . . . , dn . but. instead of applying GEPP to diag(d i )-1 ×A, they apply it to A and choose as pivot row at the kth stage a row r forwhichis maximal. This type of scaling has the sole effect of influencing the choice of pivots.
There is little justification for using it, and the bestbound for the growth factor is 2n -1 multiplied by a product of terms di1 /d i2that can be large.There is, however. one situation in which a form of implicit row scalingis beneficial. Consider the pivoting strategy that selects as the kth pivot anelementfor which(9.21)A result of Peña [825, 1994] shows that if there exists a permutation matrix Psuch that PA has an LU factorization PA = LU with |PA| = |L||U|, then sucha factorization will be produced by the pivoting scheme (9.21). This meansthat, unlike for partial pivoting, we can use the pivoting scheme (9.21) withimpunity on totally nonnegative matrices and their inverses, row permutationsof such matrices, and any matrix for which some row permutation has the“|PA| = |L||U|” property.
However, this pivoting strategy is as expensive ascomplete pivoting to implement, and for general A it is not guaranteed toproduce a factorization as stable as that produced by partial pivoting.9.8. A Posteriori Stability TestsHaving solved a linear system by LU factorization we can compute the componentwise or normwise backward error at the cost of evaluating one or twomatrix-vector products (see Theorems 7.1 and 7.3). In some situations,9.8 A P OSTERIORI S TABILITY T ESTS193though, we may wish to assess the stability of a computed LU factorization before using it to solve one or more linear systems.
One possibility isto compute the growth factor by monitoring the size of elements during theelimination, at a cost of O(n 3) comparisons. This has been regarded as ratherexpensive, and more efficient ways to estimate pn have been sought.Businger [169, 1971] describes a way to obtain an upper bound for p n inO(n2) operations. This approach is generalized by Erisman and Reid [355,1974], who apply the Holder inequality to the equationto obtain the bound(9.22)where p -1 + q-1 = 1. In practice. p = 1,2,are the values of interest.Barlow [56, 1986] notes that application of the Holder inequality instead toyields a sometimes sharper bound.It is interesting to note that in light of experience with the bound (9.22),Reid [868, 198 7] recommends computing the growth factor explicitly in thecontext of sparse matrices, arguing that the expense is justified because (9.22)can be a very weak bound.
See Erisman et al. [354, 1987] for some empiricalresults on the quality of the bound.Chu and George [209, 1985] observe that the-norm of the matrixcan be computed in O(n 2) operations without forming the matrix explicitly,sinceThus one can cheaply compute a bound onfrom the componentwisebackward error bounds in (9.6) and (9.7).All the methods discussed in this section make use of an a priori erroranalysis to compute bounds on the backward error. Because the bounds do nottake into account the statistical distribution of rounding errors, and becausethey involve somewhat pessimistic constant terms, they cannot be expectedto be very sharp. Thus it is important not to forget that it is straightforwardto compute the backward error itself: A Exact computation costs aprohibitively expensive O(n 3) operations, butcan be estimated in194LU FACTORIZATIONANDLINEAR E QUATIONSO(n 2) operations using the matrix norm estimator in Algorithm 14.4.
Anotherpossibility is to use a running error analysis, in which an error bound iscomputed concurrently with the factors (see 53.3).9.9. Sensitivity of the LU FactorizationAlthough Theorem 9.3 bounds the backward error of the computed LU factorsL and U, it does not give any indication about the size of the forward errorsL and U For most applications of the LU factorization it is thebackward error and not the forward errors that matters, but it is still of someinterest to know how big the forward errors can be. This is a quest ion ofperturbation theory and is answered by the next result.Theorem 9.14 (Barrlund, Sun). Let the nonsingular matricesand A+∆A have LU factorizations A = LU and A+∆ A = (L+∆L)(U+∆ U),and assume that ||G||2 < 1, where G = L- 1 ∆AU -1.
Then(9.23)Moreover, if< 1, where= (L + ∆L)- 1 ∆A(U + ∆U)-1, thenwhere stril(·) and triu(·) denote, respectively, the strictly lower triangular partand the upper triangular part of their matrix arguments.The normwise bounds (9.23) imply that x(A) := ||L - 1 ||2 ||U- 1 ||2 ||A|| 2 isan upper bound for the condition numbers of L and U under normwise perturbations. We haveand the ratio x(A)/κ 2 (A) can be arbitrarily large (though if partial pivotingis used then κ2 (L) < n2n - 1 ).The componentwise bounds in Theorem 9.14 are a little unsatisfactory inthat they involve the unknown mat rices ∆L and ∆U, but we can set theseterms to zero and obtain a bound correct to first order.9.10 N OTESANDR EFERENCES1959.10.
Notes and ReferencesA variant of GE was used by the Chinese around the first century AD; the JiuZhang Suanshu (Nine Chapters of the Mathematical Art) contains a workedexample for a system of five equations in five unknowns [619, 1991, pp. 156-177].[696, 1989]Gauss, who was a great statistician and numerical analyst, developed hiselimination method as a tool to help him prove results in linear regressiontheory. The first published appearance of GE is in his Theoria Motus (1809).Stewart [952, 1994] gives a survey of Gauss’s work on solving linear systems:see also the afterword in [423, 1995].The traditional form of GE, as given at the start of this chapter.
can beexpressed algorithmically asfor k = 1:nfor j = k+ 1 :nfor i = k + 1:na i j = a i j - (a i k /a k k )a k jendendendThis is identified as the kji form of GE. Altogether there are six possibleorderings of the three loops. Doolittle’s method (Algorithm 9.2) is the ijkor jik variant of Gaussian elimination. The choice of loop ordering does notaffect the stability of the computation, but can greatly affect the efficiency ofGE on a high performance computer.
For more on the different loop orderingsof GE see Chapter 12; Dongarra, Gustavson. and Karp [310, 1984]; and thebooks by Dongarra, Duff, Sorensen, and van der Vorst [315, 1991] and Goluband Van Loan [470, 1989].This chapter draws on the survey paper Higham [545, 1990]. Theorems 9.6and 9.7 are from Higham and Higham [562, 1989].Myrick Hascall Doolittle (1830-1913) was a “computer of the United Statescoast and geodetic survey” [362, 1987]. Crout’s method was published in anengineering journal in 1941 [255, 1941].GE and its variants were known by various descriptive names in the earlydays of computing. These include the bordering met hod, the escalator met hod(for matrix inversion), the square root method (Cholesky factorization), andpivotal condensation. A good source for details of these methods is Faddeeva [360, 1959].In a confidential 1948 report that “covers the general principles of both thedesign of the [Automatic Computing Engine] and the method of programmingadopted for it”, Wilkinson gives a program implementing GE with partialpivoting and iterative refinement [1080, 1948, p.
111]. This was probably the196LU F ACTORIZATIONANDLINEAR E QUATIONSfirst such program to be written and for a machine that had not yet beenbuilt!The terms “partial pivoting“ and “complete pivoting” were introduced byWilkinson in [1085, 1961]. The pivoting techniques themselves were in usein the 1940s and it is not clear who, if anyone, can be said to have inventedthem: indeed, von Neumann and Goldstine [1057, 1947, §4.2] refer to completepivoting as the “customary procedure”.There is a long history of published programs for GE.
beginning with Croutroutines of Forsythe [390, 1960], Thacher [999, 1961], McKeeman [745, 1962],and Bowdler, Martin, Peters, and Wilkinson [138, 1966], all written in Algol60 (which was the “official” language for publishing mathematical software inthe 1960s. and a strong competitor to Fort ran for practical use at that time).The GE routines in LAPACK are the latest in a lineage beginning with theFortran routines decomp and solve in Forsythe and Moler‘s book [396, 1967],and continuing with routines by Moler [766, 1972], [767, 1972] (which achieveimproved efficiency in Fortran by accessing arrays by column rather than byrow), Forsythe, Malcolm, and Moler [395, 1977] (these routines incorporatecondition estimation-see Chapter 14), and LINPACK [307, 1979].LU factorization of totally nonnegative matrices has been investigated byCryer [257, 1973], [258, 1976], Ando [21, 198 7], and de Boor and Pinkus[273, 1977].
It is natural to ask whether we can test for total nonnegativitywithout computing all the minors. The answer is yes: for an n × n matrixtotal nonnegativity can be tested in O(n 3) operations. as shown by Gascaand Peña [421, 1992]. The test involves carrying out a modified form of GEin which all the elimination operations are between adjacent rows and thenchecking whether certain pivots are positive. Note the analogy with positivedefiniteness, which holds for a symmetric matrix if and only if all the pivotsin GE are positive.The dilemma of whether to define the growth factor in terms of exact orcomputed quantities is faced by all authors; most make one choice or the other,and go on to derive, without comment, bounds that are strictly incorrect.Theorem 9.8, for example, bounds the exact growth factor; the computed onecould.
conceivably violate the bound, but only by a tiny relative amount. vanVeldhuizen [1045, 1977] shows that for a variation of partial pivoting thatallows either a row or column interchange at each stage, the growth factordefined in terms of computed quantities is at most about (1 + 3 nu)2n - 1 ,compared with the bound 2n -1 for the exact growth factor.The idea of deriving error bounds for GE by analysing the equations obtained by solving A = LU is exploited by Wilkinson [1097, 1974], who gives ageneral analysis that includes Cholesky factorization. This paper gives a concise summary of error analysis of factorization methods for linear equationsand least squares problems.Various authors have tabulated growth factors in extensive tests with ran-9.10 N OTESANDR EFERENCES197dom matrices.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
