Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 50
Текст из файла (страница 50)
Letandbe the computed block LU factors offrom Algorithm 12.2 (with Implementation 1), and letbe the computed solution to Ax = b. Under the assumptions(12.3), (12.13), and (12.14),(12.15)where the constant dn is commensurate with those in the assumptions.12.3 E RROR A NALYSISOFB LOCK LU FACTORIZATION251Proof. We omit the proof (see Demmel, Higham, and Schreiber [293,for details). It is similar to the proof of Theorem 12.3.1995]The bounds in Theorem 12.4 are valid also for other versions of block LUfactorization obtained by “block loop reordering”, such as a block gaxpy basedalgorithm.Theorem 12.4 shows that the stability of block LU factorization is determined by the ratio ||L||||U||/||A|| (numerical experiments show that thebounds are, in fact, reasonably sharp).
If this ratio is bounded by a modest function of n, then L and U are the true factors of a matrix close toA, andsolves a slightly perturbed system. However,can exceed||A|| by an arbitrary factor, even if A is symmetric positive definite or diusing theagonally dominant by rows. Indeed, ||L|| > ||L2 1 || =partitioning (12.2), and this lower bound for ||L|| can be arbitrarily large.In the following two subsections we investigate this instability more closelyand show that ||L||||U|| can be bounded in a useful way for particular classesof A. Without further comment we make the reasonable assumption that||L||U||||L||||U|| , so that these bounds maybe used in Theorem 12.4.What can be said for Implementation 2? Suppose, for simplicity, that theinverses(which are used in step 2 of Algorithm 12.2 and in the blockback substitution) are computed exactly.
Then the best bounds of the forms(12.13) and (12.14) areWorking from these results, we find that Theorem 12.4 still holds provided thefirst-order terms in the bounds in (12.15) are multiplied by max i κ(U i i ). Thissuggests that Implementation 2 of Algorithm 12.2 can be much less stablethan Implementation 1 when the diagonal blocks of U are ill conditioned, andthis is confirmed by numerical experiments.12.3.1. Block Diagonal DominanceOne class of matrices for which block LU factorization has long been knownto be stable is block tridiagonal matrices that are diagonally dominant inan appropriate block sense. A general matrixis block diagonallydominant by columns with respect to a given partitioning A = ( A i j ) and agiven norm if, for all j,(12.16)A is block diagonally dominant by rows if AT is block diagonally dominant bycolumns.
For the block size 1, the usual property of point diagonal dominance252B LOCK LU FACTORIZATIONis obtained. Note that for the 1- and co-norms diagonal dominance does notimply block diagonal dominance, nor does the reverse implication hold (seeProblem 12.2).
Throughout our analysis of block diagonal dominance we takethe norm to be an arbitrary subordinate matrix norm.First, we show that for block diagonally dominant matrices a block LUfactorization exists, using the key property that block diagonal dominance isinherited by the Schur complements obtained in the course of the factorization.In the analysis we assume that A has m block rows and columns.Theorem 12.5 (Demmel, Higham, and Schreiber). Supposeisnonsingular and block diagonally dominant by rows or columns with respect toa subordinate matrix norm in (12.16). Then A has a block LU factorization,and all the Schur complements arising in Algorithm 12.2 have the same kindof diagonal dominance as A.Proof.
This proof is a generalization of Wilkinson’s proof of the corresponding result for point diagonally dominant matrices [1085, 1961, pp. 288–289], [470, 1989, p. 120] (as is the proof of Theorem 12.6 below). We considerthe case of block diagonal dominance by columns; the proof for row-wise diagonal dominance is analogous.The first step of Algorithm 12.2 succeeds, since A 11 is nonsingular, producing a matrix that we can write asFor j = 2:m we haveusing (12.16)using (12.16),12.3 E RROR A NALYSISOF253B LOCK LU FACTORIZATION(12.17)Now ifis singular it follows thathence also A, is singular, which is a contradiction. Thusand (12.17) can be rewrittentherefore A(2), andis nonsingular,showing that A(2) is block diagonally dominant by columns.
The result followsby induction.The next result allows us to bound ||U|| for a block diagonally dominantmatrix.Theorem 12.6 (Demmel, Higham, and Schreiber). Let A satisfy the conditions of Theorem 12.5. If A(k) denotes the matrix obtained after k – 1 stepsof Algorithm 12.2, thenProof. Let A be block diagonally dominant by columns (the proof for rowdiagonal dominance is similar).
Thenusing (12.16). By induction, using Theorem 12.5, it follows thatThis yieldsB LOCK LU FACTORIZATION254The implications of Theorems 12.5 and 12.6 for stability are as follows.Suppose A is block diagonally dominant by columns. Also, assume for themoment that the (subordinate) norm has the property that(12.18)which holds for any p-norm, for example. The subdiagonal blocks in the firstblock column of L are given by Li1 =and so< 1, by(12.16) and (12.18).
From Theorem 12.5 it follows that<1 for j = 2:m. Since Uij =for j > i, Theorem 12.6 shows that ||Uij|| <2||A|| for each block of U (and ||Uij|| < ||A||). Therefore ||L|| < m and ||U|| <m2 ||A||, and so ||L||||U|| < m 3 ||A|| . For particular norms the bounds on theblocks of L and U yield a smaller bound for ||L|| and ||U||. For example, forthe 1-norm we have ||L||1 ||U||1 < 2m||A||1 and for the -normWe conclude that block LU factorization is stable if A is blockdiagonally dominant by columns with respect to any subordinate matrix normsatisfying (12.18).Unfortunately, block LU factorization can be unstable when A is blockdiagonally dominant by rows, for although Theorem 12.6 guarantees that||U i j || < 2||A||, ||L|| can be arbitrarily large.
This can be seen from theexamplewhere A is block diagonally dominant by rows in any subordinate norm forany nonsingular matrix A11. It is easy to confirm numerically that block LUfactorization can be unstable on matrices of this form.Next, we bound ||L||||U|| for a general matrix and then specialize to pointdiagonal dominance. From this point on we use the norm ||A|| := max i,j |aij|.We partition A according to(12.19)and denote by pn the growth factor for GE without pivoting. We assume thatGE applied to A succeeds.To bound ||L||, we note that, under the partitioning (12.19), for the firstblock stage of Algorithm 12.2 we have ||L2 1 || =< n pn κ(A) (seeProblem 12.4).
Since the algorithm works recursively with the Schur complement S, and since every Schur complement satisfies κ( S) < pn κ(A) (seeProblem 12.4), each subsequently computed subdiagonal block of L has normat mostSince U is composed of elements of A together with elements of Schur complements of A,||U|| < pn||A||.(12.20)12.3 E RROR ANALYSISOFB LOCK LU FACTORIZATION255Overall, then, for a general matrix(12.21)Thus, block LU factorization is stable for a general matrix A as long as GEis stable for A (that is, pn is of order 1) and A is well conditioned.If A is point diagonally dominant by columns then, since every Schurcomplement enjoys the same property, we have ||Lij || < 1 for i > j, byProblem 12.5.
Hence ||L|| = 1. Furthermore, pn < 2 (Theorem 9.8 or Theorem 12.6), giving ||U|| < 2||A|| by (12.20), and so||L||||U|| < 2||A||.Thus block LU factorization is perfectly stable for a matrix point diagonallydominant by columns.If A is point diagonally dominant by rows then the best we can do is totake pn < 2 in (12.21), obtaining||L||||U|| < 8nk(A)||A||.(12.22)Hence for point row diagonally dominant matrices, stability is guaranteed if Ais well conditioned. This in turn is guaranteed if the row diagonal dominanceamounts γj in the analogue of (12.16) for point row diagonal dominance aresufficiently large relative to ||A||, because< (minj γj )–1 (see problem 8.7(a)).12.3.2.
Symmetric Positive Definite MatricesFurther useful results about the stability of block LU factorization can bederived for symmetric positive definite matrices. First, note that the existenceof a block LU factorization is immediate for such matrices, since all theirleading principal submatrices are nonsingular. Let A be a symmetric positivedefinite matrix, partitioned asThe definiteness implies certain relations among the submatrices Aij that canbe used to obtain a stronger bound for ||L|| 2 than can be deduced for a generalmatrix (cf. Problem 12.4).Lemma 12.7. If A is symmetric positive definite thenProof.
This lemma is a corollary of Lemma 10.12, but we give a separateproof. Let A have the Cholesky factorization256BLOCK LU FACTORIZATIONTable 12.1. Stability of block and point L U factorization. p n is the growth factor forGE without pivoting.Matrix propertySymmetric positive definiteBlock column diagonally dominantPoint column diagonally dominantBlock row diagonally dominantPoint row diagonally dominantArbitraryBlock LU Point LUκ(A) 1 / 211pn11pnκ(A)1pnThe following lemma is proved in a way similar to the second inequality inProblem 12.4.Lemma 12.8.