Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 51
Текст из файла (страница 51)
If A is symmetric positive definite then the Schur complementS = A 22 –satisfies κ2 (S) < κ 2 (A).Using the same reasoning as in the last subsection, we deduce from thesetwo lemmas that each subdiagonal block of L is bounded in 2-norm by κ 2 (A) 1 / 2 .Therefore ||L|| 2 < 1 + mκ 2 (A)1/2, where there are m block stages in the algorithm. Also, it can be shown that ||U|| 2 <Hence(12.23)It follows from Theorem 12.4 that when Algorithm 12.2 is applied to a symmetric positive definite matrix A, the backward errors for the LU factorizationand the subsequent solution of a linear system are both bounded by(12.24)Any resulting bound forwill be proportional to κ 2 (A)3/2, ratherthan κ2 (A) as for a stable method.
This suggests that block LU factorizationcan lose up to 50% more digits of accuracy in x than a stable method forsolving symmetric positive definite linear systems. The positive conclusion tobe drawn, however, is that block LU factorization is guaranteed to be stablefor a symmetric positive definite matrix that is well conditioned.The stability results for block LU factorization are summarized in Table 12.1, which tabulates a bound for ||A – LU||/(cnu||A||) for block and point12.4 N OTESANDR EFERENCES257LU factorization for the matrix properties considered in this chapter.
Theconstant cn incorporates any constants in the bound that depend polynomially on the dimension, so a value of 1 in the table indicates unconditionalstability.12.4. Notes and ReferencesThe distinction between a partitioned algorithm and a block algorithm israrely made in the literature (exceptions include the papers by Schreiber [902,1988] and Demmel, Higham, and Schreiber [293, 1995]); the term “block algorithm” is frequently used to describe both types of algorithm. A partitioned algorithm might also be called a “blocked algorithm” (as is done byDongarra, Duff, Sorensen, and van der Vorst [315, 1991]), but the similarity of this term to “block algorithm” can cause confusion and so we do notrecommend this terminology. Note that in the particular case of matrix multiplication, partitioned and block algorithms are equivalent.
Our treatment ofpartitioned LU factorization has focused on the stability aspects; for furtherdetails, particularly concerning implementation on high-performance computers, see Dongarra, Duff, Sorensen, and van der Vorst [315, 1991] and Goluband Van Loan [470, 1989].Block LU factorization appears to have first been proposed for block tridiagonal matrices, which frequently arise in the discretization of partial differential equations. References relevant to this application include Isaacsonand Keller [607, 1966, p.
59], Varah [1048, 1972], Bank and Rose [53, 1977],Mattheij [737, 1984], [738, 1984], and Concus, Golub, and Meurant [235, 1985].For an application of block LU factorization to linear programming, seeEldersveld and Saunders [351, 1992].Theorem 12.3 is from Demmel and Higham [291, 1992]. The results in§12.3 are from Demmel, Higham, and Schreiber [293, 1995], which extendsearlier analysis of block LU factorization by Demmel and Higham [291, 1992].Block diagonal dominance was introduced by Feingold and Varga [366,19 6 2 ], and has been used mainly in generalizations of the Gershgorin circletheorem.
Varah [1048, 1972] obtained bounds on ||L|| and ||U|| for blockdiagonally dominant block tridiagonal matrices; see Problem 12.1.Theorem 12.5 is obtained in the case of block diagonal dominance by rowswith minj γ j > 0 by Polman [837, 1987]; the proof in [837, 1987] makes use ofthe corresponding result for point diagonal dominance and thus differs fromthe proof we have given.At the cost of a much more difficult proof, Lemma 12.7 can be strengthenedto the attainable bound< (κ 2 ( A)1/2 – κ 2 (A)-1/2)/2, as shownby Demmel [279, 1983, Thin.
4], but the weaker bound is sufficient for ourpurposes.B LOCK LU FACTORIZATION25812.4.1. LAPACKLAPACK does not implement block LU factorization, but its LU factorization(and related) routines for full matrices employ partitioned LU factorizationin order to exploit the level-3 BLAS and thereby to be efficient on highperformance machines.Problems12.1. (Varah [1048, 1972]) Suppose A is block tridiagonal and has the blockLU factorization A = LU (so that L and U are block bidiagonal and Ui,i + 1 =A i , i +1).
Show that if A is block diagonally dominant by columns thenwhile if A is block diagonally dominant by rows thenWhat can be deduced about the stability of the factorization and has the blockclasses of matrices?12.2. Show that for the 1- and co-norms diagonal dominance does not implyblock diagonal dominance, and vice versa.is symmetric, has positive diagonal elements, and is block12.3. Ifdiagonally dominant by rows, must it be positive definite?12.4. Letbe partitioned(12.25)with A11 nonsingular.
Let ||A|| := max ij |aij|. Show thatn p n κ(A). Show that the Schur complement S = A 22 –κ (S) < pn κ(A).satisfies12.5. Letbe partitioned as in (12.25), with A11 nonsingular,and suppose that A is point diagonally dominant by columns. Show that< 1.12.6. Show that under the conditions of Theorem 12.3 the computed solutionto Ax = b satisfiesPROBLEMS259and the computed solution to the multiple right-hand side system AX = B(where (12.4) is assumed to hold for the multiple right-hand side triangularsolves) satisfiesIn both cases, cn is a constant depending on n and the block size.12.7. Let X =thatwhere A is square and nonsingular. Showdet(X) = det(A) det(D – CA- 1 B ) .Assuming A, B, C, D are all m × m, give a condition under which det(X ) =det( AD – CB).PreviousHomeChapter 13Matrix InversionIt is amusing to remark that we were so involved withmatrix inversion that we probably talked of nothing else for months.Just in this period Mrs.
von Neumann acquired a big,rather wild but gentle Irish Setter puppy,which she called inverse in honor of our work!— HERMAN H. GOLDSTINE, The Computer.’From Pascal to von Neumann (1972)The most computationally intensive portionof the tasks assigned to the processors isintegrating the KKR matrix inverse over the first Brillouin zone.To evaluate the integral,hundreds or possibly thousands of complex double precision matricesof order between 80 and 300 must be formed and inverted.Each matrix corresponds to a different vertex of the tetrahedronsinto which the Brillouin zone has been subdivided.— M.
T. HEATH, G. A. GEIST, and J. B. DRAKE, Superconductivityin Early Experience with the Intel iPSC/860at Oak Ridge National Laboratory (1990)Pressto invert the matrix.Note that matrix inversion can produce erroneous resultsif you are using iii-conditioned matrices.— HEWLETT-PACKARD, HP 48G Series User’s Guide (1993)Almost anything you can do with A–1 can be done without it.— GEORGE E. FORSYTHE and CLEVE B. MOLER,Computer Solution of Linear Algebraic Systems (1967)261NextMATRIX INVERSION26213.1. Use and Abuse of the Matrix InverseTo most numerical analysts, matrix inversion is a sin. Forsythe, Malcolm,and Moler put it well when they say [395, 1977, p. 31] “In the vast majorit y of practical computational problems, it is unnecessary and inadvisable toactually compute A– 1 .” The best example of a problem in which the matrix inverse should not be computed is the linear equations problem Ax = b.Computing the solution as x = A-1 × b requires 2n 3 flops, assuming A - 1is computed by Gaussian elimination with partial pivoting (GEPP), whereasGEPP applied directly to the system costs only 2n3/3 flops.Not only is the inversion approach three times more expensive, but it ismuch less stable.
Suppose X = A–1 is formed exactly, and that the onlyrounding errors are in forming x = fl(Xb). Then= (X + ∆X) b , where|∆X| < γn|X|, by (3.10). So= A(X + ∆X)b = (I + A∆X)b, and the bestpossible residual bound isFor GEPP, Theorem 9.4 yieldsSince it is usually true thatfor GEPP, we see that thematrix inversion approach is likely to give a much larger residual than GEPPif A is ill conditioned and ifFor example, we solved 5025×25 systems Ax = b in MATLAB, where the elements of x are taken from thenormal N(0, 1) distribution and A is random with κ 2 (A) = u - 1 / 29 × 107.As Table 13.1 shows, the inversion approach provided much larger backwarderrors than GEPP in this experiment.Given the inexpedience of matrix inversion, why devote a chapter to it?The answer is twofold.
First, there are situations in which a matrix inversemust be computed. Examples are in statistics [54, 1974, §7.5], [721, 198 4,§2.3], [744, 1989, p. 342 ff], where the inverse can convey important statisticalinformation, in certain matrix iterations arising in eigenvalue-related problems[37, 1993], [174, 198 7], [566, 1994], and in numerical integrations arising inTable 13.1. Backward errors η A , bminx=A-1× bGEPPfor themax6.66e-12 1.69e-103.44e-18 7.56e-17-norm.13.1 U SEANDABUSEOF THEMATRIX INVERSE263superconductivity computations [509, 1990] (see the quotation at the start ofthe chapter). Second, methods for matrix inversion display a wide varietyof stability properties, making for instructive and challenging error analysis.(Indeed, the first major rounding error analysis to be published, that of vonNeumann and Goldstine, was for matrix inversion-see §9.6).Matrix inversion can be done in many different ways—in fact, there aremore computationally distinct possibilities than for any other basic matrixcomputation.
For example, in triangular matrix inversion different loop orderings are possible and either triangular matrix–vector multiplication, solutionof a triangular system, or a rank-1 update of a rectangular matrix can be employed inside the outer loop. More generally, given a factorization PA = LU,two ways to evaluate A–1 are.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
