Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 20
Текст из файла (страница 20)
The effect of using extended3.2 T HE P URPOSEOFR OUNDING E RROR A NALYSIS71precision inner products in an algorithm is typically to enable a factor n tobe removed from the overall error bound.Because extended precision inner product calculations are machine dependent it is difficult or impossible to write portable programs that use them.Most modern numerical codes (for example those in EISPACK, LINPACK,and LAPACK) do not use extended precision inner products. One process inwhich these more accurate products are needed is the traditional formulationof iterative refinement, in which the aim is to improve the accuracy of thecomputed solution to a linear system (see Chapter 11).We have seen that computation of an inner product is a backward stableprocess.
What can be said for an outer product A = xy T, where x, y,IRn ?The analysis is easy. We have âij = xi yi (1 + δ i j ), |δi j | < u, so = xyT + ∆,| ∆| < u|xyT |.(3.6)This is a satisfying result, but the computation is not backward stable. Infact,  = (x + ∆x)(y + ∆y)T does not hold for any ∆x and ∆y (let alone asmall ∆x and ∆y) because A is not in general a rank 1 matrix.This distinction between inner and outer products illustrates a generalprinciple: a numerical process is more likely to be backward stable when thenumber of outputs is small compared with the number of inputs, so that thereis an abundance of data onto which to “throw the backward error”.
An innerproduct has the minimum number of outputs for its 2n scalar inputs, andan outer product has the maximum number (among standard linear algebraoperations).3.2. The Purpose of Rounding Error AnalysisBefore embarking on further error analyses, it is worthwhile to consider whata rounding error analysis is designed to achieve. The purpose is to show theexistence of an a priori bound for some appropriate measure of the effectsof rounding errors on an algorithm. Whether a bound exists is the mostimportant question. Ideally, the bound is small for all choices of problemdata. If not, it should reveal features of the algorithm that characterize anypotential instability, and thereby suggest how the instability can be curedor avoided.
For some unstable algorithms, however, there is no useful errorbound. (For example, no bound is known for the loss of orthogonality due torounding error in the classical Gram-Schmidt method; see §18.7.)The constant terms in an error bound (those depending only on the problem dimensions) are the least important parts of it. As discussed in §2.6, theconstants usually cause the bound to overestimate the actual error by ordersof magnitude.
It is not worth spending much effort to minimize constantsbecause the achievable improvements are usually insignificant.BASICS72It is worth spending effort, though, to put error bounds in a concise, easilyinterpreted form. Part of the secret is in the choice of notation, which wediscuss in §3.4, including the question of what symbols to choose for variables(see the discussion in Higham [554, 1993, §3.5]).If sharp error estimates or bounds are desired they should be computeda posteriori, so that the actual rounding errors that occur can be taken intoaccount. One approach is to use running error analysis, described in the nextsection. Other possibilities are to compute the backward error explicitly, ascan be done for linear equation and least squares problems (see §§7.1, 7.2, and19.7), or to apply iterative refinement to obtain a correct ion that approximatesthe forward error (see Chapter 11).3.3.
Running Error AnalysisThe forward error bound (3.5) is an a priori bound that does not depend onthe actual rounding errors committed. We can derive a sharper, a posteriori bound by reworking the analysis. The inner product evaluation may beexpressed ass0 = 0for i = 1:nsi = si-1 + xi yiendWrite the computed partial sums ashave, using (2.5),Similarly,Hence ei = ei-1 -and letWeorwhich givesSince e0 = 0, we have |e n| < uµ n , whereµ i = µ i-1 +µ 0 = 0.Algorithm 3.2. Given x, yIRn this algorithm computes s = fl(x T y) andTa number µ such that |s - x y| < µ.3.4 N OTATIONFORERROR ANALYSIS73s =0µ =0for i = 1:nz = xi yis= s+ zµ = µ + |s| + |z|endµ = µ * µThis type of computation, where an error bound is computed concurrentlywith the solution, is called running error analysis.
The underlying idea issimple: we use the modified form (2.5) of the standard rounding error modelto write| x op y - fl(x op y)|<u|f l(x op y)|,which gives a bound for the error in x op y that is easily computed, sincefl(x op y) is stored on the computer. Key features of a running error analysisare that few inequalities are involved in the derivation of the bound and thatthe actual computed intermediate quantities are used, enabling advantageto be taken of cancellation and zero operands. A running error bound is,therefore, usually smaller than an a priori one.There are, of course, rounding errors in the computation of the runningerror bound, but their effect is negligible for nu << 1 (we do not need manycorrect significant digits in an error bound).Running error analysis is a somewhat neglected practice nowadays, but itwas widely used by Wilkinson in the early years of computing. It is applicableto almost any numerical algorithm.
Wilkinson explains [1101, 1986]When doing running error analysis on the ACE at no time did Iwrite down these expressions. I merely took an existing program(without any error analysis) and modified it as follows. As eacharithmetic operation was performed I added the absolute value ofthe computed result (or of the dividend) into the accumulatingerror bound.For more on the derivation of running error bounds see Wilkinson [1100, 1985]or [1101, 1986]. A running error analysis for Horner’s method is given in §5.1.3.4.
Notation for Error AnalysisAnother way to write (3.5) is|xT y - fl(xT y)| < nu|x|T|y| + O(u 2 ).(3.7)74BASICSIn general, which form of bound is preferable depends on the context. Theuse of first-order bounds such as (3.7) can simplify the algebra in an analysis.But there can be some doubt as to the size of the constant term concealed bythe big-oh notation. Furthermore, in vector inequalities an O(u 2) term hidesthe structure of the vector it is bounding and so can make interpretation ofthe result difficult; for example, the inequality |x - y| < nu|x| + O(u 2) doesnot imply that y approximates every element of x with the same relative error(indeed the relative error could be infinite when xi = 0, as far as we can tellfrom the bound).In more complicated analyses based on Lemma 3.1 it is necessary to manipulate the 1 + θk and γk terms. The next, lemma provides the necessaryrules.Lemma 3.3.
For any positive integer k let θ k denote a quantity boundedaccording to |θ k| < γk = ku/(1 - ku). The following relations hold:Proof. See Problem 3.4.Concerning the second rule in Lemma 3.3, we certainly havebut if we are given only the expression (1 + θ k )/(1 + θj ) and the bounds forθ k and θj , we cannot do better than to replace it by θ k+2j for j > k.Another style of writing bounds is made possible by the following lemma.L e m m a 3 . 4 . I f |δ| < u for i = 1:n and nu < 0.01, thenwhere |η n | < 1.01nu.3.4 N OTATIONFORERROR ANALYSIS75Proof.
We haveSince 1 + x < ex for x > 0, we have (1 + u) n < exp(nu), and soNote that this lemma is slightly stronger than the corresponding bound wecan obtain from Lemma 3.1: |θ n | < nu/(1 - nu) < nu/0.99 = 1.0101. . . nu.Lemma 3.4 enables us to derive, as an alternative to (3.5),| xT y -f l( x T y ) |<1 . 0 1n u|x | T |y |.(3.8)A convenient device for keeping track of powers of 1 + δ terms was introduced by Stewart [941, 1973, App. 3]. His relative error counter < k > denotesa product(3.9)The counters can be manipulated using the rules<j><k> = <j + k>,At the end of an analysis it is necessary to bound |<k > - 1|, for which anyof the techniques described above can be used.Wilkinson explained in [1100, 1985] that he used a similar notation in hisresearch, writingfor a product of r factors 1 + δi with |δi | < u.
He alsoderived results for specific values of n, before treating the general case-auseful trick of the trade!An alternative notation to fl(·) to denote the rounded value of a numberor the computed value of an expression is [·], suggested by Kahan. Thus, wewould write [a + [b * c]] instead of fl(a + fl(b * c)).BASICS76A completely different notation for error analysis has been proposed byOlver [807, 1978], and subsequently used by him and several other authors.For scalars x and y of the same sign, Olver defines the relative precision rp asfollows:yx; rp(a) means that y = eδ x, | δ| < a.Since eδ = 1 + δ + O(δ2), this definition implies that the relative error in x asan approximation to y (or vice versa) is at most a + O (a 2). But, unlike theusual definition of relative error, the relative precision possesses the propertiesofsymmetry: yx ; rp( a)xx; rp(a ) and zy; rp(β)additivity: yy; rp(a),zx : rp(a + β).Proponents of relative precision claim that the symmetry and additivity properties make it easier to work with than the relative error.Pryce [845, 1981] gives an excellent appraisal of relative precision, withmany examples.
He uses the additional notation 1(δ ) to mean a number θwith θ1; rp(δ). The 1(δ) notation is the analogue for relative precision ofStewart’s <k> counter for relative error. In later papers, Pryce extends therp notation to vector and matrices and shows how it can be used in the erroranalysis of some basic matrix computations [846, 1984], [847, 1985].Relative precision has not achieved wide use. The important thing for anerror analyst is to settle on a comfortable notation that does not hinder thethinking process. It does not really matter which of the notations describedin this section is used, as long as the final result is informative and expressedin a readable form.3.5. Matrix MultiplicationGiven error analysis for inner products it is straightforward to analyse matrixvector and matrix-matrix products.
Let AIR m × n, xIR n and y = Ax.The vector y can be formed as m inner products, yi == 1:m, whereis the ith row of A. From (3.4) we haveThis gives the backward error result(3.10)which implies the forward error bound(3.11)3.5 M ATRIX M ULTIPLICATION77Normwise bounds readily follow (see Chapter 6 for norm definitions): forexample,This inner product formation of y can be expressed algorithmically as% Sdot or inner product form.y(1:m) = 0for i = 1:mfor j = 1:ne n d y ( i ) =y ( i ) +a ( i , j) x ( j )endThe two loops can be interchanged to give% Saxpy form.y(1:m) = 0for j = 1:nfor i = 1:my (i) = y(i) +a( i,j ) x ( j )endendThe terms “sdot” and “saxpy” come from the BLAS (see §D.1).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
