Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 65
Текст из файла (страница 65)
Since these bounds involve A–1) they are nontrivialto compute. Some iterative methods automatically produce estimates of theextremal eigenvalues, and hence ofFor large,sparse symmetric positive definite matrices ||A- 1 ||2 can be cheaply estimatedusing the Lanczos method. Another possibility is to use condition estimationtechniques (Chapter 14).The discussion in this section has inevitably been very general. Otherconsiderations in practice include detecting nonconvergence of a method (dueto rounding errors or otherwise), adapting the tests to accommodate a preconditioned (the residual r provided by the method may now be that for thepreconditioned system), and using a norm that corresponds to a quantity being minimized by the method (a norm that may be nontrivial to compute).16.6 N OTESANDR EFERENCES34316.6. Notes and ReferencesThe Gauss–Seidel method was chosen by Wilkinson [1080, 1948] as an “example of coding” for the ACE. Speaking of experience at that time at theNational Physical Laboratory, he explained that “In general, direct methodshave been used on those equations which did not yield readily to relaxation,and hence those solved by direct methods have nearly always been of an illconditioned type”.Stationary iterative methods are relatively easy to program, although thereare many issues to consider when complicated data structures or parallel machines are used.
A good source of straightforward C, Fortran, and MATLABimplementations of the Jacobi, Gauss–Seidel, and SOR methods, and othernonstationary iterative methods, is the book Templates for the Solution ofLinear Systems [70, 1994]; the software is available from netlib. The bookcontains theoretical discussions of the methods, together with practical advice on topics such as data structures and the choice of a stopping criterion.The choice of stopping criteria for iterative methods is also discussed by Arioli,Duff, and Ruiz [28, 1992].An up-to-date textbook on iterative methods is Axelsson [34, 1994].Problems16.1.
Show that ifand p(B) < 1, then the seriesare both convergent, where ||·|| is any norm.16.2. (Descloux [304,where1963])andConsider the (nonlinear) iterative processsatisfies(16.34)for someG(a) = a.and where ||ek|| < a for all k. Note that a must satisfy(a) Show that(b) Show that the sequence {x k} is bounded and its points of accumulationx satisfy(c) Explain the practical relevance of the result of (b).Previous Home NextChapter 17Matrix PowersUnfortunately, the roundoff errors in the mth power of a matrix, say Bm,are usually small relative to ||B||m rather than ||Bm||.— CLEVE B. MOLER and CHARLES F. VAN LOAN,Nineteen Dubious Ways to Compute the Exponential of a Matrix (1978)It is the size of the hump that matters:the behavior of ||p(A∆ t)n|| = ||p(A∆ t)t/∆ t|| for small but nonzero t.The eigenvalues and the norm, by contrast, give sharp informationonly about the limits tor t0.— DESMOND J.
HIGHAM and LLOYD N. TREFETHEN,Stiffness of ODES (1993)345M ATRIX P OWERS346Powers of matrices occur in many areas of numerical analysis. One approachto proving convergence of multistep methods for solving differential equationsis to show that a certain parameter-dependent matrix is uniformly “powerbounded” [493, 1991 , §V.7], [862, 1992 ]. Stationary iterative methods forsolving linear equations converge precisely when the powers of the iterationmatrix converge to zero. And the power method for computing the largesteigenvalue of a matrix computes the action of powers of the matrix on a vector.It is therefore important to understand the behaviour of matrix powers, inboth exact and finite precision arithmetic.It is well known that the powers Ak oftend to zero as kifp (A) < 1, where p is the spectral radius.
However, this simple statement doesnot tell the whole story. Figure 17.1 plots the 2-norms of the first 30 powers ofa certain 3 × 3 matrix with p(A) = 0.75. The powers do eventually decay, butinitially they grow rapidly. (In this and other similar plots, k on the x -axisis plotted against ||fl(Ak)||2 on the y-axis, and the norm values are joinedby straight lines for plotting purposes.) Figure 17.2 plots the 2-norms of thefirst 250 powers of a 14 × 14 nilpotent matrix C 14 discussed by Trefethenand Trummer [1020, 198 7] (see §17.2 for details). The plot illustrates thestatement of these authors that the matrix is not power bounded in floatingpoint arithmetic, even though its 14th power should be zero.These examples suggest two important questions.• For a matrix with p(A) < 1, how does the sequence {||Ak||} behave? Inparticular, what is the size of the “hump” max k||Ak||?• What conditions on A ensure that fl(Ak)0 as kWe examine these questions in the next two sections.17.1.
Matrix Powers in Exact ArithmeticIn exact arithmetic the limiting behaviour of the powers ofisdetermined by A’s eigenvalues. As already noted, if p(A) < 1 then Akas kif p(A) > 1, ||Ak||asIn the remaining case ofp(A) = 1, ||A k ||if A has a defective eigenvaluesuch thatAk does not converge if A has a nondefective eigenvaluesuch that(although the norms of the powers may converge); otherwise, the onlyeigenvalue of modulus 1 is the nondefective eigenvalue 1 and Ak convergesto a nonzero matrix.
These statements are easily proved using the Jordancanonical form(17.1a)17.1 M ATRIX P OWERSINE XACT A RITHMETICFigure 17.1. A typical hump for a convergent, nonnormal matrix.Figure 17.2. Diverging powers of a nilpotent matrix, C 1 4 .347348MATRIX POWERSwhere X is nonsingular and(17.1b)aswhere n 1 + n2 + . .
. + ns = n. We will call a matrix for which Akk(or equivalently, p(A) < 1) a convergent matrix.The norm of a convergent matrix can be arbitrarily large, as is showntrivially by the scaled Jordan block(17.2)withand a >> 1. While the spectral radius determines the asymptoticrate of growth of matrix powers, the norm influences the initial behaviour ofthe powers. The interesting result that p(A) =for any norm(see Horn and Johnson [580, 1985, p.
299], for example) confirms the asymptotic role of the spectral radius. This formula for p(A) has actually been considered as a means for computing it; see Wilkinson [1089, 1965, pp. 615–617]and Friedland [408, 1991].An important quantity is the “hump” max k ||Ak||, which can be arbitrarilylarge for a convergent matrix.
Figure 17.1 shows the hump for the 3 × 3upper triangular matrix with diagonal entries 3/4 and off-diagonal entries 2;this matrix has 2-norm 3.57. The shape of the plot is typical of that for aconvergent matrix with norm bigger than 1. Note that if A is normal (sothat in (17.1a) J is diagonal and X can be taken to be unitary) we have||Ak|| 2 = ||||2 == p(A)k , so the problem of bounding ||Ak||is of interest only for nonnormal matrices. The hump phenomenon arises invarious areas of numerical analysis. For example, it is discussed for matrixpowers in the context of stiff differential equations by D.
J. Higham andTrefethen [529, 1993], and by Moler and Van Loan [775, 1978] for the matrixexponential eAt with tMore insight into the behaviour of matrix powers can be gained by considering the 2 × 2 matrix (17.2) withand a > 0. We haveand(17.3)17.1 MATRIX POWERSINEXACT ARITHMETIC349Figure 17.3. Infinity norms of powers of 2 × 2 matrix J in (17.2), for λ = 0.99 anda = 0 (bottom line) and a = 10 -k, k = 0:3.HenceIt follows that the norms of the powers can increase for arbitrarily many stepsuntil they ultimately decrease. Moreover, because k–1 tends to zero quiteslowly as kthe rate of convergence ofto zero can be muchslower than the convergence ofto zero (see (17.3)) when is close to 1.
Inother words, nontrivial Jordan blocks retard the convergence to zero.For this 2 × 2 matrix, the hump maxkis easily shown to be approximatelywherethis value being attained forFigure 17.3 displays the norms of the first 400 powers of the matrices withand a = 0,0.001,0.01,0.1,1. The size and location of the hump arecomplicated expressions even in this simple case. When we generalize to directsums of larger Jordan blocks and incorporate a similarity transformation,giving (17.1a), the qualitative behaviour of the powers becomes too difficultto describe precisely.In the rest of this section we briefly survey bounds for ||Ak||. First, however, we comment on the condition numberthat appears350MATRIX POWERSin various bounds in this chapter.
The matrix X in the Jordan form (17.1a)is by no means unique [413, 1959, pp. 220–221], [467, 1976]: if A has distincteigenvalues (hence J is diagonal) then X can be replaced by XD , for anynonsingular diagonal D, while if A has repeated eigenvalues then X can bereplaced by XT, where T is a block matrix with block structure conformalwith that of J and which contains some arbitrary upper trapezoidal Toeplitzblocks. We adopt the convention that κ(X) denotes the minimum possiblevalue of κ(X) over all possible choices of X. This minimum value is notknown for general A, and the best we can hope is to obtain a good estimateof it.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
