Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 64
Текст из файла (страница 64)
Because of the structure and the possibly largedimension of the coefficient matrices in these applications, iterative solutionmethods are frequently used. An important question is how the rather delicate convergence properties of the iterative methods are affected by roundingerrors. In this section we extend the analysis of stationary iteration to singularsystems.16.4.1. Theoretical BackgroundA useful tool in analysing the behaviour of stationary iteration for a singularsystem is the Drazin inverse. This can be defined, foras theunique matrix AD such thatwhere k = index(A).
The index of A is the smallest nonnegative integer ksuch that rank(Ak) = rank(A k+1); it is characterized as the dimension of thelargest Jordan block of A with eigenvalue zero. If index(A) = 1 then AD33716.4 S INGULAR S YSTEMSis also known as the group inverse of A and is denoted by A#. The Drazininverse is an “equation-solving inverse” precisely when index(A) < 1, for thenAADA = A, and so if Ax = b is a consistent system then ADb is a solution.As we will see, however, the Drazin inverse of the coefficient matrix A itselfplays no role in the analysis.
The Drazin inverse can be represented explicitlyas follows. Ifwhere P and B are nonsingular and N has only zero eigenvalues, thenFurther details of the Drazin inverse can be found in Campbell and Meyer’sexcellent treatise [180, 1979, Chap. 7].Letbe a singular matrix and consider solving Ax = b bystationary iteration with a splitting A = M – N, where M is nonsingular.First, we examine the convergence of the iteration in exact arithmetic. Sinceany limit point x of the sequence {x k } must satisfy Mx = Nx + b, or Ax = b,we restrict our attention to consistent linear systems.
(For a thorough analysisof stationary iteration for inconsistent systems see Dax [269, 1990]. ) As inthe nonsingular case we have the relation (cf. (16.4)):(16.21)mwhere G = M – 1N. Since A is singular, G has an eigenvalue 1, so G doesnot tend to zero asthat is, G is not convergent. If the iterationmust exist. Following Meyer andis to converge for all x0 thenPlemmons [752, 1977], we call a matrix B for whichexists semiconvergent.We assume from this point on that G is semiconvergent. It is easy to see[94, 1994, Lem. 6.9] that G must have the form(16.22)where P is nonsingular and< 1. HenceTo rewrite this limit in terms of G, we note that(16.23)STATIONARY ITERATIVE METHODS338and, sinceis nonsingular,(16.24)Hence(16.25)To evaluate the limit of the second term in (16.21) we note that, since thesystem is consistent, M– 1 b = M-1 Ax = (I – G)x, and soWe note in passing that the condition that G is semiconvergent is equivalentto I – G having index 1, in view of (16.23), but that this condition does notimply that A = M(I – G) has index 1.The conclusion is that if G is semiconvergent, stationary iteration converges to a solution of AX = b that depends on x 0 :(16.26)The first term in this limit is in null(I – G) and the second term is in range(I –G).
To obtain the unique solution in range(I – G) we should take for x0. anyvector in range(I – G) (x0 = 0, say).16.4.2. Forward Error AnalysisWe wish to bound em+1 = x –where x is the limit in (16.26) corresponding to the given starting vector x0. The analysis proceeds as in thenonsingular case, up to the derivation of equation (16.5):As before, the first term, Gm+ 1e 0, is negligible for large m, because itis the error after m + 1 stages of the exact iteration and this error tendsto zero. To obtain a useful bound for the second term, we cannot simplytake norms or absolute values, becausegrows unfoundedly with m(recall that G has an eigenvalue 1).
Our approach is to split the vectorsaccording towhereand16.4 S INGULAR SYSTEMS339null(I – G); this is a well-defined splitting because range(I – G) and null(I – G)are complementary subspaces (since index(I – G) = 1, or equivalently, G issemiconvergent). Using the properties of the splitting the error can be writtenasWe achieve the required splitting for ξi via the formulaewhereHence the error can be written as(16.27)Clearly, as mthe final term in this expression can become unbounded,but since it grows only linearly in the number of iterations it is unlikely tohave a significant effect in applications where stationary iteration convergesquickly enough to be of practical use.Now we bound the term(16.28)Using inequality (16.2) and the definition of γx in (16.7) and θx in (16.9), wehave(16.29)The convergence of the two infinite sums is assured by the result of Problem 16.1, since by (16.22)-(16.24),G i E = Gi (I – G)D (I – G)(16.30)340S TATIONARY I TERATIVE M ETHODSwhereWe conclude that we have the normwise error bound(16.31)On setting E = I we recover the result (16.8) for the nonsingular case.
If weassume that Γ is diagonal, so that P in (16.30) is a matrix of eigenvectors ofG, thenThis bound shows that a small forward error is guaranteed ifO(1) and the second largest eigenvalue of G is not too close to 1. (It is thissubdominant eigenvalue that determines the asymptotic rate of convergenceof the iteration.)Turning to the componentwise case, we see from (16.24) and (16.30) thatBecause of the form of the sum in (16.29), this prompts us to define the scalarc(A) > 1 byin terms of which we have the componentwise error bound(16.32)Again, as a special case we have the result for nonsingular A, (16.13).To what should we compare this bound? A perturbation result for Ax = bis given in [564, 1993] that projects the perturbations of A and b into range (I –G) and thus can be thought of as gauging the effect of perturbations to the“nonsingular part of the system”.
For perturbations of orderit gives anexpressionHence we can deduce conditions on a stationary iterative method that ensureit is componentwise forward stable, in the sense of yielding a solution whose16.5 S TOPPINGANI TERATIVE M ETHOD341error is no larger than the uncertainty in x caused by rounding the data. Theconstants θx and c(A) should be bounded by d n, where dn denotes a slowlyshould hold, as itgrowing function of n; the inequality |M| + |N| <does for the Jacobi method and for the SOR method whenwhereβ is positive and not too close to zero; and the “exact error” Gm + 1e 0 mustdecay quickly enough to ensure that the term (m + 1) | (I – E) M-1| does notgrow too large before the iteration is terminated.Numerical results given in [564, 1993] show that the analysis can correctlypredict forward and backward stability, and that for certain problems lineargrowth of the component of the error in null(A) can indeed cause an otherwiseconvergent iteration to diverge, even when starting very close to a solution.16.5.
Stopping an Iterative MethodWhat convergence test should be used to stop an iterative linear equationsolver? In this section we explain how backward errors and condition numbershelp answer this question. Note first that most iterative methods for solvingAx = b compute all or part of a matrix–vector product w = Aυ on eachiteration, and in floating point arithmetic we havewhere m is the maximum number of nonzeros per row of A. The methodtherefore cannot distinguish between A and A + ∆A where |∆ A| < γm |A|, andso there is no point in trying to achieve a componentwise relative backwarderror less than γm .
Of course, instability of a method (or simply lack ofconvergence) may pose further restrictions on how small a backward errorcan be achieved.It is worth keeping in mind throughout this discussion that in practicalapplications accuracy and stability requirements are often quite modest because of large errors or uncertainties in the data, or because the iteration isan “inner step” of some “outer” iterative process. Indeed, one of the advantages of iterative methods is that the slacker the convergence tolerance theless computational effort is required, though the relation between toleranceand work depends very much on the method.Natural stopping criteria for an iterative method are that some measureof backward error or forward error does not exceed a tolerance,We willassume that the residual r = b – Ay is available for each iterate y , and thatnorms of y, r, and A can be computed or estimated.
If r is not computeddirectly, but is recurred by the method, as, for example, in the conjugategradient method, then the norm of the computed residual may differ fromthat of the true residual by several orders of magnitude; clearly, this affectsthe way that the stopping tests are interpreted.342S TATIONARY I TERATIVE M ETHODSFrom Theorem 7.1 we have the following equivalences, for any subordinatematrix norm:(16.33a)(16.33b)(16.33c)These inequalities remain true with norms replaced by absolute values (Theorem 7.3), but to evaluate (16.33b) and (16.33c) a matrix–vector product |A||y|must be computed, which is a nontrivial expense in an iterative method.Of these tests, (16.33c) is preferred in general, assuming it is acceptableto perturb both A and b. Note the importance of including both ||A|| and||y|| in the test on ||r||; a test ||r|| <though scale independent, doesnot bound any relative backward error.
Test (16.33a) is commonly used inexisting codes, but may be very stringent, and possibly unsatisfiable. To seewhy, note that the residual of the rounded exact solution fl(x) = x + ∆ x ,|∆x| < u|x|, satisfies, for any absolute norm,andIf A is ill conditioned and x is a large-normed solution (that is, ||x||||A - 1 ||||b||), so that ||b|| is close to its lower bound, then (16.33a) is muchharder to satisfy than (16.33c).If the forward error is to be bounded, then, for a nonsingular problem, testscan be derived involving the residual and A–1: the equality x – y = A – 1 rleads to normwise and componentwise forward error bounds, such as ||x –y||/||y|| < ||A- 1 ||||r||/||y||.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
