Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 55
Текст из файла (страница 55)
Letbean unreduced upper Hessenbergmatrix (h i+ 1 , i0 for all i) and writeThe matrix H 1 is H with the first row cyclically permuted to the bottom, sodet(H 1 ) = (–1) n -1 det(H). Since T is a nonsingular upper triangular matrix,we have the LU factorization(13.35)from which we obtain det(H 1) = det(T )(η – hTT– 1 y ). Thereforedet(H) = (–1) n -1det(T)(η – hTT - 1y).(13.36)13.6 N OTESANDR EFERENCES283Hyman's method consists of evaluating (13.36) in the natural way: by solvingthe triangular system Tx = y by back substitution, then forming η – hTx andits product with det(T ) = h 2 1 h 32 . .
. hn,n- 1 .For the error analysis we assume that no underflows or overflows occur.From the backward error analysis for substitution (Theorem 8.5) and for aninner product ((3.4)) we have immediately, on defining µ := η – hTT– 1 y ,whereSince fl(det(T)) = det(T)(1 + δ1). . .(1 + δn -1), |δ i | < u, the computed determinant satisfieswhere |θ n| < γn . This is not a backward error result, because only one ofthe two Ts in this expression is perturbed.
However, we can write det( T ) =det(T + ∆T)(1 + θ( n -1)2), so thatWe conclude that the computed determinant is the exact determinant of aperturbed Hessenberg matrix in which each element has undergone a relativeperturbation not exceeding γn2 - n + 1n 2 u , which is a very satisfactory result.In fact, the constant γn 2 - n+1 may be reduced to γ2 n -1 by a slightly modifiedanalysis; see Problem 13.14.13.6.
Notes and ReferencesClassic references on matrix inversion are Wilkinson’s paper Error Analysisof Direct Methods of Matrix Inversion [1085, 1961] and his book RoundingErrors in Algebraic Processes [1088, 1963]. In discussing Method 1 of §13.2.1,Wilkinson says “The residual of X as a left-hand inverse may be larger thanthe residual as a right-hand inverse by a factor as great as ||L|| ||L - 1 || . .
. Weare asserting that the computed X is almost invariably of such a nature thatXL – I is equally small” [1088, 1963, p. 107]. Our experience concurs withthe latter statement. Triangular matrix inversion provides a good example ofthe value of rounding error analysis: it helps us identify potential instabilities,even if they are rarely manifested in practice, and it shows us what kind ofstability is guaranteed to hold.In [1085, 1961] Wilkinson identifies various circumstances in which triangular matrix inverses are obtained to much higher accuracy than the bounds284M ATRIX I NVERSIONof this chapter suggest.
The results of §8.2 provide some insight. For example,if T is a triangular M-matrix then, as noted after Corollary 8.10, its inverseis computed to high relative accuracy, no matter how ill conditioned L maybe.Sections 13.2 and 13.3 are based closely on Du Croz and Higham [322,1992 ].Method D in §13.3.4 is used by the Hewlett-Packard HP-15C calculator,for which the method’s lack of need for extra storage is an important property [523, 1982].Higham [560, 1995] gives error analysis of a divide-and-conquer methodfor inverting a triangular matrix that has some attractions for parallel computation.GJE is an old method. It was discovered independently by the geodesistWilhelm Jordan (1842-1899) (not the mathematician Camille Jordan (18381922)) and B.-I. Clasen [12, 1987].An Algol routine for inverting positive definite matrices by GJE was published in the Handbook for Automatic Computation by Bauer and Reinsch [83,1971].
As a means of solving a single linear system, GJE is 50% more expensivethan GE when cost is measured in flops; the reason is that GJE takes O ( n 3 )flops to solve the upper triangular system that GE solves in n 2 flops. However, GJE has attracted interest for vector computing because it maintainsfull vector lengths throughout the algorithm. Hoffmann [577, 1987] found thatit was faster than GE on the CDC Cyber 205, a now-defunct machine with arelatively large vector startup overhead.Turing [1027, 1948] gives a simplified analysis of the propagation of errorsin GJE, obtaining a forward error bound proportional to κ(A). Bauer [80,1966] does a componentwise forward error analysis of GJE for matrix inversion and obtains a relative error bound proportional tofor symmetricpositive definite A.
Bauer’s paper is in German and his analysis is not easy tofollow. A summary of Bauer’s analysis (in English) is given by Meinguet [746,19 6 9 ].The first thorough analysis of the stability of GJE was by Peters andWilkinson [828, 1975]. Their paper is a paragon of rounding error analysis.Peters and Wilkinson observe the connection with GE, then devote their attention to the “second stage” of GJE, in which Ux = y is solved. They showthat each component of x is obtained by solving a lower triangular system, andthey deduce that, for each i, (U + ∆U i )x( i ) = y + ∆y i , where |∆Ui | < γn|U|and |∆yi | < γn |y|, and where the ith component of x( i ) is the ith componentofThey then show thathas relative error bounded bybut thatdoes not necessarily have a small backward error. The more directapproach used in our analysis is similar to that of Dekker and Hoffmann [276,1989], who give a normwise analysis of a variant of GJE that uses row pivoting (elimination across rows) and column interchanges.
Our componentwisePROBLEMS285bounds (13.29)–(13.32) are new.Goodnight [473, 1979] describes the use of GJE in statistics for solvingleast squares problems.Error analysis of Hyman’s method is given by Wilkinson [1084, 1960],[1088, 196 3, pp.
147-154], [1089, 196 5, pp. 426-431]. Although it datesfrom the 1950s, Hyman’s method is not obsolete: it has found use in methods designed for high-performance computers; see Ward [1065, 1976], Li andZeng [701, 1992], and Dubrulle and Golub [324, 1994].13.6.1. LAPACKRoutine xGETRI computes the inverse of a general square matrix by MethodB using an LU factorization computed by xGETRF . The corresponding routinefor a symmetric positive definite matrix is xPOTRI , which uses Method D,with a Cholesky factorization computed by xPOTRF .
Inversion of a symmetricindefinite matrix is done by xSYTRI . Triangular matrix inversion is done byxTRTRI , which uses Method 2C. None of the LAPACK routines compute thedeterminant, but it is easy for the user to evaluate it after computing an LUfactorization.Problems13.1. Reflect on this cautionary tale told by Acton [4,1970,p. 246].“It was 1949 in Southern California. Our computer was a very new CPC(model 1, number 1) —a 1-second-per-arithmetic-operation clunker that washolding the computational fort while an early electronic monster was beingcoaxed to life in an adjacent room, From a nearby aircraft company therearrived one day a 16 × 16 matrix of 10-digit numbers whose inverse was desired.
. . We labored for two days and, after the usual number of glitches thataccompany any strange procedure involving repeated handling of intermediatedecks of data cards, we were possessed of an inverse matrix. During thechecking operations . . . it was noted that, to eight significant figures, theinverse was the transpose of the original matrix! A hurried visit to the aircraftcompany to explore the source of the matrix revealed that each element hadbeen laboriously hand computed from some rather simple combinations ofsines and cosines of a common angle. It took about 10 minutes to prove thatthe matrix was, indeed, orthogonal!”13.2. Rework the analysis of the methods of §13.2.2 using the assumptions(12.3) and (12.4), thus catering for possible use of fast matrix multiplicationtechniques.MATRIX INVERSION28613.3.
Show that for any nonsingular matrix A,This inequality shows that the left and right residuals of X as an approximation to A–1 can differ greatly only if A is ill conditioned.13.4. (Mendelssohn [748, 1956]) Find parametrized 2 × 2 matrices A and Xsuch that the ratio ||AX – I||/||XA – I|| can be made arbitrarily large.13.5. Let X and Y be approximate inverses ofthat satisfyandShow thatDerive forward error bounds forandInterpret all these bounds.13.6. What is the relation between the matrix on the front cover of the LAPACK Users’ Guide [17, 1995] and that on the back cover? Answer the samequestion for the LINPACK Users’ Guide [307, 1979].13.7. Show that for any matrix having a row or column of 1s, the elementsof the inverse sum to 1.13.8. Let X = A + iBbe nonsingular.
Show that X-1expressed in terms of the inverse of the real matrix of order 2 n,canbeShow that if X is Hermitian positive definite then Y is symmetric positivedefinite. Compare the economics of real versus complex inversion.13.9. For a given nonsingularand XA -1, it is interesting to askwhat is the smallest ε such that X + ∆X = (A + ∆A)-1 with ||∆ X|| <andWe require (A + ∆A)(X + ∆X) = I, orA ∆X + ∆AX + ∆A∆X = I – AX.It is reasonable to drop the second-order term, leaving a generalized Sylvesterequation that can be solved using singular value decompositions of A and X(cf.
§15.2). Investigate this approach computationally for a variety of A andmethods of inversion.PROBLEMS28713.10. For a nonsingularand given integers i and j, under whatconditions is det(A) independent of aij ? What does this imply about thesuitability of det(A ) for measuring conditioning?13.11. Prove Hadamard’s inequality forwhere ak = A(:, k). When is there equality? (Hint: use the QR factorization.)13.12. (a) Show that if AT = QR is a QR factorization then the Hadamardcondition numberwhere pi = ||R(:, i)||2.
(b) Evaluatefor A = U(1) (see (8.2)) and for the Pei matrix A = (a – 1)I + eeT.13.13. (Guggenheimer, Edelman, and Johnson [486,a nonsingular1995])(a) Prove that for(Hint: apply the arithmetic-geometric mean inequality to the n numberswhere the σ i are the singular values of A.) (b)Deduce that if A has rows of unit 2-norm thenwhereis the Hadamard condition number.13.14. Show that Hyman’s method for computing det(H), whereis an unreduced upper Hessenberg matrix, computes the exact determinantof H + ∆H where |∆ H| < γ2 n - 1 |H|, barring underflow and overflow. Whatis the effect on the error bound of a diagonal similarity transformation HD - 1 HD, where D = diag(di ), di0?13.15. What is the condition number of the determinant?13.16. (RESEARCH PROBLEM) Obtain backward and forward error boundsfor GJE applied to a diagonally dominant matrix.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
