Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 76
Текст из файла (страница 76)
The lemma is their Lemma 8.2.11, and can be shown to beequivalent to a result of Stewart [943, 1977, Thin. 5.3].Other methods for solving the LS problem not considered in this chapterinclude those of Peters and Wilkinson [826, 1970], Cline [215, 1973], andPlemmons [835, 1974], all of which begin by computing an LU factorizationof the rectangular matrix A. Error bounds for these methods can be derivedusing results from this chapter and Chapters 9 and 18.In this chapter we have not specifically treated LS problems whose coefficient matrices have rows varying greatly in norm, and we have not consideredweighted LS problems minx ||D(b — Ax)||2, where D = diag(d i ).
Error analysis for Householder QR factorization with column pivoting applied to badlyrow-scaled problems is given by Powell and Reid [840, 1969]. Methods anderror analysis for weighted LS problems are given by Barlow [59, 1988], Bar-412T HE LEAST S QUARES P ROBLEMlow and Handy [64, 1988], Barlow and Vemulapati [66, 1992], Gulliksson andWedin [489, 1992], Gulliksson [487, 1994], [488, 1995], and Van Loan [1042,1985].Another topic not considered here is the constrained LS problem, wherex is required to satisfy linear equality and/or inequality constraints. Numerical methods are described by Lötstedt [713, 1984] and Golub and VanLoan [470, 1989, 12. 1], and perturbation theory is developed by Eldén [350,1980], Lötstedt [712, 1983], and Wedin [1070, 1985].19.9.1. LAPACKDriver routine xGELS solves the full rank LS problem by Householder QRfactorization.
It caters for multiple right-hand sides, each of which defines aseparate LS problem. Thus, xGELS solves min{ ||B – AX||F : X},where AThis routine does not return anyerror bounds, and iterative refinement is not supported for LS problems inLAPACK.Driver routines xGELSX and xGELSS solve the rank-deficient LS problemwith multiple right-hand sides, using, respectively, a complete orthogonal factorization (computed via QR factorization with column pivoting) and theSVD.LAPACK also contains routines for solving the linearly constrained LSproblem (xGGLSE) and a generalized form of weighted LS problem (xGGGLM).Problems19.1.
Show that any solution to the LS problem minx ||b – Ax||2 satisfies thenormal equations A T Ax = A T b. What is the geometrical interpretation ofthese equations?19.2. Prove Theorem 19.3.19.3. The pseudo-inverse Xcan be defined as theunique matrix satisfying the four Moore–Penrose conditions(i) AXA = A,(iii) AX = (AX)T,(ii) XAX = X,(iv) XA = (XA) T .Let A = UΣ VT be an SVD, with Σ = diag(σi ) and let r = rank(A ).
Showthat X = V diag(σ1-l, . . . ,σr-l. 0, . . . , 0)UT satisfies (i)-(iv) and hence is thepseudo-inverse of A. Show that (A+)+ = A.19.4. Show that the pseudo-inverse A+ of Asolves the problemP ROBLEMS413Is the solution unique?19.5. Prove a result analogous to Theorem 19.3 for the MGS method, asdescribed in 19.3.19.6. Consider the LS problem minx||b – Ax||2, whereLet bethe computed LS solution obtained from the normal equations method andx the exact solution, and defineUsing (19.12) and(19.13) show that a bound holds of the form19.7. Prove (19.18) and (19.19).19.8. (Waldén, Karlson, and Sun [1060, 1995]) Partially complete the gap inTheorem 19.5 by evaluating η F (0) for the casethat is, ∆b = 0.19.9. Prove (19.21).Previous Home NextChapter 20Underdetermined SystemsI’m thinking of two numbers.Their average is 3.What are the numbers?— CLEVE B.
MOLER, The World’s Simplest Impossible Problem (1990)This problem arises in important algorithmsused in mathematical programming . . .In these cases, B is usually very large and sparse and,because of storage difficulties,it is often uneconomical to store and access Q1 . . .Sometimes it has been thought that [the seminormal equations method]could be disastrously worse than [the Q method] . .
.It is the purpose of this note to show that such algorithms arenumerically quite satisfactory.— C. C. PAIGE, An Error Analysis of aMethod for Solving Matrix Equations (1973)415416U NDERDETERMINED S YSTEMSHaving considered well-determined and overdetermined linear systems, wenow turn to the remaining class of linear systems: those that are underdetermined.20.1. Solution MethodsConsider the underdetermined system Ax = b, where AThe system can be analysed using a QR factorizationwith m < n.(20.1)where Qis orthogonal and Ris upper triangular.
(Wecould, alternatively, use an LQ factorization of A, but we will keep to thestandard notation.) We haveb = Ax = [RT 0]Q T x = R T y 1 ,(20.2)whereIf A has full rank then yl = R-Tb is uniquely determined and all solutions ofAx = b are given byThe unique solution xLS that minimizes ||x||2 is obtained by setting y2 = 0.We have(20.3)(20.4)where A + = A T (AA T ) – 1 is the pseudo-inverse of A. Hence x L S can becharacterized as xLS = ATy, where y solves the normal equations AATy = b.Equation (20.3) defines one way to compute xLS . We will refer to thismethod as the “Q method”. When A is large and sparse it is desirable toavoid storing and accessing Q, which can be expensive.
An alternative method20.2 P ERTURBATION T HEORY417with this property uses the QR factorization (20.1) but computes x LS asx LS = A T y, where(20.5)RTRy = b(cf. (20.4)). These latter equations are called the seminormal equations (SNE).As the “semi” denotes, however, this method does not explicitly form AAT,which would be undesirable from the standpoint of numerical stability. Notethat equations (20.5) are different from the equations RTRx = ATb for anoverdetermined least squares (LS) problem, where A = Q[ RT 0] Twith m > n, which are also called seminormal equations (see §19.6).20.2. Perturbation TheoryA componentwise perturbation result for the minimum 2-norm solution to anunderdetermined system is readily obtained.Theorem 20.1 (Demmel and Higham). Let Arank andSuppose ||A+ ∆A||2 < 1 and(m < n) be of fullIf x and y are the minimum 2-norm solutions to Ax = b and (A + ∆A)y =b + ∆b, respectively, then, for any monotonic norm,(20.6)For any Hölder p-norm, the bound is attainable to within a constant factordepending on n.Proof.
The perturbed matrix A + ∆A = A(I + A+ ∆A) has full rank, sowe can manipulate the equationy = (A+ ∆A)T((A + ∆A)(A + ∆A)T)-l(b + ∆b)to obtainy – x = (I – A+ A) ∆AT(AAT)-lb + A+(∆b – ∆Ax) + O( 2 )= (I - A + A) ∆ATA+Tx + A+(∆b – ∆Ax) + O(ε2).(20.7)The required bound follows on using absolute value inequalities and takingnorms. That the bound is attained to within a constant factor depending onn for Holder p-norms is a consequence of the fact that the two vectors on theright-hand side of (20.7) are orthogonal.418U NDERDETERMINED S YSTEMSTwo special cases are worth noting, for later use. We will use the equality||I - A+A||2 = min{l, n - m}, which can be derived by consideration of the QRfactorization (20.1), for example.
If E = |A|H, where H is a given nonnegativematrix, and f = |b|, then we can put (20.6) in the form(20.8)wherecond2 ( A) = || |A+||A| ||2.Note that cond2 (A ) is independent of the row scaling of A (cond2(D A) =cond 2 (A) for nonsingular diagonal D). If E = ||A||2emeTn and f = ||b||2em,where em denotes the m-dimensional vector of 1s, then(20.9)The following analogue of Lemma 19.6 will be needed for the error analysis in the next section, It says that if we perturb the two occurrences ofA in the normal equations AATx = b differently, then the solution of theperturbed system is the solution of normal equations in which there is onlyone perturbation of A and, moreover, this single perturbation is no larger, inthe normwise or componentwise sense, than the two perturbations we startedwith.Lemma 20.2full rank and suppose(A + ∆A 1 )and Schwetlick).
Let A= b,(m < n) be of= (A + ∆A 2 ) TAssume that 3 max (||Α+∆A1||2, ||A+ ∆A2||2) < 1. Then there is a vectorand a perturbation ∆A with∆A = ∆AIGI + ∆A 2 G 2 ,Gl + G2 = I,such thatthat is,is the minimum 2-norm solution to (A + ∆A)x = b.Proof. The proof is similar to that of Lemma 19.6, but differs in somedetails. If = (A+ ∆A 2 ) T = 0 we take ∆A = ∆A2. otherwise, we set∆ A := ∆AIP + ∆A2(I – P) =: ∆A2 + HP,20.3 E RROR ANALYSISwhere419and H = ∆A1 - ∆A2. We havewhere β =which shows that we need to setcheck that (A + ∆A ) = b, we evaluateas required. The vectorTois undefined if β = 0. Butwhich is positive if 3 max (||A+ ∆A 1 || 2 , ||A + ∆A 2 || 2 ) < 1.Note that in Lemma 20.2 we have the normwise bound20.3.
Error AnalysisWe now consider the stability of the Q method and the SNE method. For bothmethods we assume that the QR factorization is computed using Householderor Givens transformations.Before presenting the results we define a measure of stability. The componentwise backward error for a minimum-norm underdetermined systemAx = b is defined ass.t. y is the min. norm solution to (A+ ∆A)y = b + ∆b },Note the requirement in this definition that y be the minimum norm solution;the usual componentwise backward error(see (7.6)) is a generallysmaller quantity. Let us say that a method is row-wise backward stable if itproduces a computed solution for which the componentwise backward error420U NDERDETERMINED S YSTEMSis of order u, where E = |A|emeTn and f = |b|.
This condition requiresthat solve a perturbed minimum norm problem in which the perturbationsto the ith row of A are small compared with the norm of that row (similarlyfor b); cf. the discussion of componentwise backward errors in 7.2.Theorem 20.3. Let Awith rank(A) = m < n, and assume thata rendition of the form cond2 (A)mnγ c n < 1 holds. Suppose the underdetermined system AX = b is solved in the minimum 2-norm sense using the Qmethod. Then the computed solution is the minimum 2-norm solution to(A+ ∆A)x = b, whereand||G||F = 1.Proof. The Q method solves the triangular system R T y l = b and thenforms x =Assuming the use of Householder QR factorization,from Theorem 18.4 we have thatfor some orthogonal matrix Q, where ||∆A0||F < mγ cn||A||Fwith ||G0 ||F = 1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
