Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 24
Текст из файла (страница 24)
With compensated summation the errorsen are much less affected by rounding errors and do not grow in the range ofn shown (for n = 108, en is about 10 times larger than it would be in exactarithmetic). Plots of U-shaped curves showing total error against stepsizeare common in numerical analysis textbooks (see, e.g., Forsythe, Malcolm,and Moler [395, 1977, p. 119] and Shampine [910, 1994, p. 259]), but thetextbooks rarely point out that the “U” can be flattened out by compensatedsummation.96Figure 4.2. Errors |y(1) compensated summation.SUMMATIONfor Euler’s method with (“×”) and without (“o”)The cost of applying compensated summation in an ordinary differentialequation solver is almost negligible if the function f is at all expensive to evaluate.
But, of course, the benefits it brings are noticeable only when a vastnumber of integration steps are taken. Very long-term integrations are undertaken in celestial mechanics, where roundoff can affect the ability to trackplanetary orbits.
Researchers in astronomy use compensated summation, andother techniques, to combat roundoff. An example application is a 3 millionyear integration of the planets in the solar system by Quinn, Tremaine, andDuncan [854, 1991]; it used a linear multistep method of order 13 with aconstant stepsize of 0.75 days and took 65 days of machine time on a SiliconGraphics 4D-25 workstation. See also Quinn and Tremaine [853, 1991] andQuinlan [851, 1994].Finally, we describe an even more ingenious algorithm called doubly compensated summation, derived by Priest [844, 1992] from a related algorithmof Kahan.
It is compensated summation with 2 extra applications of thecorrection process11 and it requires 10 instead of 4 additions per step. Thealgorithm is tantamount to simulating double precision arithmetic with single precision arithmetic; it requires that the summands first be sorted into11The algorithm should perhaps be called triply compensated summation, but we adoptPriest’s terminology.4.4 O THER S UMMATION M ETHODS97decreasing order, which removes the need for certain logical tests that wouldotherwise be necessary.Algorithm 4.3 (doubly compensated summation). Given floating point numbers x 1,.
. . , xn , this algorithm forms the sumby doubly compensated summation. All expressions should be evaluated in the order specifiedby the parentheses.Sort the xi so that |x i| > |x 2| > . . . > |x n |.s 1 = x1; c1 = 0for i = 2:nyk = ck-1 + xku k = xk - (yk - ck- 1 )tk = yk + sk- 1υ k = yk - (tk - sk-l )zk = u k + υ ksk = tk + zkck = zk - (sk - tk)endPriest [844, 1992, §4.1] analyses this algorithm for t-digit base β arithmeticthat satisfies certain reasonable assumptions-ones which are all satisfied byIEEE arithmetic.
He shows that if n < β t -3 then the computed sumsatisfiesthat is, the computed sum is accurate virtually to full precision.4.4. Other Summation MethodsWe mention briefly two further classes of summation algorithms. The firstbuilds the sum in a series of accumulators, which are themselves added to givethe sum. As originally described by Wolfe [1107, 1964] each accumulator holdsa partial sum lying in a different interval. Each term xi is added to the lowestlevel accumulator; if that accumulator overflows it is added to the next-highestone and then reset to zero, and this cascade continues until no overflow occurs.Modifications of Wolfe’s algorithm are presented by Malcolm [723, 1971] andRoss [880, 1965].
Malcolm [723, 1971] gives a detailed error analysis to showthat his method achieves a relative error of order u. A drawback of thealgorithm is that it is strongly machine dependent. An interesting and crucialfeature of Malcolm’s algorithm is that on the final step the accumulatorsare summed by recursive summation in order of decreasing absolute value,which in this particular situation precludes severe loss of significant digitsand guarantees a small relative error.98SUMMATIONAnother class of algorithms, referred to as “distillation algorithms” byKahan [633, 1987], work as follows: given xi = fl( x i ), i = 1:n, they iterativelyconstruct floating point numberssuch thatterminating whenapproximateswith relative error at most u.Kahan states that these algorithms appear to have average run times of orderat least n log n.
See Bohlender [130, 1977], Kahan [633, 1987], Leuprecht andOberaigner [700, 1982], Pichat [830, 1972], and Priest [844, 1992, pp. 66-69)for further details and references.4.5. Statistical Estimates of AccuracyThe rounding error bounds presented above can be very pessimistic, becausethey account for the worst-case propagation of errors.
An alternative wayto compare summation methods is through statistical estimates of the error,which may be more representative of the average case. A statistical analysisof three summation methods has been given by Robertazzi and Schwartz [874,1988] for the case of nonnegative xi . They assume that the relative errors infloating point addition are statistically independent and have zero mean andfinite varianceTwo distributions of nonnegative xi are considered: theuniform distribution on [0,2µ ].
and the exponential distribution with meanµ. Making various simplifying assumptions Robertazzi and Schwartz estimate the mean square error (that is, the variance of the absolute error) ofthe computed sums from recursive summation with random, increasing anddecreasing orderings, and from insertion summation and pair-wise summation(with the increasing ordering).
Their results for the summation of n numbersare given in Table -1.1.The results show that for recursive summation the ordering affects onlythe constant in the mean square error, with the increasing ordering havingthe smallest constant and the decreasing ordering the largest; since the xi arenonnegative, this is precisely the ranking given by the rounding error bound(4.3). The insertion and pairwise summation methods have mean squareerrors proportional to n 2 rather than n 3 for recursive summation, and theinsertion method has a smaller constant than pairwise summation. This is alsoconsistent with the rounding error analysis, in which for nonnegative xi theinsertion met hod satisfies an error bound no larger than pairwise summationand the latter method has an error bound with a smaller constant than forrecursive summation (log2 n versus n).4.6.
Choice of MethodThere is a wide variety of summation methods to choose from. For eachmet hod the error can vary greatly with the data, within the freedom afforded4.6 C HOICEOFM ETHOD99Table 4.1. Mean square errors for nonnegative xi .Distrib.Unif(0,2µ)EXP(µ)Increasing0.20µ 2 n 3 σ20.13µ 2 n 3 σ2Random0.33µ 2 n 3 σ20.33µ 2 n 3 σ2Decreasing0.53µ 2 n 3 σ20.63µ 2 n 3 σ2Insertion2.6µ 2 n 2 σ22.6µ 2 n 2 σ2Pairwise2.7µ 2 n 2 σ24.0µ 2 n 2 σ2by the error bounds; numerical experiments show that, given any two of themethods, data can be found for which either method is more accurate thanthe other [553, 1993]. However, some specific advice on the choice of met hodcan be given.1. If high accuracy is important, consider implementing recursive summation in higher precision; if feasible this may be less expensive (and moreaccurate) than using one of the alternative methods at the working precision.
What can be said about the accuracy of the sum computed athigher precision? If Sn =xi is computed by recursive summationat double precision (unit roundoff u 2) and then rounded to single precision, an error bound of the formholds.Hence a relative error of order u is guaranteed if nuPriest [844, 1992, pp. 62-63] shows that if the xi are sorted in decreasing order of magnitude before being summed in double precision, thenholds provided only that n < β t -3 for t-digit baseβ arithmetic satisfying certain reasonable assumptions.
Therefore thedecreasing ordering may be worth arranging if there is a lot of cancellation in the sum. An alternative to extra precision computation is doublycompensated summation, which is the only other method described herethat guarantees a small relative error in the computed sum.2. For most of the methods the errors are, in the worst case, proportionalto n. If n is very large, pairwise summation (error constant log2 n) andcompensated summation (error constant of order 1) are attractive.3. If the xi all have the same sign then all the methods yield a relative errorof at most nu and compensated summation guarantees perfect relativeaccuracy (as long as nu < 1). For recursive summation of one-signeddata, the increasing ordering has the smallest error bound (4.3) andthe insertion method minimizes this error bound over all instances ofAlgorithm 4.1.4.
For sums with heavy cancellationrecursivesummation with the decreasing ordering is attractive, although it cannotbe guaranteed to achieve the best accuracy.SUMMATION100Considerations of computational cost and the way in which the data aregenerated may rule out some of the met hods. Recursive summation in thenatural order, pairwise summation, and compensated summation can be implemented in O(n) operations for general xi , but the other methods are moreexpensive since they require searching or sorting.
Furthermore, in an application such as the numerical solution of ordinary differential equations, wherexk is not known untilxi has been formed. sorting and searching maybe impossible.4.7. Notes and ReferencesThis chapter is based on Higham [553, 1993].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
