Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 33
Текст из файла (страница 33)
Suppose the upper triangular matrixsatisfies(8.4)Then the unit upper triangular matrix W = |U- 1 ||U| satisfies wij < 2 j-i forall j > i.T RIANGLLAR S YSTEMS156Proof. We can write W = |V- 1 ||V| where V = D -1 U and D = diag(u ii ).The matrix V is unit upper triangular with) |υij| < 1, and it is easy to showthat |(V - 1 ) ji| < 2 j - i -1 for j > i. Thus.
for j > i,Theorem 8.7. Under the conditions of Lemma 8.6, the computed solutionto Ux = b obtained by substitution satisfiesProof. From Theorem 8.5 we haveUsing Lemma 8.6 we obtainLemma 8.6 shows that for matrices satisfying (8.4), cond(T) is boundedfor fixed n, no matter how large κ(T). The bounds forin Theorem 8.7,although large if n is large and i is small. decay exponentially with increasingthus, later components of x are always computed to high accuracy relativeto the elements already computed.Analogues of Lemma 8.6 and Theorem 8.7 hold for lower triangular Lsatisfying|lii | > |lij| for all j < i.(8.5)Note, however, that if the upper triangular matrix T satisfies (8.4) then TTdoes not necessarily satisfy (8.5).
In fact, cond(T T) can he arbitrarily large,as shown by the exampleAn important conclusion is that a triangular system Tx = b can be muchmore or less ill conditioned than the system TTy = c, even if T satisfies (8.4).Theorem 8.7, or its lower triangular analogue, is applicable to1578.2 F ORWARD E RROR A NALYSIS• the lower triangular matrices from Gaussian elimination with partialpivoting or complete pivoting;• the upper triangular matrices from Gaussian elimination with completepivoting;• the upper triangular matrices from the Cholesky and QR factorizationswith complete pivoting and column pivoting, respectively.Next, we consider triangular T satisfyingtii > 0, tij < 0 for all ij.It is easy to see that such a matrix has an inverse with nonnegative elements,and hence is an hi-matrix (for definitions of an hi-matrix see Appendix B).Associated with any square matrix A is the comparison matrix:(8.6)For any nonsingular triangular T, M(T) is an M-matrix.
Furthermore, it iseasy to show that |T -1| < M(T)-1 (see Theorem 8.11).The following result shows that among all matrices R such that |R| = |T|,R = M(T) is the one that maximizes cond(R , x ).Lemma 8.8. For any triangular T,Proof. The inequality follows from |T- 1 | < M(T)-1, together with |T| =|M(T)|. Since M(T)-1 > 0, we have|M(T)- 1 ||M(T)| = M(T) -1 (2diag(|tii| ) - M(T))= 2M(T) -1diag(|tii |) - I,which yields the equality.If T=M(T) has unit diagonal then, using Lemma 8.8,This means, for example, that the system U(1)x = b (see (8.2)), where x = e,is about as ill conditioned with respect to componentwise relative perturbations in U(1) as it is with respect to normwise perturbations in U(1).158T RIANGULAR S YSTEMSThe next result gives a forward error bound for substitution that is proveddirectly, without reference to the backward error result Theorem 8.5 (indeed.
itcannot be obtained from that result!). The bound can be very weak, because||M(T) - 1 || can be arbitrarily larger than ||T- 1|| (see Problem 8.2), but ityields a pleasing bound in the special case described in the corollary.Theorem 8.9. The computed solutionobtained from substitution appliedto the triangular system Tx = b of order n satisfiesProof. Without loss of generality, suppose T = L, is lower triangular. Theproof is by induction on the components of x. The result clearly holds forthe first component.
Assume that it holds for the first n - 1 components. Ananalogue of Lemma 8.3 shows thatwheregivesfor all j. Subtracting from lnnxn =so that(8.7)WriteThen the inductive assumption can be written aswhich impliesHence (8.7) gives8.3 BOUNDSFOR THEINVERSE159Corollary 8.10. The computed solution obtained from substitution appliedto the triangular system Tx = b of order n, where T = M(T) and b > 0,satisfiesCorollary 8.10 shows that, when T is an M-matrix and the right-handside is nonnegative, the solution is obtained to high relative accuracy in everycomponent. The reason for the high accuracy is that for such a system thereare no subtractions of like-signed numbers, so that each xi is computed asa sum of nonnegative quantities.
A consequence of the corollary is that theinverse of a triangular M-matrix can be computed to high relative accuracy.Triangular systems of the type in Corollary 8.10 occur in linear equationsobtained by discretizing certain elliptic partial differential equations, such asthe Poisson equation on a rectangle, with zero boundary conditions and apositive forcing function: these problems yield symmetric positive definiteM-matrices, and the LU factors of an M-matrix are themselves M-matrices.Such systems also occur when evaluating the bounds of the next section.8.3.
Bounds for the InverseIn this section we describe bounds for the inverse of a triangular matrix andshow how they can be used to bound condition numbers. All the bounds inthis section have the property that they depend only on the absolute valuesof the elements of the matrix. The norm estimation methods of Chapter 14,on the other hand, do take account of the signs of the elements.The bounds are all based on the following theorem, whose easy proof weomit.Theorem 8.11. If U is a nonsingular upper triangular matrix then|U- 1 | < M(U)-1 < W(U)-1 < Z(U)- 1 ,where the upper triangular matrices W(U) and Z(U) are defined as follows:160T RIANGULAR S YSTEMSwhere a = mink |ukk|, β = maxi<j|uij|/|uii|.Theorem 8.11 is a special case of results in the theory of Ill-matrices.For more general results couched in terms of matrix minorants and diagonaldominance, respectively, see Dahlquist [261, 1983] and Varga [1051, 1976]; seealso Householder [587, 1964, Exercise 15, p.
58].An obvious implication of the theorem is that for any vector z and anyabsolute normBy taking z = |U|e, z = |U||x|, and z = e, respectively, we obtain upperThe cost of computing thesebounds for cond(U), cond(U,x), andbounds is just the cost of solving a triangular system with coefficient matrixM(U), W(U), or Z(U), which is easily seen to be O(n2). O(n), and O(1)flops, respectively. By comparison, computing any of these condition numbersexactly costs O(n3) flops.As an example, here is how to compute an upper bound forin n 2flops.Algorithm 8.12. Given a nonsingular upper triangular matrixthis algorithm computes µ =y n = 1 /|unn|for i = n - 1:-1:1s=1for j = i + 1:ns = s + |uij|yjendy i = y i /|uii|endµ=How good are these upper bounds? We know from Problem 8.2 that theratio ||M(T)- 1 ||/||T - 1 || can be arbitrarily large, therefore any of the upperbounds can be arbitrarily poor.
However, with suitable assumptions on T,more can be said.It is easy to show that if T is bidiagonal then |T- 1 | = M(T)- 1 . Sincea bidiagonal system can be solved in O(n) flops, it follows that the threecondition numbers of interest can each be computed exactly in O(n) flopswhen T is bidiagonal.As in the previous section, triangular matrices that result from a pivotingstrategy also lead to a special result.8.3 BOUNDSFOR THEINVERSE161Theorem 8.13.
Suppose the upper triangular matrixsatisfies|uii | > |uij | for all j > i.Then, for the 1-, 2-, and-norms,(8.8)Proof. The left-hand inequality is trivial. The right-hand inequalityfollows from the expression(see Problem 8.5),together with ||A||2 <The inequalities from the second on in (8.8) are all equalities for the matrixwith u ii = 1 and u ij = -1 (j > i). The question arises of whether equality ispossible for the upper triangular matrices arising from QR factorization withcolumn pivoting, which satisfy the inequalities (see Problem 18.5)(8.9)That equality is possible is shown by the parametrized matrix of Kahan [626,19 66](8.10)where c = cos(θ), s = sin(θ). It is easily verified that Un (θ) satisfies theinequalities (8.9)-as equalities, in fact. Prom (8.3), Un ( θ )-1 = (βi j ) is givenbyThus asand hence, for small θ,whereIt can also be verified that the matrixdefined by(-1) j-i|uij| satisfies, for small θ,while2 n - 1 /|unn|. Hence the upper bounds for ||U - 1 || can all be too big by a factorof order 2n - 1 .162T RIANGULAR S YSTEMS8.4.
A Parallel Fan-In AlgorithmSubstitution is not the only way to solve a triangular system. In this section wedescribe a different approach that has been suggested for parallel computation.Any lower triangular matrixcan be factorized L = L1 L2 . . . Ln ,where Lk differs from the identity matrix only in the kth column:(8.11)The solution to a linear system Lx = b may therefore be expressed as(8.12)where M i =When evaluated in the natural right-to-left order, thisformula yields a trivial variation of a column-oriented version of substitution.The fan-in algorithm evaluates the product (8.12) in [log(n + 1)] steps bythe fan-in operation (which is the operation used in pairwise summation: see§4.1).
For example, for n = 7 the calculation is specified bywhere all the products appearing within a particular size of parenthesis canbe evaluated in parallel. III general. the evaluation can be expressed as abinary tree of depth [log(n + 1)] + 1, with products M1b and Mi Mi-1 (i =3, 5,. . . , 2[(n - 1)/2] + 1) at the top level and a single product yielding x atthe bottom level.
This algorithm was proposed and analysed by Sameh andBrent [889, 1977], who show that it can be implemented intime steps onprocessors. The algorithm requires about n 3 /10operations and thus is of no interest for serial computation. Some pertinentcomments on the practical significance of log n terms in complexity results aregiven by Edelman [341, 1993].To derive an error bound while avoiding complicated notation that obscures the simplicity of the analysis, we take n = 7.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
