Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 54
Текст из файла (страница 54)
An observation that simplifies the analysis is thatwe can consider the algorithm as comprising two stages. The first stage isidentical to GE and forms Mn- 1 M n-2 . . . M1 A = U, where U is upper triangular. The second stage reduces U to diagonal form by elementary operations:The solution x is formed as x = D - 1 Nn . . . N2 y, where y = Mn-1 . . . M2 b.The rounding errors in the first stage are precisely the same as those in GE,so it suffices to consider the second stage. We will assume, for simplicity, thatthere are no row interchanges. As with GE, this is equivalent to assumingthat all row interchanges are done at the start of the algorithm.Define Uk+l = Nk. . .N2 U (so U2 = U) and note that Nk and Uk have theformswhere nk = [-u 1 k /u k k , .
. . . -uk- 1 ,k / u k k ] T. The computed matrices obviouslysatisfy(13.25a)(13.25b)(to be precise,is defined as Nk but withwith xk+1 = Nk . . . N2 y, we haveSimilarly,(13.26)MATRIX INVERSION278Because of the structure of∆ k , and fk , we have the useful property thatWithout loss of generality, we now assume that the final diagonal matrix D isthe identity (i.e., the pivots are all 1); this simplifies the analysis a little andhas a negligible effect on the final bounds. Thus (13.25) and (13.26) yield(13.27)(13.28)Nowand, similarly,But defining ñkwe havewhere X = U-1 + O(u) (by (13.27)). HenceCombining (13.27) and (13.28) we have, for the solution of Ux = y ,which gives the componentwise forward error bound(13.29)27913.4 G AUSS –J ORDAN E LIMINATIONTable 13.6. Gauss-Jordan elimination for Ux = b.n ηU , b16 2.0e-1432 6.4e-1064 1.7e-25.8e-117.6e-66.6e4This is an excellent forward error bound: it says that the error in isno larger than the error we would expect for an approximate solution witha tiny componentwise relative backward error.
In other words, the forwarderror is bounded in the sane way as for substitution. However, we have not,and indeed cannot, show that the method is backward stable. The best wecan do is to derive from (13.27) and (13.28) the result(13.30a)(13.30b)(13.30C)using= U + O(u). These bounds show thathas a normwisebackward error bounded approximately byHence the backward error can be guaranteed to be small only if U is well conditioned.
Thisagrees with the comments of Peters and Wilkinson [828, 1975] that “the residuals corresponding to the Gauss-Jordan solution are often larger than thosecorresponding to back-substitution by a factor of order κ.”A numerical example helps to illustrate the results. We take U to bethe upper triangular matrix with 1s on the diagonal and – 1s everywhereabove the diagonal (U = U(1) from (8.2)). This matrix has condition numberTable 13.6 reports the normwise relative backward error(see (7.2)) for b = U x, wherex = e/3. Clearly, GJE is backward unstable for these matrices—the backwarderrors show a dependence onHowever, the relative distance betweenand the computed solution from substitution (not shown in the table) is lessthanwhich shows that the forward error is bounded byconfirming (13.29).By bringing in the error analysis for the reduction to upper triangularform, we obtain an overall result for GJE.Theorem 13.5.
Suppose GJE successfully computes an approximate solutionto Ax = b, where Ais nonsingular. Then(13.31)(13.32)M ATRIX I NVERSION280whereis the factorization computed by GE.Proof. For the first stage (which is just GE), we have A + ∆A1 =by Theorems 9.3 and8.5.Using (13.30), we obtainor A = b – r, where(13.33)The bounds (13.31) and (13.32) follow easily on using (13.30).Theorem 13.5 shows that the stability of GJE depends not only on the sizeof |L||U| (as in GE), but also on the condition of U. The termisan upper bound forand if this bound is sharp then the residual boundis very similar to that for LU factorization. Note that for partial pivoting wehaveThe bounds in Theorem 13.5 have the pleasing property that they areinvariant under row or column scaling of A, though of course if we are usingpartial pivoting then row scaling can change the pivots and alter the bound.As mentioned earlier, to obtain bounds for matrix inversion we simplytake b to be each of the unit vectors in turn.
For example, the residual boundbecomesFor the special case of symmetric positive definite matrices, an informativenormwise result follows from Theorem 13.5. We make the natural assumptionthat symmetry is exploited in the elimination.Corollary 13.6. Suppose GJE successfully computes an approximate solution to Ax = b, whereis symmetric positive definite. Thenwhereis the factorization computed by symmetric GE.Proof. By Theorem 9.3 we gave A + ∆A =ic and satisfies |∆ A| <Defining D =where ∆A is symmetrwe have, by13.5 T HE D ETERMINANTsymmetry, A + ∆A =281RTR. HenceFurthermore, it is staightforward to show thatn (1- γ n ) - 1 ||A|| 2 . The required bounds follow.Corollary 13.6 shows that GJE is forward stable for symmetric positivedefinite matrices, but it bounds the backward error only by a multiple ofκ 2 (A)1/2.
Numerical experiments show that the backward error is usuallymuch less than κ 2 (A) 1 / 2 u, but (very ill-conditioned) matrices can certainly befound for which the backward error is many order of magnitude larger thanu. Hence GJE is not backward stable even for symmetric positive definitematrices.13.5. The DeterminantIt may be too optimistic to hope that determinants willfade out of the mathematical picture in a generation;their notation alone is a thing of beautyto those who can appreciate that sort of beauty.— E.
T. BELL, Review of “Contributions to the History of Determinants,1900–1920”, by Sir Thomas Muir (1931)Like the matrix inverse, the determinant is a quantity that rarely needs tobe computed. The common fallacy that the determinant is a measure of illconditioning is displayed by the observation that ifis orthogonalthen det(aQ) = andet(Q ) = ±an , which can be made arbitrarily small orlarge despite the fact that aQ is perfectly conditioned. Of course, we couldnormalize the matrix before taking its determinant and define, for example,where D –1 A has rows of unit 2-norm. This functionis called the Hadamardcondition number by Birkhoff [99, 1975], because Hadamard’s determinantalinequality (see Problem 13.11) implies that (A) > 1, with equality if andonly if A has mutually orthogonal rows.
Unless A is already row equilibrated(see Problem 13.13), (A) does not relate to the conditioning of linear systemsin any straightforward way.282M ATRIX I NVERSIONAs good a way as any to compute the determinant of a general matrixis via GEPP. If PA = LU thendet(A) = det(P)-1 det(U) = (–1) r u 11 .
. .unn,(13.34)where r is the number of row interchanges during the elimination. If weuse (13.34) then, barring underflow and overflow, the computed determinant= fl(det(A)) satisfieswhere |θ n| < γn , so we have a tiny relative perturbation of the exact determinant of a perturbed matrix A + ∆A, where ∆A is bounded in Theorem 9.3.In other words, we have almost the right determinant for a slightly perturbedmatrix (assuming the growth factor is small).However, underflow and overflow in computing the determinant are quitelikely, and the determinant itself may not be a representable number. One possibility is to compute log |det(A )log |uii |; as pointed out by Forsytheand Moler [396, 1967, p.
55], however, the computed sum may be inaccuratedue to cancellation. Another approach, used in LINPACK, is to compute andrepresent det(A) in the form y × 10e, where 1 < |y| < 10 and e is an integerexponent.13.5.1. Hyman’s MethodIn §1.16 we analysed the evaluation of the determinant of a Hessenberg matrixH by GE. Another method for computing det(H) is Hyman’s method [594,1957 ], which is superficially similar to GE but algorithmically different. Hyman’s method has an easy derivation in terms of LU factorization that leads toa very short error analysis.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
