Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 21
Текст из файла (страница 21)
Sdot standsfor (single precision) dot product, and saxpy for (single precision) a timesx plus y. The question of interest is whether (3.10) and (3.11) hold for thesaxpy form. They do: the saxpy algorithm still forms the inner productsbut instead of forming each one in turn it evaluates them all “in parallel”, aterm at a time. The key observation is that exactly the same operations areperformed, and hence exactly the same rounding errors are committed-theonly difference is the order in which the rounding errors are created.This “rounding error equivalence” of algorithms that are mathematicallyidentical but algorithmically different is one that occurs frequently in matrixcomputations. The equivalence is not always as easy to see as it is for matrixvector products.Now consider matrix multiplication: C = AB, where AIRm×n andn×pBIR .
Matrix multiplication is a triple loop procedure, with six possibleloop orderings, one of which isC(1:m,1:p) = 0for i = 1:mfor j = 1:pfor k = 1:nC ( i , j ) =C ( i , j ) +A ( i , k ) B ( k , j )78BASICSendendendAs for the matrix -vector product. all six versions commit the same roundingerrors, so it suffices to consider any one of them. The “jik” and “jki” orderingsboth compute C a column at a time: c j = Ab j, where cj = C(:,j) andb j = B(:,j). From (3.10),Hence the jth computed column of C has a small backward error: it is theexact jth column for slightly perturbed data. The same cannot be said foras a whole (see Problem 3.5 for a possibly large backward error bound).However, we have the forward error bound(3.12)and the corresponding normwise bounds includeThe bound (3.12) falls short of the ideal boundwhich saysthat each component of C is computed with high relative accuracy.
Nevertheless (3.12) is the best bound we can expect, because it reflects the sensitivity of the product to componentwise relative perturbations in the data:for any i and j we can find a perturbation ∆A with |∆A| < u|A| such that|(A+∆A)B-AB|ij =u(|A||B|)ij (similarly for perturbations in B).3.6. Complex ArithmeticTo carry out error analysis of algorithms in complex arithmetic we need amodel for the basic arithmetic operations. Since complex arithmetic must beimplemented using real arithmetic, the complex model is a consequence of thecorresponding real one.
We will assume that for complex numbers x = a + iband y = c + id we computex ± y = a + c ± i(b + d),xy = ac - bd + i(ad + bc),(3.13a)(3.13b)(3.13c)3.6 COMPLEX ARITHMETIC79Lemma 3.5. For x, ythe basic arithmetic operations computed accordingto (3.13) under the standard model (2.4) satisfyProof. Throughout the proof, δi denotes a number bounded by |δi | < u.Addition/subtraction:soas required.Multiplication:whereas required.Division:Then80B ASICSwhere, using Lemma 3.3,Using the analogous formula for the error in fl(Imx / y ) ,which completes the proof.It is worth stressing that δ in Lemma 3.5 is a complex number, so wecannot conclude from the lemma that the real and imaginary parts of fl( x op y)are obtained to high relative accuracy---only that they are obtained to highaccuracy relative to |x op y|.As explained in §25.8, the formula (3.13c) is not recommended for practicaluse since it is susceptible to overflow.
For the alternative formula (25.1), whichavoids overflow, similar analysis to that in the proof of Lemma 3.5 shows thatBounds for the rounding errors in the basic complex arithmetic operationsare rarely given in the literature. Indeed, virtually all published error analysesin matrix computations are for real arithmetic. However, because the boundsof Lemma 3.5 are of the same form as for the standard model (2.4) for realarithmetic, most results for real arithmetic (including virtually all those inthis book) are valid for complex arithmetic, provided that the constants areincreased appropriately.3.7. MiscellanyIn this section we give some miscellaneous results that will be needed in laterchapters.
The first two results provide convenient ways to bound the effectof perturbations in a matrix product. The first result uses norms and thesecond. components.Lemma 3.6. If X j + ∆ Xconsistent norm, thenjIRn×n satisfies ||∆Xj|| < dj||Xj|| for all j for a813.7 MISCELLANYProof. The proof is a straightforward induction, which we leave as anexercise (Problem 3.10).A componentwise result is entirely analogous.Lemma 3.7. If Xj + ∆XjIRn×n satisfies |∆X j| < δj |X j | for all j, thenThe final result describes the computation of the “rank 1 update” y =(I - abT )x, which is an operation arising in various algorithms, including theGram-Schmidt method and Householder QR factorization.Lemma 3.8. Let a, b, x IR n and let y = (I - ab T )x be computed asfl(x- a( b T x )).
Then= y + ∆y, whereso thatProof. Consider first the computation of w = a(bTx). We havewhereFinally,= fl(x -) satisfiesandHence= y + ∆y, where=82B ASICS3.8. Error Analysis DemystifiedThe principles underlying an error analysis can easily be obscured by thedetails. It is therefore instructive to examine the basic mechanism of forwardand backward error analyses. We outline a general framework that revealsthe essential simplicity.Consider the problem of computing z = f(a), where f : IRnIRm .Any algorithm for computing z can be expressed as follows.
Let x 1 = a andx k+l = g k (x k), k = 1:p, whereThe kth stage of the algorithm represents a single floating point operationand xk contains the original data together with all the intermediate quantitiescomputed so far. Finally, z =where is comprised of a subset of thecolumns of the identity matrix (so that each zi is a component of xp+1). Infloating point arithmetic we havewhere ∆xk+i represents the rounding errors on the kth stage and should beeasy to bound. We assume that the functions g k are continuously differentiableand denote the Jacobian of gk at a by Jk. Then, to first order,The pattern is clear: for the finalwe haveIn a forward error analysis we bound f(a) - , which requires bounds for(products of) the Jacobians Jk.
In a backward error analysis we write, againto first order,3.9 O THER A PPROACHES83where Jf is the Jacobian of f. So we need to solve, for ∆a,In most matrix problems there are fewer outputs than inputs (m < n), sothis is an underdetermined system. For a normwise backward error analysiswe want a solution of minimum norm. For a componentwise backward erroranalysis, in which we may want (for example) to minimize subject to | ∆ a| <we can writeand then we want the solution c of minimal m-norm.The conclusions are that forward error analysis corresponds to boundingderivatives and that backward error analysis corresponds to solving a largeunderdetermined linear system for a solution of minimal norm. In principal,therefore, error analysis is straightforward! Complicating factors in practiceare that the Jacobians Jk may be difficult to obtain explicitly, that an errorbound has to be expressed in a form that can easily be interpreted, and thatwe may want to keep track of higher-order terms.3.9.
Other ApproachesIn this book we do not describe all possible approaches to error analysis. Someothers are mentioned in this section.Linearized rounding error bounds can be developed by applying equationsthat describe, to first order, the propagation of absolute or relative errors inthe elementary operations +,-,*,/.
The basics of this approach are given inmany textbooks (see, for example, Dahlquist and Björck [262, 1974, §2.2] orStoer and Bulirsch [955, 1980, §1.3]), but for a thorough treatment see Stummel [963, 198 0], [964, 198 1]. Ziv [1133, 1995] shows that linearized boundscan be turned into true bounds by increasing them by a factor that dependson the algorithm.Rounding error analysis can be phrased in terms of graphs. This appearsto have been first suggested by McCracken and Dorn [743, 196 4], who use“process graphs” to represent a numerical computation and thereby to analysethe propagation of rounding errors.
Subsequent more detailed treatmentsinclude those of Bauer [82, 1974], Miller [758, 1976], and Yalamov [1117, 1994].The work on graphs falls under the heading of automatic error analysis (formore on which see Chapter 24) because processing of the large graphs requiredto represent practical computations is impractical by hand. Linnainmaa [705,1976] shows how to compute the Taylor series expansion of the forward error84B ASICSin an algorithm in terms of the individual rounding errors, and he presents agraph framework for the computation.Some authors have taken a computational complexity approach to erroranalysis, by aiming to minimize the number of rounding error terms in aforward error bound, perhaps by rearranging a computation. Because thisapproach ignores the possibility of cancellation of rounding errors, the resultsneed to be interpreted with care.
See Aggarwal and Burgmeier [6, 1979] andTsao [1024, 1983].3.10. Notes and ReferencesThe use of Lemma 3.1 for rounding error analysis appears to originate withthe original 1972 German edition of a book by Stoer and Bulirsch [955, 1980].The lemma is also used, with p i = 1, by Shampine and Allen [911, 1973,p. 18].Lemma 3.4 is given by Forsythe and Moler [396, 1967, p. 92]. Wilkinsonmade frequent use of a slightly different version of Lemma 3.4 in which theassumption is nu < 0.1 and the bound for |ηn | is 1.06nu (see, e.g., [1089,1965, p.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
