Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 39
Текст из файла (страница 39)
The performance of modern computers on two linearsystem benchmarks is summarized by Dongarra [312, 1995]; Dongarra’s reportis regularly updated and can be obtained from netlib under the benchmarkdirectory.Douglas [319, 1959] presented a forward error analysis for GE applied to9.6 H ISTORICAL P ERSPECTIVE189Table 9.1. Times for solution of a linear system of order n.MachineLogarithm tablesDesk computing equipmentHarvard Mark 1IBM 602 Calculating PunchPilot ACEPilot ACE’ACEEDSAC 2EDSAC 2 dYearc. 1884c.
19461947194919511954195819601960Timen7 weeks29 a2 weeks1845 minutes104 hours1017 over 3 hours1½ mins305 seconds304 seconds311007 minutesReference[952, 1994][403, 1987]b[1053, 1949][1110, 1958][1110, 1958][1110, 1958][73, 1960][73, 1960]aSymmetric positive definite system.p. 27], [507, 1948, p. 336].With magnetic drum store.dUsing magnetic tape for auxiliary storage.b[127, 1948,cdiagonally dominant tridiagonal systems arising in the solution of the heatequation by finite differences.
He concluded that the whole procedure ofsolving this partial differential equation “is stable against round-off error”. Itis surprising that Douglas’ paper is little known, because irrespective of thefact that his analysis can be simplified and clarified using modern techniques,his is one of the first truly positive rounding error results to be published.A major breakthrough in the error analysis of GE came with Wilkinson’spioneering backward error analysis, in which he proved Theorem 9.5 [1085,1961], [1088, 1963]. Apart from its simplicity and elegance and the realisticnature of the bounds, the main feature that distinguishes Wilkinson’s analysisfrom the earlier error analyses of GE is that it bounds the normwise backwarderror rather than the forward error.Wilkinson had been aware of the properties of the growth factor for partial pivoting long before developing his backward error analysis. In a 1954paper [1081, 1954] he noted thatAfter m reductions the largest element is at most 2”’ times as largeas the largest original coefficient.
It is possible to construct setsin which this factor is achieved but in practice an increase seldomtakes place; more frequently the coefficients become progressivelysmaller, particularly if the equations are ill-conditioned.This quote summarizes most of what we know today!Four of the first textbooks to incorporate Wilkinson’s analysis were those ofFox [400, 1964, pp.
161-174], Isaacson and Keller [607, 1966], Wendroff [1074,190LU FACTORIZATIONANDLINEAR E QUATIONS1966], and Forsythe and Moler [396, 1967, Chap. 21]. Fox gives a simplifiedanalysis for fixed-point arithmetic under the assumption that the growth factor is of order 1. Forsythe and Moler give a particularly readable backwarderror analysis that has been widely quoted.Wilkinson’s 1961 result is essentially the best that can be obtained bya normwise analysis. Subsequent work in error analysis for GE has mainlybeen concerned with bounding the backward error component wise, as in Theorems 9.3 and 9.4.
We note that Wilkinson could have given a componentwisebound for the backward perturbation ∆A, since most of his analysis is at theelement level.Chartres and Geuder [200, 1967] analyse the Doolittle version of GE. Theyderive a backward error result (A + ∆A) = b, with a componentwise boundo n ∆A; although they do not recognize it, their bound can be written in theform |∆A| < cnuReid [867, 1971] shows that the assumption in Wilkinson’s analysis thatpartial pivoting or complete pivoting is used is unnecessary. Without makingany assumptions on the pivoting strategy, he derives for LU factorization theresult LU = A + ∆A, |∆ aij| < 3.01 min(i - 1, j) u maxkAgain, this is acomponentwise bound. Erisman and Reid [355, 1974] note that for a sparsematrix, the term min(i - 1, j) in Reid’s bound can be replaced by mij , wheremij is the number of multiplications required in the calculation of lij (i > j)or u ij (i < j).de Boor and Pinkus [273, 1977] give the result stated in Theorem 9.4.They refer to the original 1972 German edition of [955, 198 0] for a proofof the result and explain several advantages to be gained by working with acomponentwise bound for ∆A, one of which is the strong result that ensues fortotally nonnegative matrices.
A result very similar to Theorem 9.4 is provedby Sautter [895, 1978].Skeel [919, 1979] carried out a detailed componentwise error analysis ofGE with a different flavour to the analysis given in this chapter. His aim wasto understand the numerical stability of GE (in a precisely defined sense) andto determine the proper way to scale a system by examining the behaviourof the backward and forward errors under scaling (see §9.7).
He later usedthis analysis to derive important results about fixed precision iterative refinement (see Chapter 11). Skeel’s work popularized the use of component wisebackward error analysis and componentwise perturbation theory.Thecomponentwise style of backward error analysis for GE isnow well known, as evidenced by its presence in the text books of Conte andde Boor [237, 198 0].
Golub and Van Loan [370, 198 9] (also the 1983 firstedition), and Stoer and Bulirsch [955, 1980].Forward error analyses have also been done for GE. The analyses are morecomplicated and more difficult to interpret than the backward error analyses.Olver and Wilkinson [810, 1982] derive a posteriori forward error bounds that9.7 S CALING191require the computation of A-1. Further results are given in a series of papersby Stummel [965, 1982], [966, 1985], [967, 1985], [968, 1985].Finally, probabilistic error analysis for GE is given by Barlow and Bareiss[63, 1985].9.7. ScalingPrior to solving a linear system Ax = b by GE we are at liberty to scale therows and the columns:(9.20)where D1 and D 2 are nonsingular diagonal matrices.
We apply GE to thescaled system A’y = c and then recover x from x = D2 y. Although scalingwas used in some of the earliest published programs for GE [396, 1967], [745,1962], how best to choose a scaling is still not well understood, and no singlescaling algorithm can be guaranteed always to perform satisfactorily.
Wilkinson’s remark “We cannot decide whet her equations are ill-conditioned withoutexamining the way in which the coefficients were derived” [1089, 1965, p. 198]sums up the problem of scaling rather well.The effect of scaling in GE without pivoting is easy to describe. If theelements of D 1 and D2 are powers of the machine base β (so that the scalingis done without error) and GE producesandsatisfying A + ∆A =then GE on A’ = D 1 A D2 producesandsatisfying A’ +D 1 ∆AD2 =In other words, the rounding errors in GEscale in the same way as A. This is a result of Bauer [78, 1963] (see [396,19 6 7 , Chap.
11] for a clear proof and discussion). With partial pivoting,however, the choice of pivots is affected by the row scaling (though not thecolumn scaling), and in a way that is difficult to predict.We can take a method-independent approach to scaling, by consideringany method for solving Ax = b that yields a solution satisfyingwith cn a constant. For the scaled system (9.20) we haveso it is natural to choose D 1 and D2 to minimizeAs we saw in§7.3 (Theorem 7.8), the minimum possible value is no larger than p ( |A- 1 ||A|) .However, a column scaling has the (usually) undesirable effect of changing thenorm in which the error is measured.
With row scaling only, the minimum192LU FACTORIZATIONANDLINEAR E QUATIONSvalue ofis cond(A) =achieved when D 1 A has rowsof unit 1-norm (see (7.12)). Thus row equilibration yields a cond-boundedforward error. For GE, though, it is possible to do even better. Skeel [919,-11979] shows that for D 1 = diag(|A||x|) , the forward error bound for GEPPis proportional to cond(A, x) =the catch is, of course,that the scaling depends on the unknown solution x! Row equilibration canbe regarded as approximating x by e in this “optimal” scaling.The LINPACK LU factorization routines do not include scaling, whilein the LAPACK driver routine xGESVX an initial scaling is optional. Onereason why scaling is not popular with numerical analysts is that a cond( A , x)bounded forward error and a small componentwise relative backward error areboth achieved by fixed precision iterative refinement (assuming it converges):see Chapter 11.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
