Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 12
Текст из файла (страница 12)
Good general34P RINCIPLESOFFINITE P RECISION C OMPUTATIONreferences on computational aspects of statistics are Kennedy and Gentle [649,1980] and Thisted [1000. 1988].The issues of conditioning and numerical stability play a role in any discipline in which finite precision computation is performed, but the understanding of these issues is less well developed in some disciplines than in others.In geometric computation, for example, there has been much interest sincethe late 1980s in the accuracy and robustness of geometric algorithms; seeMilenkovic [754, 1988], Hoffmann [576, 1989], and Priest [843, 1991], [844,1992 ].It was after discovering Lemma 1.1 that Wilkinson began to develop backward error analysis systematically in the 1950s.
He explains that in solvingeigenproblems Ax = x by deflation, the residual of the computed solution,(with the normalizationwas “always at noise level relative to A” [1101, 1986]. He continues, “After some years’ experience of thisI happened.
almost by accident, to observe that . . .Inother words and were exact for a matrixand sincethis meant that they were exact for a matrix differing from A at the noiselevel of the computer.“ For further details see [1101, 1986] or [1100, 1985].The numerical stability of Cramer’s rule for 2 × 2 systems has been investigated by Moler [768, 1974] and Stummel [964, 1981, §3.3].The example in §1.12.2 is taken from the HP-15C Advanced FunctionsHandbook [523, 1982], and a similar example is given by Kahan [629, 1980].For another approach to analysing this “innocuous calculation” see Problem 3.11.
The “f(2/3)” example in §1.13 is also taken from [629, 1980], inwhich Kahan states three “anti-theorems” that are included among our misconceptions in §1.19.The example (1.8) is adapted from an example of Sterbenz [938, 1974,p. 220], who devotes a section to discussing the effects of rerunning a computation at higher precision.The function expm1 := ex - 1 is provided in some floating point processorsand mathematics libraries as a more accurate alternative to forming ex andsubtracting 1 [991, 1992]. It is important in the computation of sinh and tanh,for example (since sinh x = e-x(e2 x - 1)/2). Algorithm 2 in §1.14.1 is due toKahan [629, 1980].The instability.
and stability of GE without pivoting applied to an upperHessenberg matrix (§1.16) was first pointed out and explained by Wilkinson [1084, 1960]; Parlett [818, 1965] also gives a lucid discussion. In the 1950sand 1960s. prior to the development of the QR algorithm, various methodswere proposed for the nonsymmetric eigenvalue problem that involved transforming a matrix to Hessenberg form H and then finding the zeros of the characteristic polynomial det(H The most successful method of this typewas Laguerre’s iteration. described by Parlett [817, 196 4], and used in conjunction with Hyman’s method for evaluating det(H Hyman’s met hod1.21 N OTESANDR EFERENCES35is described in §13.5.1.Classic papers dispensing good advice on the dangers inherent in numerical computation are the “pitfalls” papers by Stegun and Abramowitz [937,1956] and Forsythe [394, 1970].
The book Numerical Methods That Work byActon [4, 1970] must also be mentioned as a fount of hard-earned practicaladvice on numerical computation (look carefully and you will see that thefront cover includes a faint image of the word “Usually” before “Work”). Ifit is not obvious to you that the equation x2 - 10x + 1 = 0 is best thought ofas a nearly linear equation for the smaller root, you will benefit from readingActon (see p.
58). Everyone should read Acton’s “Interlude: What Not ToCompute” (pp. 245-257).Finally, we mention the paper “How to Get Meaningless Answers in Scientific Computation (and What to Do About it)” by Fox [401, 1971]. Fox, acontemporary of Wilkinson, founded the Oxford Computing Laboratory andwas for many years Professor of Numerical Analysis at Oxford. In this paperhe gives numerous examples in which incorrect answers are obtained fromplausible numerical methods (many of the examples involve truncation errorsas well as rounding errors). The section titles provide a list of reasons whyyou might compute worthless answers:• Your problem might be ill conditioned.• Your method might be unstable.• You expect too much “analysis” from the computers.• Your intuition fails you.• You accept consistency too easily.• A successful method may fail in slightly different circumstances.• Your test examples may be too special.Fox estimates [401, 1971, p.
296] that “about 80 per cent of all the resultsprinted from the computer are in error to a much greater extent than the userwould believe.”8This reason refers to using an inappropriate convergence test in an iterative process.36P RINCIPLESOFFINITE P RECISION C OMPUTATIONProblemsThe road to wisdom?Well, it’s plain and simple to express:Errand errand err againbut lessand lessand less.-PIET HEIN, Grooks (1966)1.1. In error analysis it is sometimes convenient to boundinstead ofObtain inequalities betweenand1.2. (Skeel and Keiper [923, 1993, §1.2]) The number y =was evaluated at t digit precision for several values of t, yielding the values shown in thefollowing table.
which are in error in at most one unit in the least significantdigit (the first two values are padded with trailing zeros):t1015202530y2625374 12600000000262537412640769000262537412640768744.00262537412640768744.0000000262537412640768743.999999999999Does it follow that the last digit before the decimal point is 4?1.3. Show how to rewrite the following expressions to avoid cancellation forthe indicated arguments.1.4. Give stable formulae for computing the square root x + iy of a complexnumber a + ib.1.5. [523, 1982] By writing (1 + 1/n) n = exp(nlog(1 + 1/n)), show how tocompute (1 + 1/n) n accurately for large n. Assume that the log function iscomputed with a relative error not exceeding u. (Hint: adapt the techniqueused in §1.14.1.)37P ROBLEMS1.6.
(Smith [928, 1975]) Type the following numbers into your pocket calculator, and look at them upside down (you or the calculator):07734380793188083500757738.57734 × 10403331The famous “_world” programObjectNameAdjectiveExclamation on finding a bugA high quality floating point arithmeticFallen tree trunks1.7.
A condition number for the sample variance (1.4), here denoted by V(x) :can be defined byShow thatThis condition number measures perturbations in x componentwise. A corresponding normwise condition number isShow that1.8. (Kahan, Muller, [781,the recurrence1989],Francois and Muller [406,x k+1 = 111 - (1130 - 3000/x k- 1 )/x k,x0 = 11/2,1991])Considerx 1 = 61/11.In exact arithmetic the xk form a monotonically increasing sequence that converges to 6.
Implement the recurrence on your computer or pocket calculatorand compare the computed x 34 with the true value 5.998 (to four correctsignificant figures). Explain what you see.The following questions require knowledge of material from later chapters.1.9. Cramer’s rule solves a 2 × 2 system Ax = b according tod = a 1 1 a 22 - a2 1 a 1 2 ,x1 = (b 1 a22 - b2 a 12)/d ,x2 = (a 11b2 - a 2 1 b 1 )/d.38P RINCIPLESOFFINITE P RECISION C OMPUTATIONShow that, assuming d is computed exactly (this assumption has little effecton the final bounds), the computed solutionsatisfieswhere γ3 = 3u/(1-3u), cond(A,x) =and cond(A) =This forward error bound is as small as that for a backwardstable method (see §7.2, §7.6) so Cramer’s rule is forward stable for 2 × 2systems.1.10. Show that the computed sample variancethe two-pass formula (1.4) satisfies= fl(V(x)) produced by(Note that this error bound does not involve the condition numbers KC orfrom Problem 1.7, at least in the first-order term.
This is a rare instanceof an algorithm that determines the answer more accurately than the datawarrants!)KNPreviousHomeChapter 2Floating Point ArithmeticFrom 1946-1948 a great deal of quite detailed coding was done.The subroutines for floating-point arithmetic were . . .produced by Alway and myself in 1947 . . .They were almost certainly the earliest floating-point subroutines.-J. H. WILKINSON, Turing’s Work at theNational Physical Laboratory .