Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 19
Текст из файла (страница 19)
What is 00 in IEEE arithmetic?2.16. Evaluate these expressions in any IEEE arithmetic environment available to you. Are the values returned what you would expect? (None of theresults is specified by the IEEE standard.)2.17. In the course of solving ax2 - 2bx + c = 0 for x, the expressionmust be computed. Can the true value of b 2 - ac be nonnegative and yet itscomputed value be negative?2.18.
Can Theorem 2.4 be strengthened to say that fl (x - y) is computedexactly whenever the exponents of x > 0 and y > 0 differ by at most l?2.19. Two requirements that we might ask of a routine for computingfloating point arithmetic are that the identitiessatisfied. Which, if either, of these is a reasonable requirement?in2.20. Are there any floating point values of x and y (excepting values both0, or so huge or tiny to cause overflow or underflow) for which the computedvalue ofexceeds 1?2.21. (Kahan) A natural way to compute the maximum of two numbers xand y is with the code% max(x, y)if x > y thenm a x =xelsemax = yendDoes this code always produce the expected answer in IEEE arithmetic?2.22. Prove that Kahan’s formula (2.7) computes the area of a triangle accurately if a guard digit is used in subtraction.
(Hint: you will need oneinvocation of Theorem 2.5.)65P ROBLEMS2.23. (Kahan) Describe the result of the computation y = (x + x) - x on abinary machine with a guard digit and one without a guard digit.2.24. (Kahan) Let f(x) = (((x - 0.5) + x) - 0.5) + x. Show that if f isevaluated as shown in single or double precision binary IEEE arithmetic thenf (x)0 for all floating point numbers x.2.25. (Kahan) Consider a machine that can perform a fused multiply-addoperation with just a single rounding error:f l(x + y * z) = (x + y * z)(1 + δ),| δ | < u.Show that, on such a machine, the algorithmw = bce = w - b* cx = (a * d - w) + ecomputes x =with high relative accuracy.2.26.
Derive Newton’s method for solving f(x) = a - 1/x = 0. This methodwas used on early computers (and is still used on some Cray computers, forexample) to implement reciprocation in terms of multiplication and thencedivision as a/b = a * (1/b); see, e.g., [506, 1946].2.27. Suppose we have an iterative algorithm for computing z = x/y. Derivea convergence test that terminates the iteration (only) when full accuracy hasbeen achieved. Assume the use of IEEE arithmetic with gradual underflow(use (2.8)).Previous HomeChapter 3BasicsA method of inverting the problem of round-off error is proposedwhich we plan to employ in other contexts andwhich suggests that it may be unwise toseparate the estimation of round-off errorfrom that due to observation and truncation.-WALLACE J.
GIVENS, Numerical Computation of theCharacteristic Values of a Real Symmetric Matrix (1954)The enjoyment of one’s tools is an essential ingredient of successful work.-DONALD E. KNUTH, The Art of Computer Programming,Volume 2, Seminumerical Algorithms (1981)The subject of propagation of rounding error,while of undisputed importance in numerical analysis,is notorious for the difficulties which it presents when it is to betaught in the classroom in such a manner that the student isneither insulted by lack of mathematical contentnor bored by lack of transparence and clarity.-PETER HENRICI, A Model for the Propagationof Rounding Error in Floating Arithmetic (1980)The two main classes of rounding error analysis are not,as my audience might imagine, ‘backwards’ and ‘forwards’,but rather ‘one’s own’ and ‘other people’s’.One’s own is, of course, a model of lucidity;that of others serves only to obscure theessential simplicity of the matter in hand.-J.
H. WILKINSON, The State of the Art in Error Analysis (1985)67Next68B ASICSHaving defined a model for floating point arithmetic in the last chapter, wenow apply the model to some basic matrix computations, beginning with innerproducts. This first application is simple enough to permit a short analysis,yet rich enough to illustrate the ideas of forward and backward error. It alsoraises the thorny question of what is the best notation to use in an erroranalysis.
We introduce the “γ n ” notation, which we use widely, though notexclusively, in the book. The inner product analysis leads immediately toresults for matrix-vector and matrix-matrix multiplication.In the last two sections we determine a model for rounding errors in complex arithmetic and derive some miscellaneous results of use in later chapters.3.1. Inner and Outer ProductsConsider the inner product sn = xT y, where x, yIRn . Since the order ofevaluation of sn = x1 y1 + .
. . + xnyn affects the analysis (but not the final errorbounds), we will assume that the evaluation is from left to right. (The effectof particular orderings is discussed in detail in Chapter 4, which considersthe special case of summation.) In the following analysis, and throughout thebook, a hat denotes a computed quantity.Let si = x1 y1 + . . . + xi yi denote the ith partial sum. Using the standardmodel (2.4), we have(3.1)where |δi | < u, i = 1:3.
For our purposes it is not necessary to distinguishbetween the different δi terms, so to simplify the expressions let us drop thesubscripts on the δi and write 1 + δi1 ± δ. ThenThe pattern is clear. Overall, we have(3.2)There are various ways to simplify this expression. A particularly elegant wayis to use the following result.3.1 INNERANDOUTER PRODUCTS69Lemma 3.1. If |δi | < u and pi = ±1 for i = 1:n, and nu < 1, thenwhereProof. See Problem 3.1.ClThe θn , and γn notation will be used throughout this book.
Whenever wewrite γn there is an implicit assumption that nu < 1, which is true in virtually any circumstance that might arise with IEEE single or double precisionarithmetic.Applying the lemma to (3.2) we obtain(3.3)This is a backward error result and may be interpreted as follows: the computed inner product is the exact one for a perturbed set of data x 1,. . .
, xn ,y1(1 + θn ),y 2(1 + θ' n ),. . . , yn (1 + θ2) (alternatively, we could perturb the xiand leave the yi alone). Each relative perturbation is certainly bounded byγn = nu/(l - nu), so the perturbations are tiny.The result (3.3) applies to one particular order of evaluation.
It is easy tosee that for any order of evaluation we have, using vector notation,(3.4)where |x| denotes the vector with elements |x i | and inequalities between vectors (and, later, matrices) hold componentwise.A forward error bound follows immediately from (3.4):(3.5)If y = x, so that we are forming a sum of squares x Tx, this result shows thathigh relative accuracy is obtained.
However, in general, high relative accuracyis not guaranteed if |xT y| << |x| T |y|.It is easy to see that precisely the same results (3.3)-(3.5) hold when weuse the no-guard-digit rounding error model (2.6). For example, expression(3.1) becomes= x1 y1(1 + δ1)(1 + δ3) + x2 y2(1 + δ2)(1 + δ4), where δ4 hasreplaced a second occurrence of δ3, but this has no effect on the error bounds.It is worth emphasizing that the constants in the bounds above can bereduced by focusing on particular implementations of an inner product.
Forexample, if n = 2m and we computeBASICS70s1 = x(1: m) T y(1: m )s 2 = x(m + 1:n) Ty(m + 1:n)s n = s1 + s 2thenBy accumulating the inner product in twopieces we have almost halved the error bound. This idea can be generalized by breaking the inner product into k pieces, with each mini innerproduct of length n/k being evaluated separately and the results summed.The error bound is now γn / k + k- 1 |x T ||y|, which achieves its minimal value of(or, rather, we should take k to be the nearestinteger toBut it is possible to do even better by using pairwise summation of the products xi yi (this method is described in §4.1).
With pairwisesummation, the error bound becomesSince many of the error analyses in this book are built upon the error analysisof inner products, it follows that the constants in these higher level boundscan also be reduced by using one of these nonstandard inner product implementations. The main significance of this observation is that we should notattach too much significance to the precise values of the constants in errorbounds.Inner products are amenable to being calculated in extended precision.If the working precision involves a t-digit mantissa then the product of twofloating point numbers has a mantissa of 2t - 1 or 2t digits and so can berepresented exactly with a 2t-digit mantissa.
Some computers always form the2t-digit product before rounding to t digits. thus allowing an inner product tobe accumulated at 2t-digit precision at little or no extra cost, prior to a finalrounding.The extended precision computation can be expressed as fl( f l e( x T y )),where fle denotes computations with unit roundoff u e (ue < u). Defining= fl e(xT y), the analysis above holds with u replaced by ue in (3.3) (andwith the subscripts on the θi , reduced by 1 if the multiplications are doneexactly). For the final rounding we haveand so, overall,Hence, as long as nu e(|x|T|y| < u|xT y|, the computed inner product is aboutas good as the rounded exact inner product.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
