Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 62
Текст из файла (страница 62)
The basic question that our analysis attempts to answeris “What is the limiting accuracy of a method in floating point arithmetic?”Specifically, how small can we guarantee that the backward or forward errorwill be over all iterations k = 1, 2 ,. . .? Without an answer to this questionwe cannot be sure that a convergence test of the form(say)will ever be satisfied, for any given value ofAs an indication of the potentially devastating effects of rounding errorswe present an example constructed and discussed by Hammarling and Wilkinson [497, 1976]. Here, A is the 100 × 100 lower bidiagonal matrix with aii1.5and ai,i– 11, and bi2.5.
The successive overrelaxation (SOR) methodis applied in MATLAB with parameter w = 1.5, starting with the roundedversion of the exact solution x, given by xi = 1 – (–2/3) i . The forward errorsand the co-norm backward errorsare plotted inFigure 16.1. The SOR method converges in exact arithmetic, since the iteration matrix has spectral radius 1/2, but in the presence of rounding errors itdiverges. The iteratehas a largest element of order 1013,for16.1 S URVEYOFERROR ANALYSIS327Figure 16.1.
SOR iteration.k > 238, and for k > 100,The divergenceis not a result of ill conditioning of A, sinceThe reason for theinitial rapid growth of the errors in this example is that the iteration matrixis far from normal; this allows the norms of its powers to become very largebefore they ultimately decay by a factor 1/2 with each successive power.The effect of rounding errors is to cause the forward error curve in Figure 16.1to level off near k = 100, instead of decaying to zero as it would in exact arithmetic.
More insight into the initial behaviour of the errors can be obtainedusing the notion of pseudo-eigenvalues; see §17.3.16.1. Survey of Error AnalysisBefore analysing stationary iterative methods, we briefly survey the publishederror analysis for iterative methods.
For symmetric positive definite systems,Golub [468, 196 2] derives both statistical and nonstatistical bounds for theforward error and residual of the Richardson method. Benschop and Ratz [92,1971] give a statistical analysis of the effect of rounding errors on stationaryiteration, under the assumption that the rounding errors are independentrandom variables with zero mean. Lynn [719, 196 4] presents a statisticalanalysis for the SOR method with a symmetric positive definite matrix.Hammarling and Wilkinson [497, 1976] give a normwise error analysis for328S TATIONARY I TERATIVE M ETHODSthe SOR method. With the aid of numerical examples, they emphasize thatwhile it is the spectral radius of the iteration matrix M – 1 N that determinesthe asymptotic rate of convergence, it is the norms of the powers of this matrixthat govern the behaviour of the iteration in the early stages.
This point isalso explained by Trefethen [1017, 1992], using the tool of pseudospectra.Dennis and Walker [302, 198 4] obtain bounds forfor stationary iteration as a special case of error analysis of quasi-Newtonmethods for nonlinear systems. The bounds in [302, 198 4] do not readilyyield information about normwise or componentwise forward stability.Bollen [133, 1984] analyses the class of “descent methods” for solving Ax =b, where A is required to be symmetric positive definite; these are obtained byiteratively using exact line searches to minimize the quadratic function F(x) =(A-1 b – x)TA(A-1 b – x). The choice of search direction pk = b – Axk =: rkyields the steepest descent method, while pk = ej (unit vector), where |rk|j =gives the Gauss-Southwell method.
Bollen shows that both methodsare normwise backward stable as long as a condition of the form cn κ(A)u < 1holds. If the pk are cyclically chosen to be the unit vectors e 1, e2, . . ., en thenthe Gauss–Seidel method results, but unfortunately no results specific to thismethod are given in [133, 1984].Wozniakowski [1112, 1977] shows that the Chebyshev semi-iterative methodis normwise forward stable but not normwise backward stable, and in [1113,1978] he gives a normwise error analysis of stationary iterative methods.
Someof the assumptions in [1113, 1978] are difficult to justify, as explained byHigham and Knight [563, 1993].In [1114, 198 0] Wozniakowski analyses a class of conjugate gradient algorithms (which does not include the usual conjugate gradient method). Heobtains a forward error bound proportional to κ(A)3/2 and a residual boundproportional to K(A), from which neither backward nor forward normwise stability can be deduced. We note that as part of the analysis in [1114, 198 0]Wozniakowski obtains a residual bound for the steepest descent method thatis proportional to K(A), and is therefore much weaker than the bound obtainedby Bollen [133, 1984].Zawilski [1125, 1991] shows that the cyclic Richardson method for symmetric positive definite systems is normwise forward stable provided the parameters are suitably ordered.
He also derives a sharp bound for the residualthat includes a factor κ(A), and which therefore shows that the method is notnormwise backward stable.Arioli and Romani [29, 1992] give a statistical error analysis of stationaryiterative methods. They investigate the relations between a statistically defined asymptotic stability factor, ill conditioning of M- 1 A, where A = M – Nis the splitting, and the rate of convergence.Greenbaum [479, 1989] presents a detailed error analysis of the conjugategradient method, but her concern is with the rate of convergence rather than16.2 FORWARD ERROR ANALYSIS329the attainable accuracy. An excellent survey of work concerned with the effectsof rounding error on the conjugate gradient method (and the Lanczos method)is given by Greenbaum and Strakos in the introduction of [480, 1992]; see alsoGreenbaum [481, 1994].
Notay [799, 1993] analyses how rounding errors influence the convergence rate of the conjugate gradient method for matrices withisolated eigenvalues at the ends of the spectrum. Van der Vorst [1058, 1990]examines the effect of rounding errors on preconditioned conjugate gradientmethods with incomplete Cholesky preconditioners.The analysis given in the remainder of this chapter is from Higham andKnight [563, 1993], [564, 1993], wherein more details are given.
Error analysisof Kaczmarz’s row-action method is given by Knight [663, 1993].16.2. Forward Error AnalysisA stationary iterative method has the formwhereis nonsingular and M is nonsingular. We assumethat the spectral radius p(M - 1 N) < 1, so that in exact arithmetic the iteration converges for any starting vector x0. We are not concerned with the sizeof constants in this analysis, so we denote by cn a constant of order n.The computed vectorssatisfy an equality of the formwhich we write as(16.1)whereWe will assume that M is triangular (as is the case for the Jacobi, Gauss–Seidel, SOR, and Richardson iterations), so thatand f kaccounts solely for the errors in formingHence(16.2)Solving the recurrence (16.1) we obtain(16.3)where G = M– 1 N.
Since the iteration is stationary at x,(16.4)S TATIONARY I TERATIVE M ETHODS330and so the error em+1 :=satisfies(16.5)We have(16.6)where µ k is the bound for ξk defined in (16.2). The first term, |Gm+ 1e0 |, isthe error of the iteration in exact arithmetic and is negligible for large m. Theaccuracy that can be guaranteed by the analysis is therefore determined bythe last term in (16.6), and it is this term on which the rest of the analysisfocuses.At this point we can proceed by using further componentwise inequalitiesor by using norms. First we consider the norm approach.
By taking norms in(16.6) and defining(16.7)we obtain(16.8)where the existence of the sum is assured by the result of Problem 16.1.If= q < 1 then (16.8) yieldsThus if q is not too close to 1 (q < 0.9, say), and γx andare not toolarge, a small forward error is guaranteed for sufficiently large m.Of more interest is the following componentwise development of (16.6).Defining(16.9)so thatfor all k, we have from (16.2),(16.10)16.2 FORWARD ERROR ANALYSIS331Hence (16.6) yields(16.11)where, again, the existence of the sum is assured by the result of Problem 16.1.Since A = M – N = M(I – M- 1 N) we haveThe sum in (16.11) is clearly an upper bound for |A - 1 |. Defining c(A) > 1 by(16.12)we have our final bound(16.13)An interesting feature of stationary iteration methods is that if the elements of M and N are multiples of the elements in the corresponding positionsof A, then any scaling of the formdiagonal) leaves the eigenvalues of M– 1 N unchanged; hence the asymptoticconvergence rate is independent of row and column scaling.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
