Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 79
Текст из файла (страница 79)
To gain insight, suppose the points aredistinct and consider the case n = 3. We haveV ANDERMONDE S YSTEMS436(21.22)There is no subtractive cancellation in this product as long as each matrixhas the alternating (checkerboard) sign pattern defined, for A = (aij ), by(–l)i+jaij > 0.
This sign pattern holds for the matrices Li if the pointsα i are arranged in increasing order. The matrices Ui have the required signpattern provided that (in general)θi>0,γi>0foralli,andβ i - αk < 0 for all i + k < n – 1.In view of Table 21.2 we have the following result.Corollary 21.5. If 0 < α0 < α1 < · · · < α n then for the monomials, or theChebyshev, Legendre, or Hermite polynomials,Corollary 21.5 explains the high accuracy observed by Björck and Pereyra.Note that if|P - T ||f| < t n |P - T f | = t n |a|then, under the conditions of the corollary,which showsthat the componentwise relative error is bounded by c(n, u)tn. For the problem of Björck and Pereyra it can be shown that tnn 4/24.
Another factorcontributing to the high accuracy in this problem is that many of the subtractions αj – αj-k-1 = 1/(j + 3) – 1/(j – k + 2) are performed exactly, inview of Theorem 2.5.Note that under the conditions of the corollary P-T has the alternatingsign pattern, since each of its factors does. Thus if (–1)i fi > 0 then tn =1, and the corollary implies thatis accurate essentially to full machineprecision, independent of the condition numberIn particular, takingf to be a column of the identity matrix shows that we can compute the inverseof P to high relative accuracy, independent of its condition number.21.3 S TABILITY43721.3.2.
ResidualNext we look at the residual, r = f –Rewriting (21.21),(21.23)From the proof of Theorem 21.4 and the relation (5.9) it is easy to show that( L k + ∆Lk)-l = L k-l + Ek,|E k | < [(1 - u)-3 - 1]|Lk-1 |Strictly, an analogous bound for (Uk + ∆Uk)-l does not hold, since ∆U kcannot be expressed in the form of a diagonal matrix times U k. However,it seems reasonable to make a simplifying assumption that such a bound isvalid, say(U k + ∆Uk)-l = Uk-l + Fk,|F k | < [(1 - u)-4 - 1]|Uk-1|.(21.24)Then, writing (21.23) asand invoking Lemma 3.7, we obtain the following result.Theorem 21.6. Under the assumption (21.24), the residual of the computedsolution from Algorithm 21.2 is bounded by(21.25)where d(n, u) := (1 – u) – 7n – 1 = 7nu + O(u2).For the monomials, with distinct, nonnegative points arranged in increasing order, the matrices Li and U i are bidiagonal with the alternating signproperty, as we saw above.
Thus Li- 1 > 0 and U i -1 > 0, and since PT =we obtain from (21.25) the following pleasing result, which guarantees a tiny componentwise relative backward error.Corollary 21.7. Let 0 < α0 < α1 < · · · < αn, and consider Algorithm 21.2for the monomials. Under the assumption (21.24), the computed solutionsatisjiesV ANDERMONDE S YSTEMS438Table 21.3. Results for dual Chebyshev-Vandermonde-like system.n1020-14030-72.5 x 10 26.3 X 106.0 X 10-131.1 x 10-74.7 X 10-25.3 x 10-31.8 X 1038.3 X 10-221.3.3.
Dealing with InstabilityThe potential instability of Algorithm 21.2 is illustrated by the following example. Take the Chebyshev polynomials Ti with the points αi =(the extrema of Tn ), and define the right-hand side by fi = (–1) i . The exactsolution to PTa = f is the last column of the identity matrix.
Relative errors-16and residuals are shown in Table 21.3 (u = 10 ). Even though k2 (P) < 2for all n (see Problem 21.7), the forward errors and relative residuals are largeand grow with n. The reason for the instability is that there are large intermediate numbers in the algorithm (at the start of stage II for n = 40,is of order 1015); hence severe cancellation is necessary to produce the finalanswer of order 1.
Looked at another way, the factorization of PT used bythe algorithm is unstable because it is very sensitive to perturbations in thefactors.How can we overcome this instability? There are two possibilites: prevention and cure. The only means at our disposal for preventing instability is toreorder the points αi . The ordering is arbitrary subject to condition (21.5)being satisfied.
Recall that the algorithms construct an LU factorization ofPT in factored form, and note that permuting the rows of PT is equivalentto reordering the points αi . A reasonable approach is therefore to take whatever ordering of the points would be produced by Gaussian elimination withTpartial pivoting (GEPP) applied to P . Since the diagonal elements of U inTP = LU have the formwhere ha depends only on the θi , and since partial pivoting maximizes thepivot at each stage, this ordering of the α i can be computed at the start of thealgorithm in n2 flops (see Problem 21.8).
This ordering is essentially the Lejaordering (5.13) (the difference is that partial pivoting leaves α0 unchanged).To attempt to cure observed instability we can use iterative refinementin fixed precision. Ordinarily, residual computation for linear equations is21.3 S TABILITY439trivial, but in this context the coefficient matrix is not given explicitly andcomputing the residual turns out to be conceptually almost as difficult, andcomputationally as expensive, as solving the linear system!To compute the residual for the dual system we need a means for evaluating ψ(t) in (21.7) and its first k < n derivatives, where k = max i (i – r(i))is the order of confluence.
Since the polynomials pj satisfy a three-term recurrence relation we can use an extension of the Clenshaw recurrence formula(which evaluates ψ but not its derivatives). The following algorithm implements the appropriate recurrences, which are given by Smith [927, 1965] (seeProblem 21.9).Algorithm 21.8 (extended Clenshaw recurrence). This algorithm computesthe k + 1 values yj = ψ (j)(x), 0 < j < k, where ψ is given by (21.7) and k < n.It uses a work vector z of order k.Set y(0: k) = z(0: k) = 0y0 = anfor j = n – 1: – 1:0temp = y0y0 = θj (x – β j )z0 = tempfor i = 1: min(k, n – j)temp = yiyi = θj((x – β j)yi + zi-1) – γj + l z izi = tempendendm=1for i = 2:km=m*iyi = m * yiendComputing the residual using Algorithm 21.8 costs between 3n2 flops (forfull confluence) and 6n2 flops (for the nonconfluent case); recall that Algorithm 21.2 costs at most 9n 2 /2 flops!The use of iterative refinement can be justified with the aid of Theorem 11.3.
For (confluent) Vandermonde matrices, the residuals are formedusing Homer’s rule and (11 .7) holds in view of (5.3) and (5.7). Hence forstandard Vandermonde matrices, Theorem 11.3 leads to an asymptotic componentwise backward stability result. A complete error analysis of Algorithm 21.8 is not available for (confluent) Vandermonde-like matrices, butit is easy to see that (11.7) will not always hold. Nevertheless it is clear that440V ANDERMONDE S YSTEMSa normwise bound can be obtained (see Oliver [805, 1977] for the special caseof the Chebyshev polynomials) and hence an asymptotic normwise stabilityresult can be deduced from Theorem 11.3. Thus there is theoretical backingfor the use of iterative refinement with Algorithm 21.8.Numerical experiments using Algorithm 21.8 in conjunction with the partial pivoting reordering and fixed precision iterative refinement show that bothtechniques are effective means for stabilizing the algorithms, but that iterativerefinement is likely to fail once the instability y is sufficiently severe.
Becauseof its lower cost, the reordering approach is preferable.Two heuristics are worth noting. Consider a (confluent) Vandermondelike system Px = b Heuristic 1: it is systems with a large normed solutionthat are solved to high accuracy by the fast algorithms.To produce a large solution the algorithms must sustain little cancellation,and the error analysis shows that cancellation is the main cause of instability.Heuristic 2: GEPP is unable to solve accurately Vandermonde systems witha very large normed solutionThe pivots for GEPPwill tend to satisfyso that the computed solution will tendto satisfyA consequence of these two heuristics is thatfor Vandermonde(-like) systems with a very large-normed solution the fastalgorithms will be much more accurate (but no more backward stable) thanGEPP.
However, we should be suspicious of any framework in which such systems arise; although the solution vector may be obtained accurately (barringoverflow), subsequent computations with numbers of such a wide dynamicrange will probably themselves be unstable.21.4. Notes and ReferencesThe formulae (21. 1) and (21.2), and inversion methods based on these formulae, have been discovered independently by many authors. Traub [1013, 1966,14] gives a short historical survey, his earliest reference being a 1932 bookby Kowalewski. There does not appear to be any published error analysisfor Algorithm 21.1 (see Problem 21.3).
There is little justification for usingthe output of the algorithm to solve the primal or dual linear system, as isdone in [842, 1992, §2.8]; Algorithms 21.2 and 21.3 are more efficient and almost certainly at least as stable. Calvetti and Reichel [179, 1993] generalizeAlgorithm 21.1 to Vandermonde-like matrices, but they do not present anyerror analysis. Gohberg and Olshevsky [456, 1994] give another O(n2) flopsalgorithm for inverting a Chebyshev–Vandermonde matrix.The standard condition number K(V) is not an appropriate measure ofsensitivity when only the points αi are perturbed, because it does not reflectthe special structure of the perturbations. Appropriate condition numberswere first derived by Higham [533, 1987] and are comprehensively investigatedPROBLEMS441by Bartels and D.
J. Higham [76, 1992]; see Problem 21.10.Methods for solving the dual and primal Vandermonde systems have an interesting history. The earliest algorithm was derived by Lyness and Moler [718,1966] via Neville interpolation; it solves the dual system in O(n3) flops. Thefirst O(n2 ) algorithm was obtained by Ballester and Pereyra [52, 1967]; itcomputes the LU factors of the Vandermonde matrix and requires O(n2) elements of storage. Björck and Pereyra [121, 1970] derived the specializationof Algorithms 21.2 and 21.3 to nonconfluent Vandermonde matrices; thesealgorithms require no storage other than that for the problem data. The algorithms of Björck and Pereyra were generalized by Björck and Elfving toconfluent systems [117, 1973], and by Higham to Vandermonde-like systems[536, 1988] and confluent Vandermonde-like systems [547, 1990].
The erroranalysis in this chapter is taken from [547, 1990]. Tang and Golub [992, 1981]give a block algorithm that requires only real arithmetic to solve a Vandermonde system in which all the points appear in complex conjugate pairs.Other O(n2 ) algorithms for solving Chebyshev-Vandermonde systems aregiven by Reichel and Opfer [865, 1991] and Calvetti and Reichel [178, 1992].The former algorithms are progressive, in that they allow the solution to beupdated when a new point αi is added; they generalize progressive algorithmsof Björck and Pereyra [121, 1970]. Boros, Kailath, and Olshevsky [135, 1994]use the concept of displacement structure to derive further O(n2) algorithmsfor solving Vandermonde and Chebyshev–Vandermonde systems.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
