Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 73
Текст из файла (страница 73)
The conclusionis that, unless A is very ill conditioned, the residual b –will not exceedthe larger of the true residual r = b – AX and a constant multiple of the errorin evaluating fl(r)—a very satisfactory result.19.3. Solution by the Modified Gram–Schmidt MethodThe modified Gram-Schmidt (MGS) method can be used to solve the LSproblem. However, we must not compute x from x = R–1(QTb), because thelack of orthonormality of the computedwould adversely affect the stability.Instead we apply MGS to the augmented matrix [A b]:19.4 T HE N ORMAL E QUATIONS397We haveSince qn+l is orthogonal to the columns of Q1, ||b – Ax|| 22 = ||Rx – z|| 22 + ρ2 ,so the LS solution is x = R-lZ.
Of course, z = QT1b, but z is now computedas part of the MGS procedure instead of as a product between QT and bBjörck [107, 1967] shows that this algorithm is forward stable, in theis as small as that for a normsense that the forward error ||x –wise backward stable method. It has recently been shown by Björck andPaige [119, 1992] that the algorithm is, in fact, normwise backward stable(see also Björck [115, 1994]), that is, a normwise result of the form in Theorem 19.3 holds.
Moreover, a componentwise result of the form in Theorem 19.3holds too—see Problem 19.5. Hence the possible lack of orthonormality ofdoes not impair the stability of the MGS method as a means for solving theLS problem.19.4. The Normal EquationsThe oldest method of solving the LS problem is to form and solve the normalequations, A T Ax = A T b. Assuming that A has full rank, we can use thefollowing procedure:Form C = ATA and c = ATb.Compute the Cholesky factorization C = RTR.Solve RTy = c, Rx = y.Cost: : n 2(m + n/3) flops.If m >> n, the normal equations method requires about half as manyflops as the Householder QR factorization approach (or the MGS method).However, it has less satisfactory numerical stability properties. There aretwo problems.
The first is that information may be lost when= fl(ATA) isformed-essentially because forming the cross product is a squaring operationthat increases the dynamic range of the data. A simple example is the matrixfor whichEven though A is distance approximately ε from a rank-deficient matrix, andhence unambiguously full rank if εthe computed cross product is398T HE LEAST S QUARES P ROBLEMsingular. In general, whenever K2(A) > u –1/2 we can expectto be singularor indefinite, in which case Cholesky factorization is likely to break down(Theorem 10.7).The second weakness of the normal equations method is more subtle andis explained by a rounding error analysis. In place of C = ATA and c = ATbwe computeBy Theorems 10.3 and 10.4, the computed Cholesky factorsatisfyand solution(19.11)Overall,(19.12)By boundingwith the aid of (19.11), we find that(19.13a)(19.13b)These bounds show that we have solved the normal equations in a backwardstable way, as long as ||A||2||b||2 ||ATb||2 .
But if we try to translate thisresult into a backward error result for the LS problem itself, we find thatthe best backward error bound contains a factor K2(A) [569, 1987]. The bestforward error bound we can expect, in view of (19.13), is of the form(19.14)(since K2 (A T A) = K2 (A)2). Now we know from Theorem 19.1 that the sensitivity of the LS problem is measured by K2(A)2 if the residual is large, butby K2(A) if the residual is small. It follows that the normal equations methodhas a forward error bound that can be much larger than that possessed by abackward stable method.A mitigating factor for the normal equations method is that, in view ofTheorem 10.6, we can replace (19.14) by the (not entirely rigorous) bound19.5 I TERATIVE R EFINEMENT399where A = BD, with D = diag(||A(:, i)||2), so that B has columns of unit2-norm.
Van der Sluis’s result (Theorem 7.5) shows thatHence the normal equations method is, to some extent, insensitive to poorcolumn scaling of A.Although numerical analysts almost invariably solve the full rank LS problem by QR factorization, statisticians frequently use the normal equations(though perhaps less frequently than they used to, thanks to the influence ofnumerical analysts). The normal equations do have a useful role to play. Inmany statistical problems the regression matrix is contaminated by errors ofmeasurement that are very large relative to the roundoff level; the effects ofrounding errors are then likely to be insignificant compared with the effectsof the measurement errors, especially if IEEE double precision (as opposed tosingle precision) arithmetic is used.The normal equations (NE) versus (Householder) QR factorization debatecan be summed up as follows.●The two methods have a similar computational cost if m n, but theNE method is up to twice as fast for m >> n.
(This statement assumesthat A and b are dense; for details of the storage requirements andcomputational cost of each method for sparse matrices, see, for example,Björck [116, 1996] and Heath [510, 1984].)●The QR method is always backward stable. The NE method is guaranteed to be backward stable only if A is well conditioned.●The forward error from the NE method can be expected to exceed thatfor the QR method when A is ill conditioned and the residual of the LSproblem is small.●The QR method lends itself to iterative refinement, as described in thenext section.
Iterative refinement can be applied to the NE method, butthe rate of convergence inevitably depends on K2(A)2 instead of K2(A).19.5. Iterative RefinementAs for square linear systems, iterative refinement can be used to improvethe accuracy and stability of an approximate LS solution. However, for theLS problem there is more than one way to construct an iterative refinementscheme.By direct analogy with the square system case, we might consider theschemeT HE LEAST S QUARES P ROBLEM4001. Compute r = b –2.
Solve the LS problem mind ||Ad–r||2 .3. Update y =+ d.(Repeat from step 1 if necessary, withreplaced by y).This scheme is certainly efficient-a computed QR factorization (for example)of A can be reused at each step 2. Golub and Wilkinson [466, 1966] investigated this scheme and found that it works well only for nearly consistentsystems.An alternative approach suggested by Björck [106, 1967] is to apply iterative refinement to the augmented system (19.3), so that both x and r arerefined simultaneously.
Since this is a square, nonsingular system, existing results on the convergence and stability of iterative refinement can be applied,and we would expect the scheme to work well. To make precise statementswe need to examine the augmented system method in detail.For the refinement steps we need to consider an augmented system withan arbitrary right-hand side:r + Ax = f,(19.15a)ATr = g.(19.15b)If A has the QR factorizationwhere Rthen (19. 15) transforms toThis system can be solved as follows:h = R - T g,x = R –l (d l – h).The next result describes the effect of rounding errors on the solution process.19.5 I TERATIVE R EFINEMENT401Theorem 19.4. Let Abe of full rank n < m and suppose the augmented system (19.15) is solved using a Householder QR factorization of Aas described above.
The computed and satisfywherewith ||G||F = 1, ||Hi ||F = 1, i = 1:3.Proof. The proof does not involve any new ideas and is rather tedious,so we omit it.Consider first fixed precision iterative refinement. Theorem 19.4 impliesto (19.15) satisfiesthat the computed solutionUnfortunately, this bound is not of a form that allows us to invoke Theorem 11.4. However, we can apply Theorem 11.3, which tells us that thecorrected solutionobtained after one step of iterative refinementsatisfies(19.16)Here,to the originalbound (19.16).be included inand sobe writtendenotes the residual of the augmented system correspondingcomputed solution.
We will make two simplifications to the= O(u), the first term in the bound mayFirst, since2the O(u ) term. Second, (19.16) yieldsO(u)With these two simplifications, (19.16) may402T HE LEAST S QUARES P ROBLEMIn view of the Oettli–Prager result (Theorem 7.3) this inequality tells us that,asymptotically, the solution produced after one step of fixed precision iterative refinement has a small componentwise relative backward error withrespect to the augmented system. However, this backward error allows thetwo occurrences of A in the augmented system coefficient matrix to be perturbed differently, and thus is not a true componentwise backward error forthe LS problem. Nevertheless, the result tells us that iterative refinement canbe expected to produce some improvement in stability.
Note that the bound inTheorem 19.2 continues to hold if we perturb the two occurrences of A in theaugmented system differently. Therefore the bound is applicable to iterativerefinement (with E = |A|, f = |b|), and so we can expect iterative refinement to mitigate the effects of poor row and column scaling of A. Numericalexperiments show that these predictions are borne out in practice [549, 1991].Turning to mixed precision iterative refinement, we would like to applythe analysis of $11.1, with “Ax = b“ again identified with the augmentedsystem. However, the analysis of 11.1 requires a backward error result inwhich only the coefficient matrix is perturbed (see (11.1)).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
