Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 46
Текст из файла (страница 46)
LAPACKDriver routines xPOSV (simple) and xPOSVX (expert) use the Cholesky factorization to solve a symmetric (or Hermitian) positive definite system oflinear equations with multiple right-hand sides. (There are corresponding227PROBLEMSroutines for packed storage, in which one triangle of the matrix is stored ina one-dimensional array: PP replaces PO in the names.) The expert driverincorporates iterative refinement, condition estimation, and backward andforward error estimation and has an option to scale the system AX = Bt o (D – 1 AD– 1 )DX = D – 1 B, where D = diagModulo the roundingerrors in computing and applying the scaling, the scaling has no effect onthe accuracy of the solution prior to iterative refinement, in view of Theorem 10.6.
The Cholesky factorization is computed by the routine xPOTRF ,which uses a partitioned algorithm that computes R a block row at a time.The drivers xPTSV and xPTSVX for symmetric positive definite tridiagonal matrices use LDLT factorization. LAPACK does not currently contain a routinefor Cholesky factorization of a positive semidefinite matrix, but there is sucha routine in LINPACK (xCHDC ).Driver routines xSYSV (simple) and xSYSVX (expert) use the block LDLTfactorization (computed by the diagonal pivoting method) with partial pivoting to solve a symmetric indefinite system of linear equations with multiple right-hand sides. For Hermitian matrices the corresponding routines arexHESV (simple) and xHESVX (expert). (Variants of these routines for packedstorage have names in which SP replaces SY and HP replaces HE.) The expertdrivers incorporate iterative refinement, condition estimation, and backwardand forward error estimation.
The factorization is computed by the routinexSYTRF or xHETRF.Problems10.1. Show that ifis symmetric positive definite thenWhat does this statement imply about maxi,j|aij|?10.2. If A is a symmetric positive definite matrix, how would you computex T A- 1 x?10.3. Let y =any order. Show thatwherebe evaluated in floating point arithmetic infor all i, and |θk + 1| < γ k +1 .10.4.
Letbe symmetric positive definite. Show that the reducedsubmatrix B of order n—1 at the end of the first stage of GE is also symmetric228positive definite. Deduce that 0 <that the growth factor pn = 1.CHOLESKY FACTORIZATION= akk and hence10.5. Show that the backward error result (10.6) for the solution of a symmetric positive definite linear system by Cholesky factorization implieswhere ||A||M = maxi,j |aij| (which is not a consistent matrix norm—see §6.2).The significance of this result is that the bound for ||∆ A||M/||A||M contains alinear polynomial in n, rather than the quadratic that appears for the 2-normin (10.7).be positive semidefinite of rank r and suppose it10.6. Let A = cp(A)has the Cholesky factorization (10.11) with Π = I.
Show that Z = [W , –I] Tis a basis for the null space of A, where W =10.7. Prove that (10.13) holds for the Cholesky decomposition with completepivoting.10.8. Give an example of a symmetric matrixfor which the leadingprincipal submatrices Ak satisfy det(A k) > 0, k = 1:n, but A is not positivesemidefinite (recall that det(A k) > 0, k = 1:n, implies that A is positivedefinite).
State a condition on the minors of A that is both necessary andsufficient for positive semidefiniteness.10.9. Suppose the outer product Cholesky factorization algorithm terminatesat the (k+1)st stage (see (10.15)), with a negative pivot in the (k + 1, k + 1)position. Show how to construct a direction of negative curvature for A (avector p such that pTAp < 0).10.10. What is wrong with the following argument? A positive semidefinitematrix is the limit of a positive definite one as the smallest eigenvalue tends tozero.
Theorem 10.3 shows that Cholesky factorization is stable for a positivedefinite matrix, and therefore, by continuity, it must be stable for a positivesemidefinite matrix, implying that Theorem 10.14 is unnecessarily weak (since||W|| 2 can be large).10.11. Consider the diagonal pivoting method applied to a symmetric matrix.
Show that with complete pivoting or partial pivoting any 2 × 2 pivotis indefinite. Hence give a formula for the inertia in terms of the block sizesof the block diagonal factor. Show how to avoid overflow in computing theinverse of a 2 × 2 pivot.10.12. Describe the effect of applying the diagonal pivoting method withpartial pivoting to a 2 × 2 symmetric matrix.10.13.
What factorization is computed if the diagonal pivoting method withpartial pivoting is applied to a symmetric positive definite matrix?PROBLEMS22910.14. (Sorensen and Van Loan; see [315, 1991, §5.3.2]) Suppose the partialpivoting strategy for the diagonal pivoting method is modified by redefining(thus “σ n e w = m a x (σo l d ,|a rr | )”). Show that the same growth factor boundholds as before and that for a positive definite matrix no interchanges aredone and only 1 × 1 pivots are used.10.15. Letwhere 0 << 1, and suppose the diagonal pivoting method is applied toA, yielding a factorization PAP T = LDL T.
Show that with partial pivotingis unbounded aswhereas with complete pivotingisbounded independently of10.16. Letbe nonsymmetric positive definite. Show that the Schur complement S =is also positive definite. In other words, show that GEpreserves positive definiteness.10.17.
A matrix of the formwhereandare symmetric positive definite has beencalled a symmetric quasidefinite matrix by Vanderbei [1047, 1995]. Show that(a) A is nonsingular, (b) for any permutation Π, Π T AΠ has an LU factorization, (c) AS is nonsymmetric positive definite, where S = diag(I, –I). (Thislast property reduces the question of the stability of an LDLT factorization ofA to that of the stability of the LU factorization of a nonsymmetric positivedefinite matrix, for which see §10.5. This reduction has been pointed out andexploited by Gill, Saunders, and Shinnerl [448, 1996].)10.18. (RESEARCH PROBLEM) Is the growth factor bound (2.57)n –1 for thediagonal pivoting method with partial pivoting attainable? If not, how bigcan the growth factor be? Similarly, what is a sharp bound for the completepivoting growth factor?10.19. (RESEARCH PROBLEM) Bound the growth factor for Aasen’s method[1, 1971].Previous Home NextChapter 11Iterative RefinementThe ILLIAC’s memory is sufficient to accommodate a system of 39 equationswhen used with Routine 51.The additional length of Routine 100 restricts to 37the number of equations that it can handle.With 37 equations the operation time of Routine 100 is about4 minutes per iteration.—JAMES N.
SNYDER, On the improvement of the Solutions to a Set ofSimultaneous Linear Equations Using the ILLIAC (1955)In a short mantissa computing environmentthe presence of an iterative improvement routine cansignificantly widen the class of solvable Ax = b problems.— GENE H. GOLUB and CHARLES F. VAN LOAN,Matrix Computations (1989)Most problems involve inexact input data andobtaining a highly accurate solution to animprecise problem may not be justified.— J. J. DONGARRA, J. R. BUNCH, C. B. MOLER, and G. W.
STEWART,LINPACK Users’ Guide (1979)231232I TERATIVE R EFINEMENTIterative refinement is an established technique for improving a computedsolutionto a linear system Ax = b. The process consists of three steps:1. Compute r = b – A2. Solve Ad = r.3. Update y =+ d.(Repeat from step 1 if necessary, withreplaced by y).If there were no rounding errors in the computation of r, d, and y, then y wouldbe the exact solution to the system. The idea behind iterative refinement isthat if r and d are computed accurately enough then some improvement inthe accuracy of the solution will be obtained. The economics of iterativerefinement are favorable for solvers based on a factorization of A, becausethe factorization used to computecan be reused in the second step of therefinement.Traditionally, iterative refinement is used with Gaussian elimination (GE),and r is computed in extended precision before being rounded to working precision. Iterative refinement for GE was used in the 1940s on desk calculators,but the first thorough analysis of the method was given by Wilkinson in 1963[1088, 196 3].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
