Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 22
Текст из файла (страница 22)
113]).A straight forward notation for rounding errors that is subsumed by thenotation described in this chapter is suggested by Scherer and Zeller [899,1980].Ziv [1131, 1982] proposes the relative error measurefor vectors x and y and explains some of its favourable properties for erroranalysis.Wilkinson [1089, 1965, p. 447] gives error bounds for complex arithmetic;Olver [808, 198 3] does the same in the relative precision framework. Demmel [280, 1984] gives error bounds that extend those in Lemma 3.5 by takinginto account the possibility of underflow.Henrici [521, 1980] gives a brief, easy to read introduction to the use ofthe model (2.4) for analysing the propagation of rounding errors in a generalalgorithm. He uses a set notation that is another possible notation to add tothose in §3.4.The perspective on error analysis in §3.8 was suggested by J.
W. Demmel.Problems3.1. Prove Lemma 3.1.3.2. (Kielbasiniski and Schwetlick [658, 1988], [659, 1992]) Show that if pi 1in Lemma 3.1 then the stronger bound |φ n | < nu/(1-½nu) holds for nu < 2.P ROBLEMS853.3. One algorithm for evaluating a continued fractionisq n+1 = an +1for k = n:-1:0qk = ak + bk/qk+end1Derive a running error bound for this algorithm.3.4.
Prove Lemma 3.3.3.5. (Backward error result for matrix multiplication.) Let AIRn×n andn×nBIRboth be nonsingular. Show that fl(AB) = (A + ∆A)B, where| ∆A| < γn |A||B||B -1|, and derive a corresponding bound in which B is perturbed.3.6. (Backward error definition for matrix multiplication.) Let AIRm × nn×pand BIRbe of full rank and suppose CAB. Define the componentwise backward errorwhere E and F have nonnegative entries. Show thatwhere R = C - AB and G = EF.
Explain why the definition of w makessense only when A and B have full rank. Define a mixed backward/forwarderror applicable to the general case.3.7. Give analogues of the backward error results (3.4) and (3.10) for complexx, y, and A.3.8. Let Al,. . . ,AkIRn × n .
Show that3.9. Which is the more accurate way to compute x2 - y2: as x2 - y2 or as(x+y)(x-y)?(Assume the use of a guard digit. Note that this computationarises when squaring a complex number.)86BASICS3.10. Prove Lemma 3.6.3.11. (Kahan [629, 1980]) Consider this MATLAB function, which returns theabsolute value of its first argument xIRn :function z = absolute(x , m )y = x.-2;for i=l:my = sqrt(y);endz = y;for i=l:m-1z = z.
-2;endHere is some output from a 486DX workstation that uses IEEE standarddouble precision arithmetic:> > x = [ . 2 5 . 5 . 7 5 1 . 2 5 1 . 5 2 ] ; z = a b s o l u t e ( x , 5 0 ) ; [ x ; z]ans =0.25000.50000.75001.25001.50002.00000.25280.50280.77881.28401.45502.1170Give an error analysis to explain the results.The same machine produced this output:> > x = [ . 2 5 .
5 . 7 5 1 . 2 5 1 . 5 2 ] ; z = a b s o l u t e ( x , 7 5 ) ; [ x ; z]ans =0.25000.50000.75001.2500 1.50002.00001.00001.00001.00001.0000 1.00001.0000But a Sun SPARCstation, which also uses IEEE standard double precisionarithmetic, producedans =0.25000.50000.75001.2500 1.50002.00000001.0000 1.00001.0000Explain these results and why they differ (cf. §1.12.2).3.12. Consider the quadrature rulewhere the weights wi and nodes xi are assumed to be floating point numbers.Assuming that the sum is evaluated in left-to-right order and thatobtain and interpret a bound forwherePrevious Home NextChapter 4SummationI do hate sums.There is no greater mistake than to call arithmetic an exact science.There are .
. . hidden laws of Numberwhich it requires a mind like mine to perceive.For instance, if you add a sum from the bottom up,and then again from the top down,the result is always different.-MRS. LA TOUCHE10Joseph Fourier introduced this delimited z-notation in 1820,and it soon took the mathematical world by storm.-RONALD L. GRAHAM, DONALD E. KNUTH, andOREN PATASHNIK, Concrete Mathematics (1989)One of the major difficulties in a practical [error] analysisis that of description.An ounce of analysis follows a pound of preparation.-BERESFORD N. PARLETT, Matrix Eigenvalue Problems (1965)10Quoted in Mathematical Gazette [730,1924 ].8788SUMMATIONSums of floating point numbers are ubiquitous in scientific computing.
Theyoccur when evaluating inner products, means, variances, norms, and all kindsof nonlinear functions. Although at first sight summation might appear tooffer little scope for algorithmic ingenuity, the usual “recursive summation”(with various orderings) is just one of a variety of possible techniques. Wedescribe several summation methods and their error analyses in this chapter. No one method is uniformly more accurate than the others, but someguidelines can be given on the choice of method in particular cases.4.1. Summation MethodsIn most circumstances in scientific computing we would naturally translate asuminto code of the forms =0for i = 1:ns=s+ xendiThis is known as recursive summation.
Since the individual rounding errorsdepend on the operands being summed, the accuracy of the computed sumvaries with the ordering of the xi . (Hence Mrs. La Touche, quoted at thebeginning of the chapter, was correct if we interpret her remarks as applyingto floating point arithmetic.) Two interesting choices of ordering are theincreasing order |x1| < | x 2 | < . . .
< |x n |, and the decreasing order |x 1| >|x2| > . . . > |xn|.Another method is pairwise summation (also known as cascade summation, or fan-in summation), in which the xi are summed in pairs accordingtoand this pairwise summation process is repeated recursively on the yi, i =1 : [ ( n+1)/2]. The sum is obtained in [log2 n] stages. For n = 6, for example,pairwise summation formsS 6 = (( x 1 + x 2 ) + ( x 3 + x 4 ) ) + ( x 5 + x 6 ) .Pairwise summation is attractive for parallel computing, because each of the[log2 n] stages can be done in parallel [573, 1988, §5.2.2].A third summation method is the insertion met hod.
First, the x i aresorted by order of increasing magnitude (alternatively, some other orderingcould be used). Then x1 + x2 is formed, and the sum is inserted into thel i s t x2,. . . ,x n , maintaining the increasing order. The process is repeatedrecursively until the final sum is obtained.
In particular cases the insertion4.2 ERROR ANALYSIS89method reduces to one of the other two. For example, if x i = 2 i -1, theinsertion method is equivalent to recursive summation, since the insertion isalways to the bottom of the list:12483487815.On the other hand, if 1 < x1 < x2 < ...
< xn < 2, every insertion is to theend of the list, and the method is equivalent to pairwise summation if n is apower of 2; for example, if 0 < < 1/2,To choose between these met hods we need error analysis, which we developin the next section.4.2. Error AnalysisError analysis can be done individually for the recursive, pairwise and insertion summation methods, but it is more profitable to recognize that each is aspecial case of a general algorithm and to analyse that algorithm.Algorithm 4.1. Given numbers x1,. . .
,xn this algorithm computes Sn=L e t S = {x 1,...,xn }.repeat while S contains more than one elementRemove two numbers x and y from S,and add their sum x + y to S.endAssign the remaining element of S to S n .Note that since there are n numbers to be added and hence n - 1 additionsto be performed, there must be precisely n - 1 executions of the repeat loop.First, let us check that the previous methods are special cases of Algorithm 4.1.
Recursive summation (with any desired ordering) is obtained bytaking x at each stage to be the sum computed on the previous stage of thealgorithm. Pairwise summation is obtained by [log2 n] groups of executionsof the repeat loop, in each group of which the members of S are broken intopairs, each of which is summed. Finally, the insertion method is, by definition,a special case of Algorithm 4.1.Now for the error analysis.
Express the ith execution of the repeat loopas Ti = xi1 + yil. The computed sums satisfy (using (2.5))(4.1)90SUMMATIONThe local error introduced in formingThe overall error is the sumof the local errors (since summation is a linear process), so overall we have(4.2)The smallest possible error bound is therefore(4.3)(This is actually in the form of a running error bound, because it contains thecomputed quantities-see §3.3.) It is easy to see thatfor each i, and so we have also the weaker bound(4.4)This is a forward error bound. A backward error result showing thatisthe exact sum of termswithcan be deduced from (4.1),using the fact that no number xi takes part in more than n - 1 additions.The following criterion is apparent from (4.2) and (-1.3):In designing or choosing a summation met hod to achieve high accuracy, the aim should be to minimize the absolute values of theintermediate sums Ti .The aim specified in this criterion is surprisingly simple to state.
When weconsider specific methods, however, we find that the aim is difficult to achieve.Consider recursive summation, for whichIdeally,we would like to choose the ordering of the xi to minimizeThisis a combinatorial optimization problem that is too expensive to solve in thecontext of summation. A reasonable compromise is to determine the orderingsequentially by minimizing, in turn, |x 1|, |S 2|, . . . , |S n-1|. This orderingstrategy, which we denote by Psum, can be implemented with O(n log n)comparisons. If we are willing to give up the property that the ordering isinfluenced by the signs of the xi we can instead use the increasing ordering,which in general will lead to a larger value ofthan that for thePsum ordering.
If all the x i have the same sign then all these orderingsare equivalent. Therefore when summing nonnegative numbers by recursivesummation the increasing ordering is the best ordering, in the sense of havingthe smallest a priori forward error bound.4.2 ERROR ANALYSIS91How does the decreasing ordering fit into the picture? For the summationof positive numbers this ordering has little to recommend it. The bound (4.3)is no smaller, and potentially much larger, than it is for the increasing ordering.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
