Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 36

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 36)

Applying the lemma, wededuce that, no matter what the-ordering of the inner products, the computedmatrices L and U satisfy (with lkk := 1)These inequalities constitute a backward error result for LU factorization.Theorem 9.3. If GE applied tothe computed LU factorsand(m > n) runs to completion thensatisfy(9.6)With only a little more effort, a backward error result can be obtained forthe solution of Ax = b.Theorem 9.4. Letand suppose GE produces computed LU factorsand a computed solutionto Ax = b.

Then(9.7)Proof. From Theorem 9.3,Theorem 8.5, substitution producesByandsatisfyingThuswhereWe need a constant 2γ n instead ofAlthough it is not usually worth expending effort reducing constantsin error bounds (see the Wilkinson quotation at the start of Chapter 10), wewill optimize constants in this important case. Consideration of Lemma 8.4shows that we actually have176LU FACTORIZATIONANDLINEAR E QUATIONSso thatusing Lemma 3.3, which gives the required constant.How do we interpret Theorem 9.4? Ideally, we would like |∆ A| < u|A|,which corresponds to the uncertainty introduced by rounding the elements ofA, but because each element of A undergoes up to n arithmetic operations wecannot expect better than a bound |∆ A| < cnu|A|, where cn is a constant oforder n.

Such a bound holds ifand satisfywhich certainlyholds ifandare nonnegative. because then (9.6) gives(9.8)Substituting into (9.7). we obtainThis result says thathas a small componentwise relative backward error.One class of matrices that has nonnegative LU factors is defined as follows.is totally positive (nonnegative) if the determinant of every squaresubmatrix is positive (nonnegative). In particular, this definition requires thata ij and det(A) be positive or nonnegative. Some examples of totally nonnegative matrices are given in Chapter 26. If A is totally nonnegative then it hasan LU factorization A = LU in which L and U are totally nonnegative. sothat L > 0 and U > 0 (see Problem 9.6); moreover, the computed factorsandare nonnegative for sufficiently small values of the unit roundoff u [273,1977 ].

Inverses of totally nonnegative matrices also have the property that|A| = |L||U| (see Problem 9.7). Note that the property of a matrix or itsinverse being totally nonnegative is generally destroyed under row permutations. Hence for totally- nonnegative matrices and their inverses it is best touse Gaussian elimination without pivoting.One important fact that follows from (9.6) and (9.7) is that the stabilityof GE is determined not by the size of the multipliers but by the size of thematrixThis matrix can be small when the multipliersare large, andlarge when the multipliers are of order 1 (as we will see in the next section).To understand the stability of GE further we turn to norms. For GE without pivoting. the ratio || |L||U|||/||A|| can be arbitrarily large.

For example,for the matrixthe ratio is of orderAssume then that partial pivoting is used. Then |lij | < 1 for all i > j, since the lij are the multipliers.9.3 T HE G ROWTH FACTOR177And it is easy to show by induction that |uij| < 2 i -1 maxk < i |akj|. Hence, forpartial pivoting, L is small and U is bounded relative to A.Traditionally, backward error analysis for GE is expressed in terms of thegrowth factorwhich involves all the elementstion. Using the boundclassic theorem.Theorem 9.5 (Wilkinson). Letoting produces a computed solution(k = 1:n) that occur during the eliminap n maxi,j|aij| we obtain the followingand suppose GE with partial pivto Ax = b.

Then(9.9)We hasten to admit to using an illicit manoeuvre in the derivation of thistheorem: we have used bounds for L and U that strictly are valid only for theexact L and U. We could instead have defined the growth factor in terms ofthe computedbut then any bounds for the growth factor would involvethe unit roundoff (similarly, we can only guarantee that< 1 + u). Ourbreach of correctness is harmless for the purposes to which we will put thetheorem.The assumption in Theorem 9.5 that partial pivoting is used is not necessary: essentially the same result holds for GE without pivoting (see Problem 9.8). The normwise backward stability of GE with or without pivoting istherefore governed by the growth factor, to which we now turn our attention.9.3.

The Growth FactorIt is easy to show that pn < 2 n -1 for partial pivoting. Wilkinson notes thatthis upper bound is achieved for matrices of the form illustrated for n = 4 byFor these matrices, no interchanges are required by partial pivoting, and thereis exponential growth of elements in the last column of the reduced matrices.In fact, this is just one of a nontrivial class of matrices for which partial178LU FACTORIZATIONANDLINEAR E QUATIONSpivoting achieves maximal growth. When necessary in the rest of this chapter,we denote the growth factor for partial pivoting byand that for completepivoting byTheorem 9.6 (Higham and Higham).

All real n × n matrices A for which2n -1 are of the formwhere D = diag(±1). M is unit lower triangular with mij = -1 for i > j,T is an arbitrary nonsingular upper triangular matrix of order n -1, d =(1,2,4,,.... 2n -1)T, and a is a scalar such that a := |a1n| = maxi,j |aij|.Proof. GE with partial pivoting applied to a matrix A gives a factorizationB := PA = LU, where P is a permutation matrix. It is easy to show that|uij | < 2i -1 maxr< i |brj|, with equality for i = s only if there is equality fori = 1:s - 1.

Thus pn = 2n -1 implies that the last column of U has the formaDd, and also that |b1 n| = maxi,j |bij|. By considering the final column of B,and imposing the requirement that |lij| < 1, it is easy to show that the unitlower triangular matrix L must have the form L = DMD. It follows thatat each stage of the reduction every multiplier is ±1; hence no interchangesare performed, that is, P = I.

The only requirement on T is that it benonsingular, for if t i i = 0 then the ith elimination stage would be skippedbecause of a zero pivot column and no growth would be produced on thatstage.Note that by varying the elements mij (i > j) and the vector d in Theorem 9.6 we can construct matrices for whichachieves any desired valuebetween 1 and 2n - 1 .Despite the exist once of matrices for which pn is large with partial pivoting, the growth factor is almost invariably small in practice. For example,Wilkinson says “It is our experience that any substantial increase in size ofelements of successive Ar is extremely uncommon even with partial pivotingNo example which has arisen naturally has in my experience given anincrease by a factor as large as 16” [1089, 1965, pp.

213-214].Until recently. there were no reports in the literature of large growth factorsbeing observed in practical applications. However, Wright [1116, 1993] hasfound a class of two-point boundary value problems that. when solved by themultiple shooting met hod, yield a linear system for which partial pivotingsuffers exponential growth.

The matrix is block lower bidiagonal, except fora nonzero block in the top right-hand corner. Furthermore, Foster [399, 1994]shows that a quadrature met hod for solving a practically occurring Volterra9.3 T HE G ROWTH FACTOR179integral equation gives rise to linear systems for which partial pivoting againgives large growth factors.There exist some well-known matrices that give unusually large, but notexponential, growth. They can be found using the following theorem, which isapplicable whatever the strategy for interchanging rows and columns in GE.Theorem 9.7 (Higham and Higham). Letbe nonsingular and seta = max i,j |aij|, β = maxi,j |(A- 1 )ij |, and θ = (aβ)-1.

Then θ < n, and forany permutation matrices P and Q such that PAQ has an LU factorization.the growth factor pn for GE without pivoting on PAQ satisfies p n > θ.Proof. The inequality θ < n follows froman LU factorization PAQ = LU computed by GE. We haveConsider(9.10)Henceand the result follows.Note that θ -1 = aβ satisfiesClearly, Ahas to be very well conditioned for the theorem to provide a lower bound θnear the maximum of n.We apply the theorem to three noncontrived matrices that appear in practical applications.(1) The matrix(9.11)is the symmetric, orthogonal eigenvector matrix for the second difference matrix (the tridiagonal matrix with typical row (-1, 2, -1)-see §26.5); it arises,for example, in the analysis of time series [19, 1971, §6.5]. Theorem 9.7 givesp n (S N) > (n + 1)/2.(2) A Hadamard matrix H n is an n × n matrix with elements h ij = ±1and for which the rows of Hn are mutually orthogonal.

Hadamard mat ricesexist only for certain n; a necessary condition for their existence if n > 2 isthat n is a multiple of 4. For more about Hadamard mat rices see Hall [494,19 6 7 , Chap. 14], Wallis [1062, 1993 ], and Wallis, Street, and Wallis [1063,1972 ]. We have= nI, and soTheorem 9.7 givespn > n.(3) The next matrix is a complex Vandermonde matrix based on the rootsof unity, which occurs in the evaluation of Fourier transforms (see §23.1):(9.12)LU FACTORIZATION180ANDLINEAR E QUATIONSSinceTheorem 9.7 gives pn (Vn ) > n.Note that each of these matrices is orthogonal or unitary (to within a rowscaling in the case of the Hadamard matrix), so it is not necessary to apply GEto them! This may explain why growth factors of order n for these matriceshave not been reported in the literature as occurring in practice.To summarize, although there are practically occurring matrices for whichpartial pivoting yields a moderately large, or even exponentially large, growthfactor, the growth factor is almost invariably found to be small.

Характеристики

Список файлов книги