Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 8
Текст из файла (страница 8)
For the numericalexample above, (1.6) produces the exact answer. The updating formulae (1.6)are numerically stable, though their error bound is not as small as the onefor the two-pass formula (it is proportional to the condition number KN inProblem 1.7).The problem of computing the sample variance illustrates well how mathematically equivalent formulae can have different numerical stability properties.1.10.
Solving Linear EquationsFor an approximate solution y to a linear system Ax = b ( AIRn × n , bnIR ) the forward error is defined as ||x-y||/||x||, for some appropriate norm.14PRINCIPLES OF FINITE PRECISION COMPUTATIONAnother measure of the quality of y, more or less important depending onthe circumstances, is the size of the residual r = b - Ay . When the linearsystem comes from an interpolation problem, for example, we are probablymore interested in how closely Ay represents b than in the accuracy of y. Theresidual is scale dependent: multiply A and b by a and r is multiplied bya. One way to obtain a scale-independent quantity is to divide by || A|| ||y||,yielding the relative residualThe importance of the relative residual is explained by the following resuit, which was probably first proved by Wilkinson (see the Notes and References).
We use the 2-norm, defined by ||x||2 = (xT x)½ and ||A||2 =Lemma 1.1. With the notation above, and for the 2-norm,Proof. If (A+∆A)y = b then r := b-Ay = ∆Ay, so ||r| | 2 < | |∆A| | 2 | |y | | 2 ,giving(l.7)On the other hand, (A +∆A)y = b for ∆A = ryT /(yT y) and ||∆A||2 =|| r||2 /||y||2, so the bound (1.7) is attainable.Lemma 1.1 says that p(y) measures how much A (but not b) must beperturbed in order for y to be the exact solution to the perturbed system,that is, p(y) equals a normwise relative backward error.
If the data A and bare uncertain and p(y) is no larger than this uncertainty (e.g., p(y) = O(u))then the approximate solution y must be regarded as very satisfactory. Forother problems the backward error may not be as easy to compute as it is fora general linear system, as we will see for the Sylvester equation (§15.2) andthe least squares problem (§19.7).To illustrate these concepts we consider two specific linear equation solvers:Gaussian elimination with partial pivoting (GEPP) and Cramer’s rule.1.10.1.
GEPP Versus Cramer’s RuleCramer’s rule says that the components of the solution to a linear systemAx = b are given by x i = det(A i (b))/det(A), where Ai ( b) denotes A with itsith column replaced by b. These formulae are a prime example of a method1.10 SOLVING LINEAR EQUATIONS15that is mathematically elegant, but useless for solving practical problems.The two flaws in Cramer’s rule are its computational expense and its numerical instability. The computational expense needs little comment, and is,fortunately, explained in most modern linear algebra textbooks (for example,Strang [961, 1993] cautions the student “it would be crazy to solve equationsthat way”).
The numerical instability is less well known, but not surprising.It is present even for n = 2, as a numerical example shows.We formed a 2 × 2 system Ax = b with condition number K2 (A) =||A ||2|| A -1||21013, and solved the system by both Cramer’s rule and GEPPin MATLAB (unit roundoff u1.1 × 10 -16 ). The results were as follows,where r = b - A :Cramer’s rule1.00002.00011.5075 × 10-71.9285 × 10-7GEPP1.00022.0004-4.5689 × 10 -17-2.1931 × 10 -17The scaled residual for GEPP is pleasantly small-of order the unit roundoff.
That for Cramer’s rule is ten orders of magnitude larger, showing that thecomputed solutionfrom Cramer’s rule does not closely satisfy the equations,or, equivalently, does not solve a nearby system. The solutions themselves aresimilar, both being accurate to three significant figures in each component butincorrect in the fourth significant figure. This is the accuracy we would expectfrom GEPP because of the rule of thumb “forward errorbackward error ×condition number”. That Cramer’s rule is as accurate as GEPP in this example, despite its large residual, is perhaps surprising, but it is explained bythe fact that Cramer’s rule is forward stable for n = 2; see Problem 1.9. Forgeneral n, the accuracy and stability of Cramer’s rule depend on the methodused to evaluate the determinants, and satisfactory bounds are not knowneven for the case where the determinants are evaluated by GEPP.The small residual produced by GEPP in this example is typical: erroranalysis shows that GEPP is guaranteed to produce a relative residual oforder u when n = 2 (see §9.2).
To see how remarkable a property this is,consider the rounded version of the exact solution: z = fl(x) = x + ∆x,where ||∆x ||2 <u ||x ||2. The residual of z satisfies ||b- A z||2=||- A ∆ x || 2 <u||A|| 2 ||x||2u||A||2 ||z||2. Thus the computed solution from GEPP has aboutas small a residual as the rounded exact solution, irrespective of its accuracy.Expressed another way, the errors in GEPP are highly correlated so as toproduce a small residual. To emphasize this point, the vector [1.0006,2.0012],which agrees with the exact solution of the above problem to five significantfigures (and therefore is more accurate than the solution produced by GEPP),has a relative residualof order 10-6 .16PRINCIPLES OF FINITE PRECISION COMPUTATIONTable 1.1.
Computed approximationsn10110210310 410 510 610 72.5937432.7048112.7170512.7185972.7219622.5952273.293968= fl((1+1/ n ) n ) to e = 2.71828. . . .1.251.351.233.153.681.235.76×××××××10 -110 -210 -310 -410 -310 -110 -11.11. Accumulation of Rounding ErrorsSince the first electronic computers were developed in the 1940s, commentsalong the following lines have often been made: “The enormous speed ofcurrent machines means that in a typical problem many millions of floatingpoint operations are performed.
This in turn means that rounding errors canpotentially accumulate in a disastrous way.” This sentiment is true, but misleading. Most often, instability is caused not by the accumulation of millionsof rounding errors. but by the insidious growth of just a few rounding errors.As an example, let us approximate e = exp(1) by taking finite n in thedefinition e :=Table 1.1 gives results computed in Fortran 90 in single precision (u6 × 10-8).The approximations are poor, degrading as n approaches the reciprocalof the machine precision. For n a power of 10, l/n has a nonterminatingbinary expansion. When 1+1/n is formed for n a large power of 10, onlya few significant digits from l/n are retained in the sum.
The subsequentexponentiation to the power n, even if done exactly, must produce an inaccurate approximation to e(indeed, doing the exponentiation in double precisiondoes not change any of the numbers shown in Table 1.1). Therefore a singlerounding error is responsible for the poor results in Table 1.1.There is a way to compute (1+1/n) n more accurately, using only singleprecision arithmetic; it is the subject of Problem 1.5.Strassen’s method for fast matrix multiplication provides another example of the unpredictable relation between the number of arithmetic operationsand the error. If we evaluate fl(AB) by Strassen’s method, for n×n matricesA and B, and we look at the error as a function of the recursion thresholdn 0 <n. we find that while the number of operations decreases as n 0 decreasesfrom n to 8, the error typically increases; see §22.2.2.1.12 INSTABILITY WITHOUT CANCELLATION171.12.
Instability Without CancellationIt is tempting to assume that calculations free from subtractive cancellationmust be accurate and stable, especially if they involve only a small numberof operations. The three examples in this section show the fallacy of thisassumption.1.12.1. The Need for PivotingSuppose we wish to compute an LU factorizationClearly, u 11 =and u 22 = 1 - l2 1 u12 =In floating point arithmetic, ifis sufficiently small thenevaluates toAssuming l21 is computed exactly, we then haveThus the computed LU factors fail completely to-reproduce A. Notice thatthere is no subtraction in the formation of L and U. Furthermore, the matrixA is very well conditionedThe problem, of course, is withthe choice ofas the pivot. The partial pivoting strategy would interchangethe two rows of A before factorizing it, resulting in a stable factorization.1.12.2. An Innocuous Calculation?For any x>0 the following computation leaves x unchanged:for i = 1:60endfor i = 1:60x = x2endSince the computation involves no subtractions and all the intermediate numbers lie between 1 and x, we might expect it to return an accurate approximation to x in floating point arithmetic.On the HP 48G calculator, starting with x = 100 the algorithm producesx = 1.0.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
