Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 45
Текст из файла (страница 45)
The conclusions are that the computedfactors satisfyand the computed solution to a linear system Ax = b satisfieswhere p1 and p 2 are linear polynomials. For the partial pivoting strategy,Higham shows that if linear systems involving 2 × 2 pivots are solved by GEPPor by use of the explicit inverse, then the computed solutions do indeed havea small componentwise relative backward error, and that, moreover,|| |L||D||LT| ||M < 36n p n||A||M ,where ||A||M = max i,j |aij|.
Thus the diagonal pivoting method with partialpivoting is stable if the growth factor is small.10.5. Nonsymmetric Positive Definite MatricesThe notion of positive definiteness can be extended to nonsymmetric matrices.A nonsymmetric matrixis positive definite if xTAx > 0 for allx0. This is equivalent to the condition that the symmetric part AS ofA is positive definite, where A A S + A K w i t h A S = ( A + A T )/2 andA K = ( A – A T )/2. A positive definite matrix clearly has nonsingular leadingprincipal submatrices, and so has an LU factorization, A = LU. It can evenbe shown that pivots uii are positive.
However, there is no guarantee that thefactorization is stable without pivoting, as the exampleshows. Thestandard error analysis for LU factorization applies (Theorems 9.3 and 9.4),and so the question is whethercan be suitably bounded. Golub andVan Loan [469, 1979] show that, for the exact LU factors,(10.29)Letwhich is just κ 2 (A) when A is symmetric.
Mathias [731, 1992] shows that || |L||U| ||F (involving now the computed LU factors) is at most a factor 1 + 30un3 / 2 X(A) times larger than the224C HOLESKY FACTORIZATIONupper bound in (10.29), and that the LU factorization (without pivoting)succeeds if 24n3 / 2 X (A )u < 1.These results show that it is safe not to pivot provided that the symmetricpart of A is not too ill conditioned relative to the norm of the skew-symmetricpart.
If A is symmetric (AK = 0) then we recover the results for symmetricpositive definite matrices.10.6. Notes and ReferencesAndré-Louis Cholesky (1875-1918) was a French military officer involved ingeodesy and surveying in Crete and North Africa. In some books his nameis misspelled “Choleski”. Details of Cholesky’s life-and a discussion aboutthe pronunciation of his name!-can be found in the electronic mail magazineNA-Digest, volume 90, 1990, issues 7, 8, 10–12, and 24; see, in particular,the biography [22, 1922].
Cholesky’s work was published posthumously on hisbehalf by Benoit [91, 1924].The properties of the Cholesky factorization are intimately associated withthe properties of the Schur complement, as is apparent from some of the proofsin this chapter. The same is true for GE in general. An excellent survey of theSchur complement, containing historical comments, theory, and applications,is given by Cottle [248, 1974].For results on the Cholesky factorization in Hilbert space see Power [841,1986].A book by George and Liu [438, 1981] is devoted to the many practical issues in the implementation of Cholesky factorization for the solution of sparsesymmetric positive definite systems.There is no floating point error analysis of Cholesky factorization in Wilkinson’s books, but he gives a detailed analysis in [1092, 1968], showing that= A+E, with ||E||2 < 2.5n3 / 2 u||A|| 2.
It is unfortunate that this paperis in a rather inaccessible proceedings, because it is a model of how to phraseand interpret an error analysis. Meinguet [747, 1983] and Sun [973, 1992] givecomponentwise backward error bounds similar to those in Theorems 10.3 and10.4. Kielbasinski [657, 198 7] reworks Wilkinson’s analysis to improve theconstant.The fact that κ2 (H) can replace the potentially much larger κ 2 (A) inthe forward error bound for the Cholesky method was stated informally andwithout proof by Wilkinson [1092, 1968, p. 638]. Demmel [283, 1989] madethis observation precise and explored its implications; Theorems 10.5, 10.6,and 10.7 are taken from [283, 1989].The bounds in Theorem 10.8 are from Sun [971, 1991], [972, 1992].
Similarbounds are given by Stewart [944, 1977], [951, 1993], Barrlund [71, 1991],and Sun [973, 1992]. A perturbation bound that can be much smaller than10.6 N OTESANDR EFERENCES225the normwise one in Theorem 10.8 is derived and explored by Chang andPaige [198, 1995]. Perturbation results of a different flavour, including onefor structured perturbations of the form of ∆A in Theorem 10.5, are given byDrmac, Omladc, and Veselic [321, 1994].The perturbation and error analysis of §10.3 for semidefinite matrices isfrom Higham [540, 1990], where in a perturbation result for the QR factorization with column pivoting is also given.
For an application in optimizationthat makes use of Cholesky factorization with complete pivoting and the analysis of §10.3.1 see Forsgren, Gill, and Murray [384, 1995].Fletcher and Powell [382, 1974] describe several algorithms for updatingan LDLT factorization of a symmetric positive definite A when A is modifiedby a rank-1 matrix. They give detailed componentwise error analysis for someof the methods.An excellent way to test whether a given symmetric matrix A is positive(semi) definite is to attempt to compute a Cholesky factorization. This testis less expensive than computing the eigenvalues and is numerically stable.Indeed, if the answer “yes” is obtained, it is the right answer for a nearbymatrix, whereas if the answer is “no” then A must be close to an indefinitematrix.
See Higham [535, 1988] for an application of this definiteness test.An algorithm for testing the definiteness of a Toeplitz matrix is developed byCybenko and Van Loan [260, 1986], as part of a more complicated algorithm.According to Kerr [654, 1990], misconceptions of what is a sufficient conditionfor a matrix to be positive (semi) definite are rife in the engineering literature(for example, that it suffices to check the definiteness of all 2 × 2 submatrices).See also Problem 10.8. For some results on definiteness tests for Toeplitzmatrices, see Makhoul [722, 1991].A major source of symmetric indefinite linear systems is the least squaresproblem, because the augmented system is symmetric indefinite; see Chapter 19.
Other sources of such systems are interior methods for solving constrained optimization problems (see Forsgren, Gill, and Shinnerl [385, 1996],Turner [1030, 1991], and Wright [1115, 1992]) and linearly constrained optimization problems (see Gill, Murray, Saunders, and Wright [445, 1990], [446,1991 ]).The idea of using a block LDLT factorization with some form of pivotingfor symmetric indefinite matrices was first suggested by Kahan in 1965 [166,1971]. Bunch and Parlett [166, 1971] developed the complete pivoting strategyand Bunch [158, 1971] proved its stability.
Bunch [160, 1974] discusses a ratherexpensive partial pivoting strategy that requires repeated scalings. Bunch andKaufman [164, 1977] found the efficient partial pivoting strategy presentedhere, which is the one now widely used, and Bunch, Kaufman and Parlett [165,1976] give an Algol code implementing the diagonal pivoting method with thispivoting strategy.
Dongarra, Duff, Sorensen, and van der Vorst [315, 1991,§5.4.5] show how to develop a partitioned version of the diagonal pivoting226C HOLESKY FACTORIZATIONmethod with partial pivoting.Liu [709, 198 7] shows how to incorporate a threshold into the Bunch–Kaufinan partial pivoting strategy for sparse symmetric matrices; see also Duffet al. [326, 1991].
The partial pivoting strategy and variants of it describedby Bunch and Kaufman [164, 1977] do not preserve band structure, but thefill-in depends on the number of 2 × 2 pivots, which is bounded by the numberof negative eigenvalues (see Problem 10.11). Jones and Patrick [615, 1993],[616, 1994] show how to exploit this fact.The complete and partial pivoting strategies of Bunch et al.
use a fixednumber of tests to determine each pivot. Another possibility is to prescribegrowth bounds corresponding to 1 × 1 and 2 × 2 pivots and to search insome particular order for a pivot satisfying the bound. Fletcher [375, 1976]uses this approach to define a pivoting strategy that usually requires onlyO(n 2) operations. Duff, Reid, and co-workers apply the same approach to thediagonal pivoting method for sparse matrices, where sparsity considerationsalso influence the choice of pivot [331, 1979], [326, 1991]; their Fortran codesMA27 [329, 1982] and MA47 [330, 1995] implement the methods.Gill, Murray, Ponceleón, and Saunders [443, 1992] show how for sparse,symmetric indefinite systems the diagonal pivoting factorization can be usedto construct a (positive definite) preconditioned for an iterative method.Another method for solving symmetric indefinite systems is Aasen’s method[1, 1971], which employs the factorization PAPT = LTL T , where L is unitlower triangular and T is tridiagonal.
It is competitive with the diagonalpivoting method in terms of speed. Barwell and George [77, 1976] comparethe performance of Fortran codes for several methods for solving symmetric indefinite systems, including the diagonal pivoting method and Aasen’smethod.Dax and Kaniel [270, 1977] propose computing a factorization PAPT =LDL T for symmetric indefinite matrices by an extended form of Gaussianelimination in which extra row operations are used to “build up” a pivot element prior to the elimination operations; here, L is unit lower triangular andD is diagonal. A complete pivoting strategy for determining the permutationP is described in [270, 1977] and partial pivoting strategies in Dax [268, 1982].Analogues of the factorization for symmetric matrices exist for skewsymmetric matrices; see Bunch [161, 1982].Bunch [159, 1971] shows how to scale a symmetric matrix so that in everynonzero row and column the largest magnitude of an element is 1.10.6.1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
