Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 69
Текст из файла (страница 69)
The vectorwhereso thatdiffers from x only in elements i and j. We haveand similarly forHenceWe takeFor the next result we need the notion of disjoint Givens rotations. Rotations Gi1 ,j 1, . . . , Gir,jr are disjoint if is js and js jt for st . Disjoint rotations are “nonconflicting” and therefore commute; it matters neither mathematically nor numerically in which order the rotations are applied.
(Disjoint374QR F ACTORIZATIONrotations can therefore be applied in parallel, though that is not our interest here. ) Our approach is to take a given sequence of rotations and reorderthem into groups of disjoint rotations. The reordered algorithm is numericallyequivalent to the original one, but allows a simpler error analysis.As an example of a rotation sequence already ordered into disjoint groups,consider the following sequence and ordering illustrated for a 6 × 5 matrix:Here, an integer k in position (i, j) denotes that the (i, j) element is eliminatedon the kth step by a rotation in the (j, i) plane, and all rotations on the kthstep are disjoint. For an m × n matrix with m > n there are r = m + n — 2stages, and the Givens QR factorization can be written as Wr Wr-1 .
. . W1 A =R, where each Wi is a product of at most n disjoint rotations. It is easy to seethat an analogous grouping into disjoint rotations can be done for the schemeillustrated at the start of this section.Lemma 18.8. Consider the sequence of transformationsA k+1 = WkAk,k = 1:r,where A1 =and each Wk is a product of disjoint Givens rotations.Assume that the individual Givens rotations are performed using computedsine and cosine values related to the exact values defining the Wk by (18.15).Then the computed matrix Âr+1 satisfiesÂr+1 = QT(A + ∆A),where QT = Wr Wr–1 . .
. W1 and ∆A satisfies the normwise and componentwise bounds(In fact, we can take G = m - 1 eeT, where e = [1, 1, . . . , 1] T.) In the specialcase n = 1, so that A = a, we have â(r + 1) = (Q + ∆Q)T a with ||∆ Q||F < γ c r .Proof. The proof is analogous to that of Lemma 18.3, so we offer onlya sketch. First, we consider the jth column of A, aj, which undergoes the18.6 I TERATIVE R EFINEMENT375transformations= W r . . .
W1 aj. By Lemma 18.7 and the disjointnessof the rotations, we havewhere each ∆Wk depends on j and satisfies ||∆ Wk|| 2 <we obtainUsing Lemma 3.6(18.16)HenceThe inequalities (18. 16) for j = 1:n imply thatwhere ||G||F = 1. The result for n = 1 is proved as in Lemma 18.3. IIWe are now suitably equipped to give a result for Givens QR factorization.Theorem 18.9. Letbe the computed upper trapezoidal QR factor(m > n) obtained via the Givens QR algorithm, with anyofstandard choice and ordering of rotations.
Then there exists an orthogonalsuch thatwith ||∆A||F < γc(m+n)||A||F and |∆ A| < mγ c(m+n)G|A|, ||G||F = 1. (Thematrix Q is a product of Givens rotations, the kth of which corresponds to theexact application of the kth step of the algorithm to Âk.)It is interesting that the error bounds for QR factorization with Givensrotations are a factor n smaller than those for Householder QR factorization.This appears to be an artefact of the analysis, and we are not aware of anydifference in accuracy in practice.18.6. Iterative RefinementConsider a nonsingular linear system Ax = b, whereSuppose wesolve the system using a QR factorization A = QR computed using Householder or Givens transformations (thus, x is obtained by solving Rx = QTb).Theorem 18.5, and its obvious analogue for Givens rotations, show thatsatisfies(18.17)376QR F ACTORIZATIONwhere ||G||F = 1 and p is a low-degree polynomial.
Hencehas a normwisebackward error of order p(n)u. However, since G is a full matrix, (18.17) suggests that the componentwise relative backward errorneed not besmall. In fact, we know of no nontrivial class of matrices for which Householderor Givens QR factorization is guaranteed to yield a small componentwise relative backward error.Suppose that we carry out a step of fixed precision iterative refinement, toobtainThe form of the bound (18.17) enables us to invoke Theorem 11.4.We conclude that the componentwise relative backward errorafter one step of iterative refinement will be small as long as A is not too illconditioned andis not too badly scaled.
This conclusion is similar tothat for Gaussian elimination with partial pivoting (GEPP), except that forGEPP there is the added requirement that the LU factorization not sufferlarge element growth.Recall from (18.7) that we also have a backward error result in which onlyA is perturbed. The analysis of $11.1 is therefore applicable, and analoguesof Theorems 11.1 and 11.2 hold in which η =The performance of QR factorization with fixed precision iterative refinement is illustrated in Tables 11. 1–1 1.3 in §11.2. The performance is as predicted by the analysis.
Notice that the initial componentwise relative backward error is large in Table 11.2 but that iterative refinement successfullyreduces it to the roundoff level (despitebeing huge). It isworth stressing that the QR factorization yielded a small normwise relativebackward error in each examplein fact), as we know it must.18.7. Gram–Schmidt OrthogonalizationThe oldest method for computing a QR factorization is the Gram-Schmidtorthogonalization method. It can be derived directly from the equation A =(Gram-Schmidt does not computeQR, where A,andthe m x m matrix Q in the full QR factorization and hence does not providea basis for the orthogonal complement of range(A)). Denoting by a j and q jthe jth columns of A and Q, respectively, we haveyields, since Q has orthonormal columns,Premultiplying byi = 1:j — 1.
Further,18.7 G RAM -S CHMIDT O RTHOGONALIZATION377whereHence we can compute Q and R a column at a time. To ensure that rjj > 0we require that A has full rank.Algorithm 18.10 (classical Gram–Schmidt). Givenof rank nthis algorithm computes the QR factorization A = QR, where Q is m × n andR is n × n, by the Gram-Schmidt method.for j = 1:nfor i = 1:j –1endCost: 2mn2 flops (2n 3/3 flops more than Householder QR factorizationwith Q left in factored form).In the classical Gram–Schmidt method (CGS), aj appears in the computation only at the jth stage. The method can be rearranged so that as soon asqj is computed, all the remaining vectors are orthogonalized against qj. Thisgives the modified Gram-Schmidt method (MGS).Algorithm 18.11 (modified Gram-Schmidt). Givenof rank nthis algorithm computes the QR factorization A = QR, where Q is m × n andR is n × n, by the MGS method.endendCost: 2mn2flops.378QR F ACTORIZATIONIt is worth noting that there are two differences between the CGS andMGS methods.
The first is the order in which the calculations are performed:in the modified method each remaining vector is updated once on each stepinstead of having all its updates done together on one step. This is purely amatter of the order in which the operations are performed. Second, and morecrucially in finite precision computation, two different (but mathematicallyequivalent ) formulae for rkj are used: in the classical method, rkj =which involves the original vector a j, whereas in the modified method a j isreplaced in this formula by the partially orthogonalized vectorAnotherway to view the difference between the two Gram–Schmidt methods is viarepresentations of an orthogonal projection; see Problem 18.7.The MGS procedure can be expressed in matrix terms by defining Ak =MGS transforms A1 = A into An+1 = Q by thesequence of transformations Ak = Ak+ 1Rk, where Rk is equal to the identityexcept in the kth row, where it agrees with the final R.
For example, if n = 4and k = 3,Thus R = Rn . . . R 1 .The Gram-Schmidt methods produce Q explicitly, unlike the Householderand Givens methods, which hold Q in factored form. While this is a benefit,in that no extra work is required to form Q, it is also a weakness, becausethere is nothing in the methods to force the computedto be orthonormal inthe face of roundoff.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
