Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 74
Текст из файла (страница 74)
This causes nodifficulties because from Theorem 19.4 we can deduce a normwise result (cf.Problem 7.7):The theory of 11.1 then tells us that mixed precision iterative refinement willconverge as long as the condition number of the augmented system matrix isnot too large and that the rate of convergence depends on this conditionnumber. How does the condition number of the augmented system relate tothat of A? Consider the matrix that results from the scaling(α > 0):(19.17)Björck [106, 1967] shows that the eigenvalues of C(α) are(19.18)( 0’,m - n times,where σi i = 1: n, are the singular values of A, and that(19.19)with minα K2(C(α)) being achieved for α =(see Problem 19.7). HenceC (α) may be much more ill conditioned than A. However, in our analysis19.6 T HE S EMINORMAL E QUATIONS403we are at liberty to take minα K2(C(α )) as the condition number, becausescaling the LS problem according to b – Ax(b – Ax)/α does not changethe computed solution or the rounding errors in any way (at least not ifα is a power of the machine base).
Therefore it is, K2(A) that governs thebehaviour of mixed precision iterative refinement, irrespective of the size ofthe LS residual. As Björck [111, 1990] explains, this means that “in a senseiterative refinement is even more satisfactory for large residual least-squaresproblems.” He goes on to explain that “When residuals to the augmentedsystem are accumulated in precision β -t2, t2 > 2t1, this scheme gives solutionsto full single-precision accuracy even though the initial solution may have nocorrect significant figures.”Iterative refinement can be applied with the MGS method. Björck [106,1967] gives the details and shows that mixed precision refinement works justas well as it does for Householder’s method.19.6.
The Seminormal EquationsWhen we use a QR factorization to solve an LS problem minx ||b – Ax||2,the solution x is determined from the equation Rx = QTb (or via a similarexpression involving Q for the MGS method). But if we need to solve forseveral right-hand sides that are not all available when the QR factorizationis computed, we need to store Q before applying it to the right-hand sides.If A is large and sparse it is undesirable to store Q, as it can be much moredense than A. We can, however, rewrite the normal equations asR T Rx = A T b,which are called the seminormal equations. The solution x can be determinedfrom these equations without the use of Q. Since the cross product matrixATA is not formed explicitly and R is determined stably via a QR factorization, we might expect this approach to be more stable than the normalequations method.Björck [110, 1987] has done a detailed error analysis of the seminormalequations (SNE) method, under the assumption that R is computed by abackward stable method.
His forward error bound for the SNE method is ofthe same form as that for the normal equations method, involving a factor2K2 (A) . Thus the SNE method is not backward stable. Björck considersapplying one step of fixed precision iterative refinement to the SNE method,and he calls the resulting process the corrected seminormal equations (CSNE)method:R T Rx = A T br = b – Ax404T HE LEAST S QUARES P ROBLEMR T Rw = A T ry=x+wIt is important that the normal equations residual be computed as shown, asTAT(b – Ax), and not as ATb – A Ax. Björck derives a forward error boundfor the CSNE method that is roughly of the formHence, if k2(A)2U < 1, the CSNE method has a forward error bound similarto that for a backward stable method, and the bound is actually smaller thanthat for the QR method if k2(A)2U << 1 and r is small. However, the CSNEmethod is not backward stable for all A.19.7. Backward ErrorAlthough it has been known since the 1960s that a particular method forsolving the LS problem, namely the Householder QR factorization method,yields a small normwise backward error (see $19.2), it was for a long timean open problem to obtain a formula for the backward error of an arbitraryapproximate solution.
Little progress had been made towards solving thisproblem until Waldén, Karlson, and Sun [1060, 1995] found an extremelyelegant solution. We will denote by λ min and σmin the smallest eigenvalueof a symmetric matrix and the smallest singular value of a general matrix,respectively.Theorem 19.5 (Waldén, Karlson, and Sun). Let(m > n), band r = b – Ay.
The normwise backward error(19.20)is given bywhereThe backward error (19.20) is not a direct generalization of the usualnormwise backward error for square linear systems, because it minimizes19.7 B ACKWARD E RROR405Table 19.1. LS backward errors and residual for Vandermonde system.instead ofHowever, the Parameter θ allows us some flexibility: taking the limitforces ∆b = 0,giving the case where only A is perturbed.Theorem 19.5 can be interpreted as saying that if λ* > 0 then the backwarderror is essentially that given by Theorem 7.1 for a consistent system. Ifλ* < 0, however, the nearest perturbed system of which y is the LS solutionis inconsistent.
A sufficient condition for λ* < 0 is b range(A) (assumingµ 0), that is, the original system is inconsistent.The formulae given in Theorem 19.5 are unsuitable for computation because they can suffer from catastrophic cancellation when λ * < 0. Instead,the following alternative formula derived in [1060, 1995] should be used (seeProblem 19.9):(19.21a)(19.21b)To illustrate Theorem 19.5, we consider an LS problem with a 25 x 15Vandermonde matrix A =where the pi are equally spaced on [0, 1],and b = Ax with xi = i (giving a zero residual problem).
The conditionnumber k2 (A) = 1.47 x 109 . We solved the LS problem in MATLAB in fourdifferent ways: by using the NE with Cholesky factorization, via HouseholderQR factorization, and via the MGS method, using both the stable approachdescribed in 19.3 and the unstable approach in which QTb is formed asa matrix–vector product (denoted MGS(bad)). The results, including thenorms of the residuals r = b –are shown in Table 19.1. As would beexpected from the analysis in this chapter, the QR and stable MGS methodsproduce backward stable solutions, but the NE method and the unstable MGSapproach do not.As we saw in the analysis of iterative refinement, sometimes we need toconsider the augmented system with different perturbations to A and AT.The next result shows that from the point of view of normwise perturbations406T HE LEAST S QUARES P ROBLEMand componentwise perturbations bounded by |∆ A| < εG|A|, the lack of“symmetry” in the perturbations has little effect on the backward error of y.Lemma 19.6and Schwetlick).
Letconsider the perturbed augmented systemThere is a vector(m > n) andand a perturbation ∆A with∆ A = G1∆ A1 + G2 ∆ A 2 ,(19.22)such thatthat is, y solves the LS problem minx ||(A + ∆A)x – b||2.Proof. If s = b– (A + ∆A1)y = 0 we take ∆A = ∆A1. Otherwise, we set∆ A := P∆ A2 + (I – P)∆ A1 =: ∆A1 + PH,where P = ssT/sTs and H = ∆A2 – ∆A1. We havewhere β = 1 – s T Hy/s T s. ThenNote that (19.22) implies a bound stronger than just ||∆ A||2 < ||∆ A1 ||2 +||∆A 2 || 2 :p = 2, F.Turning to componentwise backward error, the simplest approach is to apply the componentwise backward errorto the augmented system (19.3),settingso as not to perturb the diagonal blocks I and O of the augmented system coefficient matrix.
However, this approach allows A and AT to undergo differentperturbations ∆A1 and ∆A2 withand thus does not give a truebackward error, and Lemma 19.6 is of no help. This problem can be overcome by using a structured componentwise backward error to force symmetry19.8 P ROOFOFW EDIN ’ S T HEOREM407of the perturbations; see Higham and Higham [527, 1992] for details. Oneproblem remains: as far as the backward error of y is concerned, the vector rin the augmented system is a vector of free parameters, so to obtain the truecomponentwise backward error we have to minimize the structure-preservingcomponentwise backward error over all r.
This is a nonlinear optimizationproblem to which no closed-form solution is known. Experiments show thatwhen y is a computed LS solution, r = b – Ay is often a good approximationto the minimizing r [527, 1992], [549, 1991].19.8. Proof of Wedin’s TheoremIn this section we give a proof of Theorem 19.1. We define PA := AA+, theorthogonal projection onto range(A).Lemma 19.7.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
