Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 34
Текст из файла (страница 34)
The result we obtain iseasily seen to be valid for all n. We will not be concerned with the precisevalues of constants, so we write cn for a constant depending on n . We assumethat the inverses Mi =are formed exactly, because the errors in formingthem affect only the constants. From the error analysis of matrix-vector andmatrix-matrix multiplication (§3.5), we find that the computed solutionsatisfies(8.13)8.4 A P ARALLEL FAN- I N A LGORITHM163wherePremultiplying (8.13) on the left by L, we find that the residual r =is a sum of terms of the formAll these terms share the same upper bound, which we derive for just one ofthem. For j = 5, k =where we have used the property that, for anyThe overall residual bound is therefore of the form(8.14)or, on taking norms,(8.15)By considering the binary tree associated with the fan-in algorithm, andusing the fact that the matrices at the ith level of the tree have at most 2i - 1nontrivial columns, it is easy to see that we can take dn = an log n , where ais a constant of order 1.It is not hard to find numerical examples where the bound in (8.15) isapproximately attained (for dn = 1) and greatly exceedswhichis the magnitude required for normwise backward stability.
One way to construct such examples is to use direct search (see Chapter 24).The key fact revealed by (8.15) is that the fan-in algorithm is only conditionally stable. In particular, the algorithm is normwise backward stable if Lis well conditioned. A special case in which (8.14) simplifies is when L is an M matrix and b > 0: Problem 8.4 shows that in this case |L - 1 ||L||x| < (2n-1)|x|,so (8.14) yields< (2n - 1)2 dnu|L||x| + O(u 2 ), and we have componentwise backward stability (to first order).164T RIANGULAR S YSTEMSWe can obtain from (8.15) the result(8.16)which was proved by Sameh and Brent (889, 1977] (with a n = ¼n 2 log n +O(n log n)). However, (8.16) is a much weaker bound than (8.14) and (8.15).In particular, a diagonal scaling Lx = b(where Dj isdiagonal) leaves (8.14) (and, to a lesser extent, (8.15)) essentially unchanged,but can change the bound (8.16) by an arbitrary amount.A forward error bound can be obtained directly from (8.13).
We find that(8.17)where M(L) is the comparison matrix (a bound of the same form as thatin Theorem 8.9 for substitution-see the Notes and References and Problem 8.10). which can be weakened to(8.18)We also have the bound(8.19)which is an immediate consequence of (8.14). Either bound in (8.18) and(8.19) can be arbitrarily larger than the other, for fixed n .
An example where(8.19) is the better bound (for large n) is provided by the matrix with lij ,1,for which |L- 1 ||L| has maximum element 2 and M(L)- 1 |L| has maximumelement 2n - 1 .8.5. Notes and ReferencesSection 8.2 is based on Higham [538, 198 9]. Many of the results presentedin §§8.2 and 8.3 have their origin in the work of Wilkinson. Indeed, thesesections are effectively a unification and extension of Wilkinson’s results in[1085, 1961], [1088, 1963], [1089, 1965].Classic references for Theorems 8.3 and 8.5 are Wilkinson [1085, 196 1,p.
294], [1088, 1963, pp. 100-102], Forsythe and Moler [396, 1967, §21], andStewart [941, 1973, pp. 150, 408-410].Analogues of Theorem 8.7 and Corollary 8.10 for matrix inversion areproved by Wilkinson in [1085, 1961, pp. 322 323], and Corollary 8.10 itself isproved in [1089, 1965, pp. 250-251].8.5 N OTESANDR EFERENCES165A result of the form of Theorem 8.9 holds for any triangular system solverthat does not rely on algebraic cancellation-in particular, for the fan-in algorithm, as already seen in (8.17).
See Problem 8.10 for a more precise formulation of this general result.The bounds in §8.3 have been investigated by various authors. The unifiedpresentation given here is based on Higham [534, 1987]. Karasalo [642, 1974]derives an O(n) flops algorithm for computing ||M(T) - 1 ||F. Manteuffel [726,1981] derives the first two inequalities in Theorem 8.11, and Algorithm 8.12.A different derivation of the equations in Algorithm 8.12 is given by Jennings [613, 1982, §9]. The formulae given in Problem 8.5 are derived directlyas upper bounds forby Lemeire [699, 1975].Thatcan be computed in O(n) flops when B is bidiagonal, aswas first pointed out by Higham [531, 1986]. Demmel andKahan [296, 1990] derive an estimate for the smallest singular value σmin ofa bidiagonal matrix B by using the inequalitywhereThey computein O(n) flops asSection 8.4 is adapted from Higham [560, 1995], in which error analysis isgiven for several parallel methods for solving triangular systems.The fan-in method is topical because the fan-in operation is a special caseof the parallel prefix operation and several fundamental computations in linearalgebra are amenable to a parallel prefix-based implementation, as discussedby Demmel [287, 1992], [288, 1993].
(For a particularly clear explanation of theparallel prefix operation see the textbook by Buchanan and Turner [154, 1992,§13.21.) The important open question of the stability of the parallel prefiximplementation of Sturm sequence evaluation for the symmetric tridiagonaleigenproblem has recently been answered by Mathias [734, 1995]. Mathiasshows that for positive definite matrices the relative error in a computed minorcan be as large as a multiple ofwhereis the smallest eigenvalue ofthe matrix; the corresponding bound for serial evaluation involvesTheanalogy with (8.19), where we also see a condition cubing effect, is intriguing.Higham and Pothen [568, 1994] analyse the stability of the “partitionedinverse method” for parallel solution of sparse triangular systems with manyright-hand sides.
This method has been studied by several authors in the1990s; see Alvarado, Pothen, and Schreiber [13, 1993] and the referencestherein. The idea of the method is to factor a sparse triangular matrixas L = L 1 L 2 . . . Ln = G 1 G 2 . . . Gm , where each Gi is a product of consecutive Lj terms and 1 < m < n, with m as small as possiblesubject to the Gi being sparse. Then the solution to Lx = b is evaluated asT RIANGULAR S YSTEMS166where eachis formed explicitly and the product is evaluated from rightto left. The advantage of this approach is that x can be computed in m serialsteps of parallel matrix-vector multiplication.8.5.1. LAPACKComputational routine xTRTRS solves a triangular system with multiple righthand sides: xTBTRS is an analogue for banded triangular matrices.
There isno driver routine for triangular systems.ProblemsBefore you start an exercise session,make sure you have a g/ass of water anda mat or towel nearby.-MARIE HELVIN, Mode/ Tips for a Healthy Future (1994)8.1. Show that under the no-guard-digit model (2.6). Lemma 8.2 remainstrue if (8.1) is changed toand that the corresponding modification of Theorem 8.5 has8.2. Show that for a triangular matrix T the ratio ||M(T) - 1||/||T - 1 || can bearbitrarily large.8.3. Suppose the unit upper triangular matrixsatisfies |uij| < 1for j > i.
By using Theorem 8.9. show that the computed solutionfromsubstitution on Ux = b satisfiesCompare with the result of applying Theorem 8.7.8.4. Letbe triangular and suppose T = M(T) and Tx = b > 0.Show that |T- 1||T||x| < (2n-1)|x|, and hence that cond(T, x) < 2n-1. Thisresult show that a triangular M-matrix system with a nonnegative righthand side is very well conditioned with respect to componentwise relativeperturbations, irrespective of the size of κ(T) (and so leads to an alternativeproof of Corollary 8.10).8.5.
Show that for a triangularboth the l- and -norms (a and β are defined in Theorem 8.11).forP ROBLEMS1678.6. Write detailed algorithms for efficiently computingand8.7. Bounds from diagonal dominance. (a) Prove the following result (Ahlbergand Nilson [8, 1963], Varah [1049, 1975]): if(not necessarily triangular) satisfies(that is, A is strictly diagonally dominant by rows), then(b) Hence show that (Varga [1051,1976])ifsatisfiesfor some positive diagonal matrix D = diag(di ) (that is, AD is strictly diagonally dominant by rows), then(c) Use part (b) to provide another derivation of the upper bound8.8. (a) Letbe nonsingular.
For a given i and j, determine, ifpossible, a ij such that A +is singular. Where is the “best” place toperturb A to make it singular?(b) Let T = U(1) in (8.2), so that, for example,Show that Tn is made singular by subtracting 22-n from a certain element ofTn.8.9. (Zha [1127, 1993]) Show that if c and s are nonnegative (with c2 + s2 = 1)then the Kahan matrix Un (θ) in (8.10) hasas its second smallestsingular value. (That there should be such an explicit formula is surprising;none is known for the smallest singular value.)168TRIANGULAR SYSTEMS8.10. Consider a method for solving triangular systems Tx = b that computesxi = fi (T, b) where, for all i, fi is a multivariate rational function in which theonly divisions are by diagonal elements of L and such that when T = M(T)and b > 0 there are no subtractions in the evaluation of fi .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
