Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 67
Текст из файла (страница 67)
On the otherhand, if the cnu||A||-pseudospectrum is wholly contained within the unit disc,each A + ∆Ai will have spectral radius less than 1 and the product can beexpected to converge. (Note, however, that if p(A) < 1 and p (B) < 1 it is notnecessarily the case that p(AB) < 1.) To make this heuristic precise, we needan analogue of Theorem 17.1 phrased in terms of the pseudospectrum ratherthan the Jordan form.Theorem 17.2 (Higham and Knight). Suppose thatizable with A =and has a unique eigenvaluemodulus. Suppose thatandis diagonalof largestwhere17.2 B OUNDSFORFINITE P RECISION A RITHMETIC357Figure 17.5.
Pseudospectra for chebspec matrices.X -1 = (yij). If< 1 for= cnu||A||2, where cn is a constant depending only on n, then, provided that a certainterm can be ignored,fl(Am ) = 0.Proof. It can be shown (see [565, 1995]) that the conditions on ||X|| 1 andimply there is a perturbation à = A + ∆A of A with ||∆ A|| 2 =such thatHence, if< 1 thenterm and rearranging givesIgnoring theUsing Theorem 17.1 we have the required result for cn = 4 n 2 (n + 2), sincet = 1.Suppose we compute the eigenvalues of A by a backward stable algorithm,that is, one that yields the exact eigenvalues of A+E, where ||E||2 < cnu||A|| 2 ,with cn a modest constant. (An example of such an algorithm is the QRalgorithm [470, 1989, §7.5]).
Then the computed spectral radius satisfies <In view of Theorem 17.2 we can formulate a rule of thumb-oneM ATRIX P OWERS358that bears a pleasing symmetry with the theoretical condition for convergence:The computed powers of A can be expected to converge to 0 if thespectral radius computed via a backward stable eigensolver is lessthan 1.This rule of thumb has also been discussed by Trefethen and Trummer [1020,19 8 7 ] and Reichel and Trefethen [866, 1992 ]. In our experience the rule ofthumb is fairly reliable whenis not too close to 1. For the matrices used inour examples we have, using MATLAB’S eig function,and we observed convergence of the computed powers for C 8 and C13 + 0.01Iand divergence for the other matrices.17.3.
Application to Stationary IterationAs we saw in the previous chapter, the errors in stationary iteration satisfye k = (M - 1N) k e0, so convergence of the iteration depends on the convergenceof (M – 1 N) k to zero as kWhile the errors in stationary iteration are notprecisely modelled by the errors in matrix powering, because matrix powersare not formed explicitly, the behaviour of the computed powers fl((M - 1N) k )can be expected to give some insight into the behaviour of stationary iteration.For the successive overrelaxation (SOR) example at the start of Chap=0.5(–1) i – j .ter 16, the matrix G = M – 1N is lower triangular with gijThe computed powers of G in MATLAB reach a maximum norm of1028at k = 99 and then decay to zero; the eventual convergence is inevitable inview of the condition (17.12), which clearly is satisfied for this triangular G.An approximation to the u||G||2-pseudospectrum is plotted in Figure 17.6,and we see clearly that part of the pseudospectrum lies outside the unit disk.These facts are consistent with the nonconvergence of the SOR iteration (seeFigure 16.1).That the pseudospectrum of G gives insight into the behaviour of stationary iteration has also been observed by Trefethen [1015, 1990], [1017, 1992],[1018] and Chatelin and Frayssé [203, 1992], but no rigorous results about theconnect ion are available.17.4.
Notes and ReferencesThis chapter is baaed closely on Higham and Knight [565,1995].P ROBLEMS359Figure 17.6. Pseudospectrum for SOR iteration matrix.The analysis for the powers of the matrix (17.2) is modelled on that ofStewart [953, 1994], who uses the matrix to construct a Markov chain whosesecond largest eigenvalue does not correctly predict the decay of the transient.For some results on the asymptotic behaviour of the powers of a nonnegative matrix, see Friedland and Schneider [409, 1980].Another application of matrix powering is in the scaling and squaringmethod for computing the matrix exponential, which uses the identity eA =(eA/m )m together with a Taylor or Padé approximation to eA/m ; see Molerand Van Loan [775, 1978].Problems17.1.
Letbe diagonalizable: A = XΛX -1, Λ = diagstruct a parametrized example to show that the boundcan be arbitrarily weak.Con-17.2. Show that p(|A|)/p(A) can be arbitrarily large for17.3. (RESEARCH PROBLEM) Explain the scalloping patterns in the curves ofnorms of powers of a matrix, as seen, for example, in Figure 17.4. (Considerexact arithmetic, as the phenomenon is not rounding error dependent. )17.4. (RESEARCH PROBLEM) obtain a sharp sufficient condition for fl(Ak)0 in terms of the Schur decomposition ofwith p(A) < 1.Previous HomeChapter 18QR FactorizationAny orthogonal matrix can be written as the product of reflector matrices.Thus the class of reflections is rich enough for all occasionsand yet each member is characterized by a single vectorwhich serves to describe its mirror.— BERESFORD N.
PARLETT, The Symmetric Eigenvalue Problem (1980)A key observation for understanding the numerical properties of themodified Gram–Schmidt algorithm is that it can be interpreted asHouseholder QR factorization applied to the matrix Aaugmented with a square matrix of zero elements on top.These two algorithms are not only mathematically . . .but also numerically equivalent.This key observation, apparently by Charles Sheffield,was relayed to the author in 1968 by Gene Golub.— AKE BJÖRCK, Numerics of Gram-Schmidt Orthogonalization (1994)The great stability of unitary transformations in numerical analysissprings from the fact that both theand the Frobenius norm are unitarily invariant.This means in practice that even when rounding errors are made,no substantial growth takes place in thenorms of the successive transformed matrices.— J. H.
WILKINSON,Error Analysis of Transformations Based on theUse of Matrices of the Form I – 2wwH (1965)361Next362QR F ACTORIZATIONThe QR factorization is a versatile computational tool that finds use in linear equation, least squares and eigenvalue problems. It can be computed inseveral ways, including by the use of Householder transformations and Givensrotations, and by the Gram–Schmidt method. We explore the numerical properties of all three methods in this chapter. We also examine the use of iterativerefinement on a linear system solved with a QR factorization and consider theinherent sensitivity of the QR factorization.18.1. Householder TransformationsA Householder matrix (also known as a Householder transformation, or Householder reflector) is a matrix of the formIt enjoys the properties of symmetry and orthogonality, and, consequently, isinvoluntary (P2 = I).
The application of P to a vector yieldsFigure 18.1 illustrates this formula and makes it clear why P is sometimescalled a Householder reflector: it reflects x about the hyperplaneHouseholder matrices are powerful tools for introducing zeros into vectors.Consider the question “given x and y can we find a Householder matrix P suchthat Px = y?” Since P is orthogonal we clearly require that ||x|| 2 = ||y||2 .Nowand this last equation has the form aυ = x – y for some a. But P is independent of the scaling of v, so we can set a = 1.With υ = x – y we haveand, since xTx = yTy,Thereforeas required.
We conclude that, provided ||x||2 = ||y||2, we can find a Householder matrix P such that Px = y. (Strictly speaking, we have to exclude thecase x = y, which would require υ = 0, making P undefined).18.2 QR FACTORIZATION363Figure 18.1. Householder matrix P times vector x.Normally we choose y to have a special pattern of zeros. The usual choiceis y = σe1 where σ = ±||x||2, which yields the maximum number of zeros iny.
ThenWe choose sign(σ) = –sign(x1) to avoid cancellation in the expression for υ 1 .18.2. QR FactorizationA QR factorization ofwith m > n is a factorizationwhereis orthogonal andis upper triangular. Thematrix R is called upper trapezoidal, since the term triangular applies only tosquare matrices. Depending on the context, either the full factorization A =QR or the “economy size” version A = Q 1 R 1 can be called a QR factorization.A quick existence proof of the QR factorization is provided by the Choleskyfactorization: if A has full rank and ATA = RTR is a Cholesky factorization,then A = AR–1 · R is a QR factorization. The QR factorization is uniqueif A has full rank and we require R to have positive diagonal elements (A =QD · DR is a QR factorization for any D = diag(±1)).The QR factorization can be computed by premultiplying the given matrix by a suitably chosen sequence of Householder matrices.
The process is364QR F ACTORIZATIONillustrated for a generic 4 × 3 matrix as follows:The general process is adequately described by the k th stage of the reduction to triangular form. With A1 = A we have, at the start of the kthstage,where Rk-1 is upper triangular. Choose a Householder matrixand embedinto an m × m matrixsuch that(18.2)Then let Ak+1 = PkAk. Overall, we obtain R = PnPn-1 .
. . P1 A =: QTA(Pn = I if m = n).To compute Ak+1 we need to formWe can writewhich shows that the matrix product can be formed as a matrix–vector product followed by an outer product. This approach is much more efficient thanformingexplicitly and doing a matrix multiplication.The overall cost of the Householder reduction to triangular form is 2n 2 (m–n/3) flops. The explicit formation of Q requires a further 4(m 2 n–mn2 + n 3/3)flops, but for many applications it suffices to leave Q in factored form.18.3. Error Analysis of Householder ComputationsIt is well known that computations with Householder matrices are very stable.
Wilkinson showed that the computation of a Householder vector, andthe application of a Householder matrix to a given matrix, are both normwisestable, in the sense that the computed Householder vector is very close to18.3 E RROR A NALYSISOFH OUSEHOLDER C OMPUTATIONS365the exact one and the computed update is the exact update of a tiny normwise perturbation of the original matrix [1089, 1965, pp. 153–162, 236], [1090,19 6 5 ]. Wilkinson also showed that the Householder QR factorization algorithm is normwise backward stable [1089, p.
236]. In this section we give acombined componentwise and normwise error analysis of Householder matrixcomputations. The componentwise bounds provide extra information overthe normwise ones that is essential in certain applications (for example, theanalysis of iterative refinement).Lemma 18.1. LetConsider the following construction ofandsuch that Px = σe1, where P = I – βvvT is a Householder matrixwith β = 2/(vTv):In floating point arithmetic the computedandandsatisfywhere |θk| < γk .Proof. We sketch the proof. Each occurrence of δ denotes a differentnumber bounded by |δ| < u. We compute fl(xTx) = (1 + θ n ) x T x, and then(the latter term1 + θn +1 is suboptimal, but our main aim is to keep the analysis simple).HenceFor notational convenience, define w = υ 1 + s. We have(essentially because there is no cancellation in the sum).HenceFor convenience we will henceforth write Householder matrices in the formI – v v T, which requires ||υ||2 =and amounts to redefining υ :=QR F ACTORIZATION366and β := 1 in the representation of Lemma 18.1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
