Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 71
Текст из файла (страница 71)
The idea of exploiting disjoint rotations in the error analysis was developed by Gentleman [436, 1975], who gave a normwiseanalysis that is simpler and produces smaller bounds than Wilkinson’s (ournormwise bound in Theorem 18.9 is essentially the same as Gentleman’s).For more details of algorithmic and other aspects of Householder andGivens QR factorization, see Golub and Van Loan [470, 1989, §5.2].The error analysis in §18.3 is a refined and improved version of analysisthat appeared in the technical report [546, 1990] and was quoted withoutproof in Higham [549, 1991].The WY representation for a product of Householder transformationsshould not be confused with a genuine block Householder transformation.Schreiber and Parlett [904, 1988] define, for a given(m > n), the“block reflector that reverses the range of Z” asIf n = 1 this is just a standard Householder transformation.
A basic task is(m > n) find a block reflector H such thatas follows: givenSchreiber and Parlett develop theory and algorithms for block reflectors, inboth of which the polar decomposition plays a key role.Sun and Bischof [977, 1995] show that any orthogonal matrix can be expressed in the form Q = I – YSYT, even with S triangular, and they explorethe properties of this representation.Another important use of Householder matrices, besides computation ofthe QR factorization, is to reduce a matrix to a simpler form prior to iterative computation of eigenvalues (Hessenberg or tridiagonal form) or singularvalues (bidiagonal form). For these two-sided transformations an analogue ofLemma 18.3 holds with normwise bounds (only) on the perturbation.
Erroranalyses of two-sided application of Householder transformations is given byOrtega [811, 1963] and Wilkinson [1086, 1962], [1089, 1965, Chap. 6].Mixed precision iterative refinement for solution of linear systems by Householder QR factorization is discussed by Wilkinson [1090, 1965, §10], who notesthat convergence is obtained as long as a condition of the form cn κ 2 (A)u < 1holds.18.9 N OTESANDR EFERENCES385Fast Givens rotations can be applied to a matrix with half the numberof multiplications of conventional Givens rotations, and they do not involvesquare roots. They were developed by Gentleman [435, 1973] and Hammarling [498, 1974].
Fast Givens rotations are as stable as conventional ones-seethe error analysis by Parlett in [820, 19 8 0 , §6.8.3], for example-but, forthe original formulations, careful monitoring is required to avoid overflow.Rath [861, 1982] investigates the use of fast Givens rotations for performingsimilarity transformations in solving the eigenproblem. Barlow and Ipsen [65,1987] propose a class of scaled Givens rotations suitable for implementationon systolic arrays, and they give a detailed error analysis.
Anda and Park [16,1994] develop fast rotation algorithms that use dynamic scaling to avoid overflow.Rice [870, 1966] was the first to point out that the MGS method producesa more nearly orthonormal matrix than the CGS method in the presenceof rounding errors. Björck [107, 1967] gives a detailed error analysis, proving(18.21) and (18.22) but not (18.23), which is an extension of the correspondingnormwise result of Björck and Paige [119, 1992]. Björck and Paige give adetailed assessment of MGS versus Householder QR factorization.The difference between the CGS and MGS methods is indeed subtle.Wilkinson [1095, 1971] admitted that “I used the modified process for manyyears without even noticing explicitly that I was not performing the classicalalgorithm.”The orthonormality of the matrixfrom Gram-Schmidt can be improvedby reorthogonalization, in which the orthogonalization step of the classicalor modified method is iterated.
Analyses of Gram–Schmidt with reorthogonalization are given by Abdelmalek [2, 1971], Ruhe [883, 1983], and Hoffmann [578, 1989]. Daniel, Gragg, Kaufman, and Stewart [263, 1976] analysethe use of classical Gram–Schmidt with reorthogonalization for updating aQR factorization after a rank one change to the matrix.The mathematical and numerical equivalence of the MGS method withHouseholder QR factorization of the matrixwas known in the 1960s(see the Björck quotation at the start of the chapter) and the mathematicalequivalence was pointed out by Lawson and Hanson [695, 1974, Ex. 19.39].A block Gram-Schmidt method is developed by Jalby and Philippe [608,1991 ] and error analysis given.
See also Björck [115, 1994 ], who gives anup-to-date survey of numerical aspects of the Gram–Schmidt method.For more on Gram-Schmidt methods, including historical comments, seeBjörck [116, 1996].One use of the QR factorization is to orthogonalize a matrix that, becauseof rounding or truncation errors, has lost its orthogonality; thus we computeA = QR and replace A by Q. An alternative approach is to replace(m > n) by the nearest orthonormal matrix, that is, the matrix Q that solves{ ||A - Q|| : QTQ = I} = min.
For the 2- and Frobenius norms, the optimal386QR F ACTORIZATIONQ is the orthonormal polar factor U of A, where A = UH is a polar decomposition:has orthonormal columns andis symmetricpositive semidefinite. If m = n, U is the nearest orthogonal matrix to A in anyunitarily invariant norm, as shown by Fan and Hoffman [361, 1955]. Chandrasekaran and Ipsen [196, 1994] show that the QR and polar factors satisfyunder the assumptions that A has full rank andcolumns of unit 2-norm and that R has positive diagonal elements. Sun [974,1995] proves a similar result and also obtains a bound for ||Q — U||F in termsof ||ATA – I||F. Algorithms for maintaining orthogonality in long products oforthogonal matrices, which arise, for example, in subspace tracking problemsin signal processing, are analysed by Edelman and Stewart [347, 1993] andMathias [733, 1995].Various iterative methods are available for computing the orthonormalpolar factor U, and they can be competitive in cost with computation of aQR factorization.
For more details on the theory and numerical methods, seeHigham [530, 1986], [539, 1989], Higham and Papadimitriou [567, 1994], andthe references therein.A notable omission from this chapter is a treatment of rank-revealing QRfactorizations-ones in which the rank of A can readily be determined fromR. This topic is not one where rounding errors play a major role, and hence itis outside the scope of this book. Pointers to the literature include Golub andVan Loan [470, 1989, §5.4], Chan and Hansen [192, 1992], and Björck [116,199 6]. A column pivoting strategy for the QR factorization, described inProblem 18.5, ensures that if A has rank r then only the first r rows of Rare nonzero.
A perturbation theorem for the QR factorization with columnpivoting is given by Higham [540, 1990]; it is closely related to the perturbation theory in §10.3.1 for the Cholesky factorization of a positive semidefinitematrix.18.9.1. LAPACKLAPACK contains a rich selection of routines for computing and manipulating the QR factorization and its variants. Routine xGEQRF computes the QRfactorization A = QR of an m × n matrix A by the Householder QR algorithm.
If m < n (which we ruled out in our analysis, merely to simplify thenotation), the factorization takes the form A = Q[ R 1 R2], where R1 is m × mupper triangular. The matrix Q is represented as a product of Householdertransformations and is not formed explicitly. A routine xORGQR (or xUNGQR inthe complex case) is provided to form all or part of Q, and routine xORMQR (orxUNMQR ) will pre- or postmultiply a matrix by Q or its (conjugate) transpose.Routine xGEQPF computes the QR factorization with column pivoting (seeProblem 18.5).An LQ factorization is computed by xGELQF . When A is m × n with m < nP ROBLEMS387it takes the form A = [L 0] Q. It is essentially the same as a QR factorizationof AT and hence can be used to find the minimum 2-norm solution to anunderdetermined system (see §20.1).LAPACK also computes two nonstandard factorization of an m × n A:where L is lower trapezoidal and R upper trapezoidal.Problems18.1.
Find the eigenvalues of a Householder matrix and a Givens matrix.18.2. Letwhere andin Lemma 18.1. Derive a bound forare the computed quantities described18.3. A complex Householder matrix has the formwhereand β = 2/υ∗υ. For givendetermine, if possible, P so that Px = y.show how to18.4. (Wilkinson [1089, 1965, p. 242]) Letand let P be a Householdermatrix such that Px = ±||x||2 e1. Let G1,2, . . . .
Gn– 1 ,n be Givens rotationssuch that Qx := G 1,2. . . Gn- 1 ,nx = ±||x||2 e1. True or false: P = Q?18.5. In the QR factorization with column pivoting, columns are interchangedat the start of the kth stage of the reduction so that, in the notation of (18.1),||xk|| 2 > ||Ck(:, j)||2 for all j. Show that the resulting R factor satisfiesso that, in particular, |r11| > |r22| > . .
. > |rnn|. (These are the sameequations as (10.13), which hold for the Cholesky factorization with completepivoting—why?)18.6. Letbe a product of disjoint Givens rotations. Show that18.7. The CGS method and the MGS method applied to( m > n)compute a QR factorization A = QR,Define the orthogonalprojection Pi =where qi = Q(:, i). Show thatQR F ACTORIZATION388Show that the CGS method corresponds to the operationswhile MGS corresponds to18.8. Let A = [a 1, a2 ]θ, 0 < θ < π/2. (Thus, cosθ :=18.9. (Björck [107,1967])and denote the angle between a 1 and a 2 byShow thatLetwhich is a matrix of the form discussed by Lauchli [692, 196 1].
Assumingthatevaluate the Q matrices produced by the CGS and MGSmethods and assess their orthonormality.18.10. Show that the matrix P in (18.19) has the formwhere Q is the matrix obtained from the MGS method applied to A.18.11. (Björck and Paige [119,1992])For any matrices satisfyingwhere both P11 and P21 have at least as many rows as columns, show thatthere exists an orthonormal Q such that A + ∆A = QR, where(Hint: use the CS decomposition P11 = UCWT, P21 = VSWT, where Uand V have orthonormal columns, W is orthogonal, and C and S are square,nonnegative diagonal matrices with C2 + S 2 = I.