Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 43
Текст из файла (страница 43)
The vector Dx =H - 1D – 1b is likely to have components that do not vary much in magnitude.Theorem 10.6 then guarantees that we obtain the components of Dx to goodrelative accuracy and this means that the components of x (which will varygreatly in magnitude) are obtained to good relative accuracy.So far, our results have all contained the proviso that Cholesky factorization runs to completion—in other words, the results assume that the argumentof the square root is always positive. Wilkinson [1092, 1968] showed that success is guaranteed if 20n3 / 2 κ 2 (A )u < 1, that is, if A is not too ill conditioned.It would be nice to replace A in this condition by H , where A = DHD. Justification for doing so is that Algorithm 10.2 is scale invariant, in the sense thatif we scale A FAF, where F is diagonal, then R scales to RF; moreover, ifF comprises powers of the machine base, then even the rounding errors scaleaccording to F.
The following theorem gives an appropriate modification ofWilkinson’s condition.Theorem 10.7 (Demmel). Let A = DHDbe symmetric positivedefinite, when D = diag(A )1/2. If min(H) > nγn +1/(1– γn +1) then Choleskyfactorization applied to A succeeds (barring underflow and overflow) and produces a nonsingularI f min(H) < –nγn +1/(1– γ n+1) then the computationis certain to fail.10.2 S ENSITIVITYOF THEC HOLESKY FACTORIZATION209Proof.
The proof is by induction. Consider Algorithm 10.2. The firststage obviously succeeds and gives> 0, since a 11 > 0. Suppose thealgorithm has successfully completed k – 1 stages, producing a nonsingularR k-1, and consider equations (10.1)-(10.3) with n replaced by k. The kthstage can be completed, but may give a pure imaginary(it will if fl(a –< 0). However, in the latter event, the error analysis of Theorem 10.5is still valid! Thus we obtainsatisfying= Ak + ∆A k, |∆Ak| < (1 –whereNow, with D k = diag(d k),we haveusing the interlacing of the eigenvalues [470, 1989, Cor.
8.1.4] and the condition of the theorem. Henceis positive definite, andtherefore so is the congruent matrix Ak + ∆A k, showing thatmust be realand nonsingular, as required to complete the induction.If Cholesky succeeds, then, by Theorem 10.5, D –1 (A + ∆A)D –1 is positivedefinite and so 0 <Hence if min(H) < -nγ n+1/(1– γ n +1) then the computation must fail.Note that, since ||H||2 > 1, the condition for success of Cholesky factorization can be written as κ 2 (H)nγn +1/(1– γ n + 1 ) < 1 .10.2.
Sensitivity of the Cholesky FactorizationThe Cholesky factorization has perturbation bounds that are similar to thosefor LU factorization, but of a simpler form thanks to the positive definiteness(||A- 1 ||2 replaces ||U- 1||2 ||L - 1 ||2 in the normwise bounds).Theorem 10.8 (Sun). Letbe symmetric positive definite with theCholesky factorization A = RTR and let AA be a symmetric matrix satisfying||A- 1 ∆A|| 2 < 1. Then A + ∆A has the Cholesky factorization A + ∆A =(R + ∆R)T(R + ∆R), where210Moreover, ifC HOLESKY< 1, whereF ACTORIZATION= ( R + ∆R) - T ∆A(R + ∆R)-1, then| ∆R| < triuwhere triu(·) denotes the upper triangular part.Note that the Cholesky factor of Ak = A(1:k,1:k) is Rk, and κ 2 (A k+l) >κ2 (A k) by the interlacing property of the eigenvalues.
Hence if Ak +1 (andhence A) is ill conditioned but Ak is well conditioned then Rk will be relativelyinsensitive to perturbations in A but the remaining columns of R will be muchmore sensitive.10.3. Positive Semidefinite MatricesIf A is symmetric and positive semidefinite (xTAx > 0 for all x) then aCholesky factorization exists, but the theory and computation are more subtlethan for positive definite A.The questions of existence and uniqueness of a Cholesky factorization areanswered by the following result.Theorem 10.9. Letbe positive semidefinite of rank r.
(a) Thereexists at least one upper triangular R with nonnegative diagonal elements suchthat A = RTR. (b) There is a permutation Π such that Π T AP has a uniqueCholesky factorization, which takes the form(10.11)where R11 is r × r upper triangular with positive diagonal elements.Proof. (a): Let the symmetric positive semidefinite square root X of Ahave the QR factorization X = QR with rii > 0. Then A = X 2 = XTX =R T Q T QR = R T R. (b): The algorithm with pivoting described below amountsto a constructive proof.Note that the factorization in part (a) is not in general unique. For example,For practical computations a factorization of the form (10.11) is needed,because this factorization so conveniently displays the rank deficiency.
Sucha factorization can be computed using an outer product Cholesky algorithm,comprising r = rank(A) stages. At the kth stage, a rank-1 matrix is subtracted from A so as to introduce zeros into positions k:n in the k th row and10.3 POSITIVE SEMIDEFINITE MATRICES211column. Ignoring pivoting for the moment, at the start of the kth stage wehave(10.12)where= [0, . . . , 0, rii , . . ., rin ]. The reduction is carried one stage furtherby computingOverall we have,To avoid breakdown whenvanishes (or is negative because of roundingerrors), pivoting is incorporated into the algorithm as follows. At the startof the kth stage an element> 0 (s > k) is selected as pivot, and rowsand columns k and s of Ak, and the kth and sth elements of ri , i = 1:k – 1,are interchanged.
The overall effect is to compute the decomposition (10.11),where the permutation Π takes account of all the interchanges.The standard form of pivoting is defined byThis is equivalent to complete pivoting in GE, since Ak is positive semidefiniteso its largest element lies on the diagonal. We note for later reference thatthis pivoting strategy produces a matrix R that satisfies (cf.
Problem 18.5)(10.13)It will be convenient to denote by cp(A) := ΠT AΠ the permuted matrixobtained from the Cholesky algorithm with complete pivoting.10.3.1. Perturbation TheoryIn this section we analyse, for a positive semidefinite matrix, the effect onthe Cholesky factorization of perturbations in the matrix.
This perturbationtheory will be used in the error analysis of the next section.C HOLESKY FACTORIZATION212Throughout this section A is assumed to be an n × n positive semidefinitematrix of rank r whose leading principal submatrix of order r is positivedefinite. For k = 1:r we will write(10.14)and other matrices will be partitioned conformably.We have the identity(10.15)where R11 is the Cholesky factor of A 11, R12 =andis the Schur complement of A11 in A. Note that Sr( A)0, so that for k = r,(10.15) is the (unique) Cholesky factorization of A. The following lemmashows how S k(A) changes when A is perturbed.Lemma 10.10. Let E be symmetric and assume A 11 + E 11 is nonsingular.Then(10.16)where W =Proof. We can expandThe result is obtained by substituting this expansion into Sk( A+E) = (A 22 +E22) – (A 12 + E1 2 ) T (A11 + E1 1 ) – 1 ( A 12 + E12), and collecting terms.Lemma 10.10 shows that the sensitivity of S k( A ) to perturbations in Ais governed by the matrix W =The question arises of whether,for a given A, the potentialmagnification of E indicated by (10.16)is attainable.
For the no-pivoting strategy, Π = I, the answer is trivially“yes”, since we can take E =with |γ| small, to obtain ||Sk (A+E) For complete pivoting, however, the answeris complicated by the possibility that the sequence of pivots will be different forA+ E than for A, in which case Lemma 10.10 is not applicable.
Fortunately, amild assumption on A is enough to rule out this technical difficulty, for small||E|| 2. In the next lemma we redefine A := cp(A ) in order to simplify thenot at ion.10.3 P OSITIVE S EMIDEFINITE M ATRICES213Lemma 10.11. Let A := cp(A). Suppose that( S i ( A ) ) 11 > (S i ( A ))j j, j = 2: n - i , i = 0: r- 1(10.17)(where S 0 (A) := A). Then, for sufficiently small ||E||2 , A+E = cp(A + E) .For E =with |γ| sufficiently small,Proof.
Note that since A = cp(A ), (10.17) simply states that thereare no ties in the pivoting strategy (since (S i ( A ) ) 11in (10.12)).Lemma 10.10 shows that Si (A+E) = Si (A) + O(||E|| 2), and so, in view of(10.17), for sufficiently small ||E||2 we haveThis shows that A + E = c p(A+E). The last part then follows fromLemma 10.10.We now examine the quantity ||W|| 2 =We show first that||W|| 2 can be bounded in terms of the square root of the condition number ofA1 1 .Lemma 10.12.
If A, partitioned as in (10.14), is symmetric positive definiteand A11 is positive definite thenProof. Writetogether withfrom the fact that the Schur complementinite.which followsis positive semidef-Note that, by the arithmetic-geometric mean inequality(x, y > 0), we also have, from Lemma 10.12,||A 2 2 || 2 )/2.The inequality of Lemma 10.12 is attained for the positive semidefinitematrixwhere Ip,q is the p × q identity matrix. This example shows that ||W||2 canbe arbitrarily large. However, for A := cp(A ), ||W||2 can be bounded solelyin terms of n and k.
The essence of the proof, in the next lemma, is thatlarge elements inare countered by small elements in A12. Hereafter weset k = r, the value of interest in the following sections.214C HOLESKY FACTORIZATIONL e m m a 1 0 . 1 3 . Let A := cp(A) and set k = r. Then(10.18)There is a parametrized family of rank-r matrices A(θ ) = cp(A(θ )),for whichProof.
The proof is a straightforward computation. The matrix A(O) :=R(θ)TR(θ), where(10.19)with c = cosθ, s = sin θ. This is the r × n version of Kahan’s matrix (8.10). Rsatisfies the inequalities (10.13) (as equalities) and so A(θ ) = cp(A(θ )).We conclude this section with a “worst-case” example for the Cholesky factorization with complete pivoting. Let U( θ) = diag(r, r–1, . . . . 1)R(θ), whereR(θ) is given by (10.19), and define the rank-r matrix C(θ) = U(θ) T U(θ).Then C(θ) satisfies the conditions of Lemma 10.11. Also,Thus, from Lemma 10.11, for E =with |γ| and θ sufficiently small,This example can be interpreted as saying that in exact arithmetic the residual after an r-stage Cholesky factorization of a semidefinite matrix A canoverestimate the distance of A from the rank-r semidefinite matrices by afactor as large as (n – r)(4r – 1)/3.10.3.2. Error AnalysisIn this section we present a backward error analysis for the Cholesky factorization of a positive semidefinite matrix.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
