Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 10
Текст из файла (страница 10)
Defining υ =OFF INITE P RECISION C OMPUTATION- 1, we haveFor small x ( y 1 )From (1.9) it follows thatapproximates f with relative error at most about3.5u.The details of the analysis obscure the crucial property that ensures itssuccess. For small x, neither - 1 nor log agrees with its exact arithmeticcounterpart to high accuracy.
But ( - 1)/log is an extremely good approximation to (y-l)/logy near y=1, because the function g(y)=(y-1)/log yvaries so slowly there (g has a removable singularity at 1 and g'(1) = 1). Inother words, the errors in- 1 and logalmost completely cancel.1.14.2. QR FactorizationAny matrix AIR m×n, m > n, has a QR factorization A = QR, where Qm×IRhas orthonormal columns and RIR n×x is upper trapezoidal (rij = 0for i > j). One way of computing the QR factorization is to premultiply A bya sequence of Givens rotations-orthogonal matrices G that differ from theidentity matrix only in a 2×2 principal submatrix, which has the formWith A1 := A, a sequence of matrices Ak satisfying Ak = GkAk-1 is generated. Each A k, has one more zero than the last, so Ap = R for p =n(m - (n + 1)/2).
To be specific, we will assume that the zeros are introduced in the order (n,1). (n - 1,1), . . . , (2,1); (n,2), . . . , (3,2); and so on.For a particular 10×6 matrix A, Figure 1.5 plots the relative errors||Ak - Ak||2/||A||2, where Ak, denotes the matrix computed in single precisionarithmetic (u6 × 10-8). We see that many of the intermediate matrices arevery inaccurate, but the final computedhas an acceptably small relativeerror, of order u. Clearly, there is heavy cancellation of errors on the lastfew stages of the computation. This matrix AIR10×6 was specially chosen,1.14 C ANCELLATIONOFR OUNDING E RRORS25Figure 1.5.
Relative errors ||A k- Â k||2/||A ||2 for Givens QR factorization. Thedotted line is the unit roundoff level.following a suggestion of Wilkinson [1100,||A||21 and A 10 has the form1985],as a full matrix such thatBecause y is at the roundoff level, the computed is the result of severe subtractive cancellation and so is dominated by rounding errors.
Consequently,the computed Givens rotations,...,whose purpose is to zero thevectorand which are determined by ratios involving the elements of , bearlittle relation to their exact counterparts, causing Âk to differ greatly fromAk for k = 11,12,. . . .To shed further light on this behaviour, we note that the Givens QR factorization is perfectly backward stable; that is, the computed R is the exactR factor of A + ∆A, where ||∆A||2 <cu||A||2, with c a modest constant depending on the dimensions (Theorem 18.9). By invoking a perturbation resultfor the QR factorization (namely (18.27)) we conclude thatisbounded by a multiple of K2 (A)u.
Our example is constructed so that K2 (A) issmall ( 24), so we know a priori that the graph in Figure 1.5 must eventuallydip down to the unit roundoff level.P RINCIPLES26OFF INITE P RECISION C OMPUTATIONWe also note thatis of order u in this example, as again we canshow it must be from perturbation theory.
Since Q is a product of Givensrotations, this means that even though some of the intermediate Givens rotations are very inaccurate, their product is highly accurate, so in the formationof Q, too, there is extensive cancellation of rounding errors.1.15. Rounding Errors Can Be BeneficialAn old method for computing the largest eigenvalue (in absolute value) ofa matrix A and the corresponding eigenvector is the power method, whichconsists of repeatedly multiplying a given starting vector by A. With scalingto avoid underflow and. overflow, the process in its simplest form is% Choose a starting vector x.while not convergedx := AxendThe theory says that if A has a unique eigenvalue of largest modulus and xis not deficient in the direction of the corresponding eigenvector υ, then thepower method converges to a multiple of υ (at a linear rate).Consider the matrixwhich has eigenvalues 0, 0.4394 and 1.161 (correct to the digits shown) and aneigenvector [l, 1, 1]T corresponding to the eigenvalue zero.
If we take [1,1,1]Tas the starting vector for the power method then, in principle, the zero vectoris produced in one step, and we obtain no indication of the desired dominanteigenvalue-eigenvector pair. However, when we carry out the computation inMATLAB, the first step produces a vector with elements of order 10-16 andwe obtain after 38 iterations a good approximation to the dominant eigenpair. The explanation is that the matrix A cannot be stored exactly in binary floating point arithmetic.
The computer actually works with A + ∆Afor a tiny perturbation ∆A, and the dominant eigenvalue and eigenvector ofA + ∆A are very good approximations to those of A. The starting vector[1,1,1]T contains a nonzero (though tiny) component of the dominant eigenv e c t o r o f A + ∆A. This component grows rapidly under multiplication byA + ∆A, helped by rounding errors in the multiplication, until convergenceto the dominant eigenvector is obtained.1.16 S TABILITYOF ANA LGORITHM D EPENDSON THEP ROBLEM27Perhaps an even more striking example of beneficial effects of roundingerrors is in inverse iteration, which is just the power method applied to theshifted and inverted matrix (A - µI )-1.
The shift µ is usually an approximateeigenvalue. The closer µ is to an eigenvalue, the more nearly singular A - µIis, and hence the larger the error in computing y = (A - µI)-1 x (which is doneby solving (A-µ I)y = x). However, it can be shown that the error lies almostentirely in the direction of the required eigenvector, and so is harmless; see, forexample, Parlett [820, 1980, §4.31] or Golub and Van Loan [470, 1989, §7.6.1].1.16. Stability of an Algorithm Depends on the ProblemAn algorithm can be stable as a means for solving one problem but unstable when applied to another problem. One example is the modified GramSchmidt method, which is stable when used to solve the least squares problembut can give poor results when used to compute an orthonormal basis of amatrix (see §§518.7 and 19.3).A lesser known and much simpler example is Gaussian elimination (GE)without pivoting for computing the determinant of an upper Hessenberg matrix.
A square matrix A is upper Hessenberg if a ij = 0 for i > j + 1. GEtransforms A to upper triangular form by n - 1 row eliminations. one for eachof the boxed entries in this 4 × 4 illustration:The determinant of A is given by the product of the diagonal elements of U.It is easy to show that this is a stable way to evaluate det(A), even thougharbitrarily large multipliers may arise during the elimination.
Note. first, that,if A( k ) denotes the matrix at the start of the kth stage (A(1) = A), thenbecause the kth row of A( k -1) is the same as the kth row of A. In floatingpoint arithmetic the model (1.1) shows that the computedsatisfy28PRINCIPLESOFFINITE P RECISION COMPUTATIONTable 1.3. Results from GE without pivoting on an upper Hessenberg matrix.ExactComputed Relative error< u, i = 1:3. This equation says that the computed diagonalwhereelementsare the exact diagonal elements corresponding not to A, but toa matrix obtained from A by changing the diagonal elements to akk ( 1 +)In otherand the subdiagonal elements towords, the computedare exact for a matrix differing negligibly from A.The computed determinant d, which is given byis therefore a tiny relative perturbation of the determinant of a matrix differingnegligibly from A.
so this method for evaluating det(A) is numerically stable(in the mixed forward backward error sense of (1.2)).However, if we use GE without pivoting to solve an upper Hessenberglinear system then large multipliers can cause the solution process to be unstable. If we try to extend the analysis above we find that the computed LUfactors (as opposed to just the diagonal of U) do not, as a whole, necessarilycorrespond to a small perturbation of A.A numerical example illustrates these ideas. LetWe took a = 10-7 and b = Ae(e = [1,1,1,1]T) and used GE without pivotingin single precision arithmetic (u6 × 10-8) to solve Ax = b and computedet(A).
The computed and exact answers are shown to five significant figuresin Table 1.3. Not surprisingly, the computed determinant is very accurate.But the computed solution to Ax = b has no correct figures in its first component. This reflects instability of the algorithm rather than ill conditioningof the problem because the condition number= 16. The source of theinstability is the large first multiplier, a 2 1 /a11 = 107.1.17 ROUNDING ERRORS ARE NOT RANDOM29Figure 1.6. Values of rational function r(x) computed by Horner’s rule (marked as"×" ), for x = 1.606 + (k - 1)2-52; solid line is the “exact” r(x).1.17. Rounding Errors Are Not RandomRounding errors, and their accumulated effect on a computation, are notrandom.
This fact underlies the success of many computations, including someof those described earlier in this chapter. The validity of statistical analysis ofrounding errors is discussed in §2.6. Here we simply give a revealing numericalexample (due to W. Kahan).Define the rational functionwhich is expressed in a form corresponding to evaluation of the quartic polynomials in the numerator and denominator by Horner’s rule. We evaluatedr(x) by Horner’s rule in double precision arithmetic for 361 consecutive floating point numbers starting with a = 1.606, namely x = a + (k - 1)2-52,k = 1:361; the function r(x) is virtually constant on this interval.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
