Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 42
Текст из файла (страница 42)
Show that forthreshold pivotingwhere µ j is the number of nonzero entries in the jth column of U. Henceobtain a bound for p n .PROBLEMS2019.13. (RESEARCH PROBLEM) Obtain sharp bounds for the growth factor forGE with partial pivoting applied to (a) a matrix with lower bandwidth pand upper bandwidth q (thus generalizing Theorem 9.10), and (b) a quasitridiagonal matrix (an n × n matrix that is tridiagonal except for nonzero(1, n) and (n, 1) elements).9.14. (R ESEARCH P ROBLEM ) Explain why the growth factor for GE withpartial pivoting is almost always small in practice.9.15. (R ESEARCH P ROBLEM ) For GE with complete pivoting what is thevalue of limng (n)/n (see (9.14))? Isequal to n for Hadamard matrices?Previous HomeChapter 10Cholesky FactorizationThe matrix of that equation system is negative definite-which is apositive definite system that has been multiplied through by – 1.For all practical geometries the common finite differenceLaplacian operator gives rise to these,the best of all possible matrices.Just about any standard solution method will succeed,and many theorems are available for your pleasure.—FORMAN S.
ACTON, Numerical Methods That Work (1970)Many years ago we made out of half a dozen transformersa simple and rather inaccurate machine forsolving simultaneous equations—the solutions beingrepresented as flux in the cores of the transformers.During the course of our experiments weset the machine to solve the equations—X + Y + Z = 1X + Y + Z = 2X + Y + Z = 3The machine reacted sharply—it blew the main fuse and put all the lights out.—B. V. BOWDEN, The Organization of a Typical Machine (1953)There does seem to be some misunderstanding about thepurpose of an a priori backward error analysis.All too often, too much attention is paidto the precise error bound that has been established.The main purpose of such an analysis is either toestablish the essential numerical stability of an algorithm or toshow why it is unstable and in doing so toexpose what sort of change is necessary to make it stable.The precise error bound is not of great importance.—J.
H. WILKINSON, Numerics/ Linear Algebra on Digits/ Computers (1974)203NextCHOLESKY FACTORIZATION20410.1. Symmetric Positive DefiniteSymmetric positive definiteness is one of the highest accolades to which amatrix can aspire. Symmetry confers major advantages and simplificationsin the eigenproblem and, as we will see in this chapter, positive definitenesspermits economy and numerical stability in the solution of linear systems.A symmetrical matrixis positive definite if xTAx > 0 for allnonzeroWell-known equivalent conditions to A = AT being positivedefinite are• det(A k) > 0, k = 1:n, where Ak = A(1:k, 1:k) is the leading principalsubmatrix of order k.•>0, k = 1:n, wheredenotes the kth largest eigenvalue.The first of these conditions implies that A has an LU factorization, A = LU(see Theorem 9.1).
Another characterization of positive definiteness is that thepivots in LU factorization are positive, since u kk = det(A k)/det(A k–1). Byfactoring out the diagonal of U and taking its square root, the LU factorizationcan be converted into a Cholesky factorization: A = RTR, where R is uppertriangular with positive diagonal elements. This factorization is so importantthat it merits a direct proof.Theorem 10.1. Ifunique upper triangularA = R T R.is symmetric positive definite then there is awith positive diagonal elements such thatProof.
The proof is by induction. The result is clearly true for n =1. Assume it is true for n – 1. The leading principal submatrix An –1 =A(1:n–1,1:n–1) is positive definite, so it has a unique Cholesky factorizationA n–1 =We have a factorization(10.1)if(10.2)(10.3)Equation (10.2) has a unique solution since Rn-1 is nonsingular. Then (10.3)gives β 2 = a – rTr. It remains to check that there is a unique real, positive βsatisfying this equation. From the equation0 < det(A) = det(R T) det(R) = det(R n- 1 ) 2 β 210.1 S YMMETRIC P OSITIVE D EFINITE M ATRICES205we see that β 2 >0, hence there is a unique β > 0.The proof of the theorem is constructive, and provides a way to computethe Cholesky factorization that builds R a column at a time.
Alternatively,we can work directly from the equationswhich follow by equating (i, j) elements in A = RTR. By solving these equations in the order (1,1), (1,2), (2,2), (1,3), (2,3), (3,3), . . . , (n,n), we obtainthe following algorithm.Algorithm 10.2. Given a symmetric positive definiterithm computes the Cholesky factorization A = RTR.this algofor j = 1:nfor i = 1:j–1endendCost: n 3/3 flops (half the cost of LU factorization).As for Gaussian elimination (GE), there are different algorithmic forms ofCholesky factorization. Algorithm 10.2 is the jik or “sdot” form.
We describethe kij, outer product form in §10.3.Given the Cholesky factorization A = RTR, a linear system Ax = b canbe solved via the two triangular systems RTy = b and Rx = y.If we define D =then the Cholesky factorization A = R T RTcan be rewritten as A = LDL , where L = RT diag(r i i )–1 is unit lowertriangular. The LDLT factorization is sometimes preferred over the Choleskyfactorization because it avoids the need to compute the n square roots thatdetermine the rii .
The LDLT factorization is certainly preferred for solvingtridiagonal systems, as it requires n less divisions than Cholesky factorizationin the substitution stage. All the results for Cholesky factorization in thischapter have analogues for the LDLT factorization. Block LDLT factorizationfor indefinite matrices is discussed in §10.4.10.1.1. Error AnalysisError bounds for Cholesky factorization are derived in a similar way to thosefor LU factorization.
Consider Algorithm 10.2. Using Lemma 8.4 we have(10.4)206C HOLESKY FACTORIZATIONFrom a variation of Lemma 8.4 in which the division is replaced by a squareroot (see Problem 10.3), we haveA backward error result is immediate.Theorem 10.3. If Cholesky factorization applied to the symmetric positiveruns to completion then the computed factor Rdefinite matrixsatisfies(10.5)Theorem 10.4.
Letbe symmetric positive definite and supposeCholesky factorization produces a computed factor R and a computed solutionto Ax = b. Then(10.6)Proof. The proof is analogous to the proof of Theorem 9.4.These results imply that Cholesky factorization enjoys perfect normwisebackward stability. The key inequality iswhose analogue for the computedis, from (10.5),Thus (10.6) implies(10.7)where for the last inequality we have assumed that nγ n +1 < 1/2. Anotherindicator of stability is that the growth factor for GE is exactly 1 (see Problem 10.4). It is important to realize that the multipliers can be arbitrarilylarge (consider, for example,as θ0).
But, remarkably, for a positivedefinite matrix the size of the multipliers has no effect on stability.Note that the perturbation ∆A in (10.6) is not symmetric, in general,because the backward error matrices for the triangular solves with R andRT are not the transposes of each other. For conditions guaranteeing that a“small” symmetric ∆A can be found, see Problem 7.11.The following rewritten version of Theorem 10.3 provides further insightinto Cholesky factorization.10.1 S YMMETRIC P OSITIVE D EFINITE M ATRICES207Theorem 10.5 (Demmel).
If Cholesky factorization applied to the symmetric positive definite matrixruns to completion then the computedfactorsatisfieswhere di =Proof. Theorem 10.3 shows that= A+∆A with |∆ A| < γn + 1Denoting bythe ith column ofwe haveso thatThen, using the Cauchy–Schwarz inequality,giving(10.8)and the required bound for ∆A.Standard perturbation theory applied to Theorem 10.4 yields a boundof the formHowever, with the aid ofTheorem 10.5 we can obtain a potentially much smaller bound.
The idea isto write A = DHD where D = diag(A )1/2, so that H has unit diagonal. vander Sluis’s result (Corollary 7.6) shows that(10.9)so D is nearly a condition-minimizing diagonal scaling. It follows that κ 2 (H) <nκ 2 (A) and that κ 2 (H) << κ 2 (A) is possible if A is badly scaled. Note that1 < ||H||2 < n, since H is positive definite with hii1.Theorem 10.6 (Demmel, Wilkinson).
Let A = DHDbe symmetricpositive definite, where D = diag(A )1/2, and suppose Cholesky factorizationsuccessfully produces a computed solutionto Ax = b. Then the scaled errorsatisfies(10.10)where= 2n(1 - γ n + l ) – 1 γ n + 1 .Proof. Straightforward analysis shows that (cf. the proof of Theorem 9.4)(A + ∆A) = b, where208C HOLESKYF ACTORIZATIONwith |∆A 1 | < (1–γn + 1 ) – 1 γ n+ 1ddT (by Theorem 10.5) and |∆ 1| < diag( γi )|∆2| < diag( γn- i+ 1)Scaling with D, we haveand standard perturbation theory gives- 1T -1But, using (10.8) and ||D dd D || 2=||eeT||2 = n, we haveusing Lemma 3.3, which yields the result.Care needs to be exercised when interpreting bounds that involve scaledquantities, but in this case the interpretation is relatively easy. Suppose thatH is well conditioned and κ 2 (D) is large, which represents the artificial illconditioning that the DHD scaling is designed to clarify.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
