Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 44
Текст из файла (страница 44)
An important consideration is that10.3 P OSITIVE S EMIDEFINITE M ATRICES215a matrix of floating point numbers is very unlikely to be ‘(exactly” positivesemidefinite; errors in forming or storing A will almost certainly render thesmallest eigenvalue nonzero, and possibly negative. Therefore error analysisfor a rank r positive semidefinite matrix may appear, at first sight, to be oflimited applicability. One way around this difficulty is to state results forA = A + ∆A, where A is the matrix stored on the computer, A is positivesemidefinite of rank r, and ∆A is a perturbation, which could represent therounding errors in storing Ã, for example.
However, if the perturbation AAis no larger than the backward error for the Cholesky factorization, then thisextra formalism can be avoided by thinking of ∆A as being included in thebackward error matrix. Hence for simplicity, we frame the error analysis fora positive semidefinite A.The analysis makes no assumptions about the pivoting strategy, so Ashould be thought of as the pre-permuted matrix Π T AΠ.Theorem 10.14. Let A be an n × n symmetric positive semidefinite matrixof floating point numbers of rank r < n. Assume that A11 = A(1:r, 1:r) ispositive definite with(10.20)where A11 = D 11 H11 D 11 and D11 = diag(A 11)1/2.
Then, in floating pointarithmetic, the Cholesky algorithm applied to A successfully completes r stages(barring underflow and overflow), and the computed r × n Cholesky factorsatisfies(10.21)where W =Proof. First, note that condition (10.20) guarantees successful completionof the first r stages of the algorithm by Theorem 10.7.Analysis very similar to that leading to Theorem 10.3 shows that(10.22)whereand(10.23)216CHOLESKY FACTORIZATIONTaking norms in (10.23) and using the inequalitywe obtainwhich implies(10.24)Our aim is to obtain an a priori bound forIt is clear from( r +1 )(10.22)-(10.24) that to do this we have only to bound ||Â||2.
To this end,we interpret (10.22) in such a way that the perturbation theory of §10.3.1 maybe applied.Equation (10.22) shows thatis the true Schur complement for thematrix A + E and thatis positive definite. Hence we canuse Lemma 10.10 to deduce thatSubstituting from (10.24) we find thatFinally, using (10.22) and (10.24), we obtainTheorem 10.14 is just about the best result that could have been expected,because the bound (10.21) is essentially the same as the bound obtained ontaking norms in Lemma 10.10. In other words, (10.21) simply reflects theinherent mathematical sensitivity of A–RTR to small perturbations in A.We turn now to the issue of stability.
Ideally, for A as defined in Theorem 10.14, the computed Cholesky factor Rr produced after r stages of thealgorithm would satisfy10.3 P OSITIVE S EMIDEFINITE M ATRICES217where cn is a modest constant. Theorem 10.14 shows that stability dependson the size of ||W||2 =(to the extent that ||W||2 appears in arealistic bound for the backward error).If no form of pivoting is used then ||W|| 2 can be arbitrarily large for fixedn (see §10.3.1) and the Cholesky algorithm must in this case be classed asunstable.
But for complete pivoting we have from Lemma 10.13 the upperbound ||W|| 2 < (1/3(n – r)(4 r – 1))1/2. Thus the Cholesky algorithm withcomplete pivoting is stable if r is small, but stability cannot be guaranteed,and seems unlikely in practice, if ||W|| 2 (and hence, necessarily, r and n) islarge.Numerical experiments show that ||W||2 is almost always small in practice(typically less than 10) [540, 1990]. However, it is easy to construct exampleswhere ||W||2 is large. For example, if R is a Cholesky factor of A from completepivoting then let C = M(R ) T M( R), where M( R) is the comparison matrix;C will usually have a much larger value of ||W|| 2 than A.An important practical issue is when to terminate the Cholesky factorization of a semidefinite matrix. The LINPACK routine xCHDC proceeds withthe factorization until a nonpositive pivot is encountered, that is, up to andincluding stage k – 1, where k is the smallest integer for which(10.25)Usually k > r + 1, due to the effect of rounding errors.A more sophisticated termination criterion is to stop as soon as(10.26)for some readily computed norm and a suitable toleranceThis criterionterminates as soon as a stable factorization is achieved, avoiding unnecessarywork in eliminating negligible elements in the computed Schur complementNote thatis indeed a reliable order-of-magnitude estimate of thetrue residual, sinceis the only nonzero block of Â(k) and, by (10.22) andwith ||E|| = O(u)(||A|| + ||Â(k)||).(10.24),Another possible stopping criterion is(10.27)This is related to (10.26) in that if A (pre-permuted) and Âk are positivesemidefinite then= maxi,j |aij|||A|| 2, and similarly maxNote that (10.27) boundssince if (10.27) holds first at thekth stage then, using Theorem 8.13,218C HOLESKY FACTORIZATIONPractical experience shows that the criteria (10.26) and (10.27) with=nu both work well, and much more effectively than (10.25) [540, 1990].
Wefavour (10.27) because of its negligible cost.10.4. Symmetric Indefinite Matrices and Diagonal Pivoting MethodLetbe symmetric but indefinite, that is, (x T A x) (y T Ay) < 0 forsome x and y. How can we solve Ax = b efficiently?Gaussian elimination with partial pivoting (GEPP) can be used to compute the factorization PA = LU, but it does not take advantage of the symmetry to reduce the cost and storage. We might try to construct a factorizationA = L D L T, where L is unit lower triangular and D is diagonal.
But thisfactorization may not exist, even if symmetric pivoting is allowed, and if itdoes exist its computation may be unstable. For example, considerThere is arbitrarily large element growth for 0 <<< 1, and the factorizationdoes not exist for = 0.The most popular approach for solving symmetric indefinite systems is touse a block LDLT factorizationP A P T = L D L T,where L is unit lower triangular and D is block diagonal with 1 × 1 or 2 × 2diagonal blocks. This factorization is essentially a symmetric block form ofGE, with pivoting.
Note that by Sylvester’s inertia theorem, A and D havethe same inertia13, which is easily determined from D (see Problem 10.11).To begin the computation of the factorization we choose a permutation Πand an integer s = 1 or 2 so thatwith E nonsingular. Then we compute the factorization13The inertia of a symmetric matrix is an ordered triple {i +, i–, i0}, where i+ is thenumber of positive eigenvalues, i– the number of negative eigenvalues, and i0 the numberof zero eigenvalues.21910.4 I NDEFINITE M ATRICESThis process is repeated on the (n – s) × (n – s) Schur complementà = B – CE- 1C T .The cost of the method is n 3/3 flops (the same as the cost of Cholesky factorization of a positive definite matrix) plus the cost of determining the permutations 17.
This method for computing the block LDLT factorization iscalled the diagonal pivoting method. It can be thought of as a generalizationof Lagrange’s method for reducing a quadratic form to diagonal form (devisedby Lagrange in 1759 and rediscovered by Gauss in 1823) [763, 1961, p. 371].One conceivable difficulty with the diagonal pivoting method can be disposed of immediately. If a nonsingular pivot matrix E of dimension 1 or 2cannot be found, then all 1 × 1 and 2 × 2 principal submatrices of the symmetric matrix A are singular, and this is easily seen to imply that A is thezero matrix.The strategy for choosing 17 is crucial for achieving stability. A suitablemodification of the error analysis for block LU factorization (Theorem 12.4)tells us that, provided linear systems involving 2 x 2 pivots are solved in anormwise backward-stable way, the condition ||L|| ||D|| ||LT|| < cn||A||, fora modest constant c n , is sufficient to ensure stability.
A key requirement,therefore, is to choose the pivot E so that the Schur complement à is suitablybounded, since D is made up of elements of Schur complements. We describetwo suitable pivoting strategies.10.4.1. Complete PivotingBunch and Parlett [166, 1971] devised the following strategy for choosing 17.It suffices to describe the interchanges for the first stage of the factorization.Let µ 0 = maxi,j |aij|, µ1 = maxi |aii|, and choose aif µ 1 > aµ 0Set s = 1, and choose Π so that |e 11| = µ 1 .elseSet s = 2, and choose Π so that |e 21| = µ 0 .end(0, 1).Note that µ 1 is the best 1 × 1 pivot under symmetric permutations andµ 0 is the pivot that would be chosen by GE with complete pivoting.
Thisstrategy says “as long as there is a diagonal pivot element not much smallerthan the complete pivot, choose it as a 1 × 1 pivot”, that is, “choose a 1 × 1pivot whenever possible”. If the strategy dictates the use of a 2 × 2 pivot thenthat pivot E is indefinite (see Problem 10.11).It remains to determine a. This is done by minimizing a bound on theelement growth. For the following analysis we assume that the interchanges220C HOLESKYF ACTORIZATIONhave already been done. If s = 1 thenNow consider the case s = 2. The (i, j) element of the Schur complementà = B – CE - 1 CT is(10.28)Nowand, using the symmetry of E,Since a(0, 1), we have |det(E)| > (1 – a 2 )ThusSince |cij| < µ 0, we obtain from (10.28)To determine a we equate the maximum growth for two s = 1 steps withthat for one s = 2 step:which reduces to the quadratic equation 4a 2 – a – 1 = 0.
We require thepositive rootThe analysis guarantees a growth factor bound of (1+a- 1 ) n -1 = (2.57) n - 1 .This bound is pessimistic, however; a much more detailed analysis by Bunch[158, 1971] shows that the growth factor is no more than 3.07( n –1)0.446 timeslarger than the bound (9.13) for LU factorization with complete pivoting—avery satisfactory result. Strictly speaking, bounding the growth factor boundsonly ||D||, not ||L||.
But it is easy to show that for s = 1 and 2 no element ofCE-1 exceeds max{1/a,1/(1–a)} in absolute value, and so ||L|| is boundedindependently of A.Since complete pivoting requires the whole active submatrix to be searchedat each stage, it needs up to n3/6 comparisons, and so the method is ratherexpensive.22110.4 I NDEFINITE M ATRICES10.4.2. Partial PivotingBunch and Kaufman [164, 1977] devised a pivoting strategy for the diagonal pivoting method that requires only O( n 2) comparisons.
At each stageit searches at most two columns and so is analogous to partial pivoting forLU factorization. The strategy contains several logical tests. As before, wedescribe the pivoting for the first stage only. Recall that s denotes the size ofthe pivot block.Choose a(0, 1).If = 0 there is nothing to do on this stage of the elimination.r := min{i > 2: |ai1| =if |a 11| >(1)s = 1, Π = Ielseif |a1 1 |σ > a 2(2)s = 1, Π = Ielse if |arr| > aσ(3)s = 1 and choose Π to swap rows and columns 1 and r.else(4)s = 2 and choose Π to swap rows and columns 2 and r,so that |(ΠAΠ T )21| =endendTo understand the algorithm it helps to consider the matrixand to note that the pivot is one of a 11, arr , and= |ar1|, this matrix with replaced by ar1 ).(or, rather, since222C HOLESKY FACTORIZATIONTo bound the element growth reconsider each case in turn, noting thatfor cases (1) and (2) the elements of the Schur complement are given by14Case (l):Case (2): Using symmetry,Case (3): The original arr is now the pivot, and |arr| > a σ, soCase (4): This is where we use a 2 x 2 pivot, which, after the interchanges,is E = (Π T AΠ)(1:2,1:2) =NowThe elements of the Schur complement à = B-CE- 1 C T a r e g i v e n b ysoThis analysis shows that the bounds for the element growth for s = 1 ands = 2 are the same as the respective bounds for the complete pivoting strategy.Hence, using the same reasoning, we again choose a =The growth factor for the partial pivoting strategy is bounded by (2.57)n - 1 .As for GEPP, large growth factors do not seem to occur in practice.
But unlike for GEPP, no example is known for which the bound is attained [164,1977]; see Problem 10.18.14We commit a minor abuse of notation, in that in the rest of this section ãij shouldreally be ãi– 1, j–1 (s = 1) or ãi- 2 ,j –2 (s = 2).10.5 N ONSYMMETRIC P OSITIVE D EFINITE M ATRICES223As noted in the previous subsection, a bound on the growth factor p n doesnot in itself ensure stability. Indeed although ||D||/||A|| is bounded for partialpivoting, ||L||/||A|| can be arbitrarily large for fixed n; see Problem 10.15.Higham [559, 1995] gives a detailed error analysis of the diagonal pivotingmethod with an arbitrary pivoting strategy, under the assumption that computed solutions to linear systems involving 2 × 2 pivots have a small componentwise relative backward error.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
