Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 77
Текст из файла (страница 77)
The computedsatisfiesFrom Lemma 18.3, the computed solutionandsatisfies(20.10)We now rewrite the latter two equations in such a way that we can applyLemma 20.2:It is straightforward to show that|∆ Ai | < mnγ´cn|A|Gi ,and ||∆Ai ||F <Lemma 20.2.mγ´||A||F, i = 1:2.||Gi ||F = 1, i = 1:2,The result follows on invocation of20.3 E RROR ANALYSIS421Theorem 20.3 says that the Q method is row-wise backward stable. This isnot altogether surprising, since (Householder or Givens) QR factorization forthe LS problem enjoys an analogous backward stability result (Theorem 19.3),albeit without the restriction of a minimum norm solution.
Applying (20.8)to Theorem 20.3 we obtain the forward error bound(20.11)The same form of forward error bound (20.11) can be derived for the SNEmethod as for the Q method [292, 1993]. However, it is not possible to obtaina result analogous to Theorem 20.3, nor even to obtain a residual bound of theform(which would imply that solved a nearbysystem, though would not necessarily be the minimum norm solution). Themethod of solution guarantees only that the seminormal equations themselveshave a small residual. Thus, as in the context of overdetermined LS problems,the SNE method is not backward stable.
A possible way to improve thestability is by iterative refinement, as shown in [292, 1993].Note that the forward error bound (20. 11) is independent of the row scalingof A, since cond2 (A) is. The bound is therefore potentially much smaller thanthe boundobtained by Paige [813, 1973] for the SNE method and by Jennings and Osborne [614, 1974] and Arioli and Laratta [26, 1985, Thin. 4] for the Q method.Finally, we mention an alternative version of the Q method that is basedon the modified Gram–Schmidt (MGS) method. The obvious approach isto compute the QR factorization AT = QR using MGSsolve RTy = b, and then form x = Qy.
Since Q is provided explicitlyby the MGS method, the final stage is a full matrix–vector multiplication,unlike for the Householder method. However, because the computed Q maydepart from orthonormality, this method is unstable in the form described.The formation of x = Qy should instead be done as follows:The recurrence can be written as x(k-l) = x(k) + ykq k – (qkTx(k))qk, and thelast term is zero in exact arithmetic if the qk are mutually orthogonal.
In finiteprecision arithmetic this correction term has the “magical” effect of making422U NDERDETERMINED S YSTEMSTable 20.1. Backward errors for underdetermzned Vandermonde system.Householder QRMGS with x := QyMGS with x formed stably (see text)SNE method (using Householder QR)9.76 X 10-184.10 x 10-42.25 X 10-171.99 x 10-4the algorithm stable, in the sense that it satisfies essentially the same resultas the Q method in Theorem 20.3; see Björck and Paige [120, 1994].A numerical example is instructive. Take the 20 x 30 Vandermonde matrixwhere the pi are equally spaced on [0, 1], and let b have elementsequally spaced on [0, 1]. The condition number k2 (A) = 4.35 x 1014 . The(standard) normwise backward errors in the 2-norm are shown in Table 20.1.For A T , thesupplied by MGS satisfies= 1.41 x 10-3, whichexplains the instability of the “obvious” MGS solver.20.4.
Notes and ReferencesThe seminormal equations method was suggested by Gill and Murray [442,1973] and Saunders [894, 1972]. Other methods for obtaining minimal 2-normsolutions of underdetermined systems are surveyed by Cline and Plemmons[219, 1976].Theorem 20.1 is from Demmel and Higham [292, 1993].
The bound (20.9)is well known; it follows from Wedin’s original version of our Theorem 19.1,which applies to minimum 2-norm underdetermined problems as well as LSproblems.Theorem 20.3 is new. Demmel and Higham [292, 1993] prove the weakerresult that from the Q method is very close to a vector that satisfies thecriterion for row-wise backward stability, and Lawson and Hanson [695, 1974,Thin. 16. 18] give a corresponding result in which satisfies the criterion forgeneral normwise backward stability.
The key to showing actual backwardstability is the use ofand Schwetlick’s lemma, which is a modification of Lemma 8.2.11 in [658, 1988] and [659, 1992] (our Lemma 19.6).Demmel and Higham [292, 1993] also give error analysis for the seminormalequations method.The new MGS algorithm for solving the minimum norm problem was firstsuggested by Björck and Paige [119, 1992]; see also Björck [115, 1994].Arioli and Laratta [27, 1986] give error analysis of QR factorization methods for solving the general problem min{ ||x – c||2 : Ax = b }, where APROBLEMS423with m < n.20.4.1.
LAPACKThe same routines that solve the (overdetermined) LS problem also solve underdetermined systems for the solution of minimal 2-norm. Thus, xGELS solvesa full rank underdetermined system with multiple right-hand sides by the Qmethod. Routines xGELSX and xGELSS solve rank-deficient problems with multiple right-hand sides, using, respectively, a complete orthogonal factorization(computed via QR factorization with column pivoting) and the singular valuedecomposition.Problems20.1. (Björck [114, 1992]) Show that the system(20.12)characterizes the solution to the following generalizations of the LS problemand the problem of finding the minimum norm solution to an underdeterminedsystem:(20.13)(20.14)20.2.
(RESEARCH PROBLEM) Find a formula for the backward error of anarbitrary approximation to the minimum 2-norm solution of an underdetermined system. That is, for Awith rank(A) = m < n, findm(y) := min{ e : y is the minimum 2-norm solution to (A + ∆A)y = b + ∆b,where ||∆A||2 < e||A||2, ||∆b||2 < e||b||2 }.Previous Home NextChapter 21Vandermonde SystemsWe began, 25 years ago,theThe original motivation came fromcomputation of Gaussmoments ofto take up [the conditioning of]class of Vandermonde matrices.unpleasant experiences with thetype quadrature rules from thethe underlying weight function.— WALTER GAUTSCHI, How (Un)stable are Vandermonde Systems? (1990)Extreme ill-conditioning of the [Vandermonde] linear systemswill eventually manifest itself as n increases by yieldingan error curve which is not sufficiently Ievelled on the current reference .
. .or more seriously fails to have the correct number of sign changes.— M. ALMACANY, C. B. DUNHAM, and J. WILLIAMS,Discrete Chebyshev Approximation by Interpolating Rationals (1984)425V ANDERMONDE S YSTEMS426A Vandermonde matrix is defined in terms of scalars α0, α1, . . . . αnbyVandermonde matrices play an important role in various problems, such asin polynomial interpolation. Suppose we wish to find the polynomial pn(x) =anxn + an-1xn–1 + · · · + a0 that interpolates to the datafor distinctpoints αi , that is, pn(α i ) = fi , i = 0:n. Then the desired coefficient vectora = [a0, a1, .
. . ,an ] T is the solution of the dual Vandermonde systemV T a = f (dual).The primal systemVx = b (primal)represents a moment problem, which arises, for example, when determiningthe weights for a quadrature rule: given moments bi find weights xi such thatBecause a Vandermonde matrix depends on only n+ 1 parameters and hasa great deal of structure, it is possible to perform standard computations withreduced complexity. The easiest algorithm to derive is for matrix inversion.21.1.
Matrix InversionAssume that V is nonsingular and let Vthe equation WV = I may be written-1= W =The ith row ofThese equations specify a fundamental interpolation problem that is solvedby the Lagrange basis polynomial:(21.1)The inversion problem is now reduced to finding the coefficients of li (x). It isclear from (21.1) that V is nonsingular iff the αi are distinct. It also followsfrom (21. 1) that V -1 is given explicitly by42721.1 M ATRIX I NVERSIONwhere σk(y1, . .
. , yn) denotes the sum of all distinct products of k of the arguments y1 , . . . , yn (that is σ k is the kth elementary symmetric function). Anefficient way to find the wij is first to form the master polynomialand then to recover each Lagrange polynomial by synthetic division:The scalars qi (α i ) can be computed by Homer’s rule as the coefficients of qiare formed.Algorithm 21.1. Given distinct scalars α0, α1, . . .
. ,αncomputes W =this algorithm% Stage 1: Construct the master polynomial.a0 = −α0; a1 = 1for k = l:nak+l = 1for j = k: –1: 1a j = a j – 1 - αk a jend% Stage 2: Synthetic division.for i = 0:nwin = l; s = lfor j = n – 1: –1:0w i j = a j + l + αi w i , j + ls = αi s + wijendw(i, : ) = w(i, : )/sendCost: 6n2 flops.The O(n2) complexity is optimal, since the algorithm has n2 output values,each of which must partake in at least one operation.Vandermonde matrices have the deserved reputation of being ill conditioned.
The ill conditioning is a consequence of the monomials being a poorV ANDERMONDE S YSTEMS428Table 21.1. Bounds and estimates forαi(VI):(V2):(V3):(V4):(V5):1/i+larbitraryαi > 0equispaced [0,1]equispaced [– 1, 1]Reference[428, 1990][1035, 1994][429, 1988][428, 1990][425, 1975](V6):(V7):Chebyshev nodes [–1, 1]roots of unity[425, 1975]well known.Bound or estimatebasis for the polynomials on the real line. A variety of bounds for the condition number of a Vandermonde matrix have been derived by Gautschi andhis co-authors. Let Vn = V(α0, α1, .
. . , α n–1 )For arbitrary distinctpoints αi ,(21.3)with equality on the right whenfor all j with a fixed θ (inparticular, when αj > 0 for all j) [424, 1962], [426, 1978]. Note that theupper and lower bounds differ by at most a factor 2 n–1. More specific boundsare given in Table 21.1, on which we now comment.Bound (V1) and estimate (V4) follow from (21.3).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
