Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 47

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 47)

The behaviour of iterative refinement for GE is usually summarized as follows: if double precision is used in the computation of r, andA is not too ill conditioned, then the iteration produces a solution correct toworking precision and the rate of convergence depends on the condition number of A. In the next section we give a componentwise analysis of iterativerefinement that confirms this summary and provides some further insight.11.1. Convergence of Iterative RefinementLetbe nonsingular and letbe a computed solution to Ax = b.Define x 1 =and consider the following iterative refinement process: ri =b – Axi (precision ), solve Adi = ri (precision u), xi+1 = xi + di (precisionu), i = 1, 2, .

. . . For traditional iterative refinement, = u 2. Note that in thischapter subscripts specify members of a vector sequence, not vector elements.We henceforth define ri , di , and xi to be the computed quantities (to avoida profusion of hats). The only assumption we will make on the solver is thatthe computed solutionto a system Ay = c satisfies(11.1)Thus the solver need not be LU factorization or even a factorization method.The page or so of analysis that follows is straightforward but tedious.

Thereader is invited to jump straight to (11.4), at least on first reading.11.1 C ONVERGENCEOF233I TERATIVE R EFINEMENTConsider first the computation of ri . There are two stages. First, si =f l( b – Axi ) = b – Axi + ∆si is formed in the (possibly) extended precision(cf. (3.10)), whereso that |∆ si | <Second, the residual is rounded to the working precision: ri = fl(s i) = si + fi ,where |fi | < u|si |.

HenceBy writing xi = x + (xi – x), we obtain the bound(11.2)For the second step we have, by (11.1), (A + ∆A i)di = ri . Now writewhere, since θi :=Hence(11.3)For the last step,Using (11.3) we haveHenceSubstituting the bound for |∆ ri | from (11.2) gives(11.4)Note that234ITERATIVE REFINEMENTAs long as A is not too ill conditioned and the solver is not too unstable, wehavewhich means that the error contracts until we reach a pointat which the gi term becomes significant. The limiting normwise accuracy,that is, the minimum size ofMoreover, iffor someµ, then we can expect to obtain a componentwise relative error of order µu ,that is, miniWe concentrate now on the case where the solver uses LU factorization. Inthe traditional use of iterative refinement,= u 2, and one way to summarizeour findings is as follows.Theorem 11.1 (mixed precision iterative refinement). Let iterative refinement be applied to the nonsingular linear system Ax = b, using LU factorization and with residuals computed in double the working precision.

Letη = u|| |A - 1||L||U| || where L and U are the computed LU factors of A.Then, provided η is sufficiently less than 1, iterative refinement reduces theerror by a factor approximately η at each stage, untilThis theorem is stronger than the standard results in the literature, whichhavein place of η. We can have η <<since η is independentof the row scaling of A (modulo changes in the pivot sequence). For example,ifthen ηcond(A)u, and cond(A) can be arbitrarily smallerthanConsider now the case where= u, which is called fixed precision iterativerefinement. We have an analogue of Theorem 11.1.Theorem 11.2 (fixed precision iterative refinement). Let iterative refinement in fixed precision be applied to the nonsingular linear system Ax = bof order n, using LU factorization.

Let η = u|| |A – 1 ||L||U| || where L andare the computed LU factors of A. Then, provided η is sufficiently less than1, iterative refinement reduces the error by a factor approximately η at eachstage, until2ncond(A, x)u.The key difference between mixed and fixed precision iterative refinementis that in the latter case a relative error of order u is no longer ensured. Butwe do have a relative error bound of order cond(A, x)u.

This is a strongerbound than holds for the original computed solutionfor which we can sayonly that(this bound is obtained by applying Theorems 7.4 and 9.4, or from (11.4) withi = 0!). In fact, a relative error bound of order cond(A, x)u is the best we canpossibly expect if we do not use higher precision, because it corresponds to the11.2 I TERATIVE R EFINEMENT I MPLIES S TABILITY235uncertainty introduced by making componentwise relative perturbations to Aof size u (again, see Theorem 7.4); this level of uncertainty is usually present,because of errors in computing A or in rounding its elements to floating pointform.The gist of this discussion is that iterative refinement is beneficial even ifresiduals are computed only at the working precision. This fact became widelyappreciated only after the publication of Skeel’s 1980 paper [920, 1980].

Onereason for the delayed appreciation may be that comments such as that madeby Forsythe and Moler, “It is absolutely essential that the residuals rk becomputed with a higher precision than that of the rest of the computation”[396, 1967, p. 49], were incorrectly read to mean that without the use of higherprecision no advantage at all could be obtained from iterative refinement. Inthe next section we will see that fixed precision iterative refinement does morethan just produce a cond(A, x)u-bounded forward error for LU factorization—it brings componentwise backward stability as well.11.2.

Iterative Refinement Implies StabilityWe saw in the last section that fixed precision iterative refinement can improvethe accuracy of a solution computed by GE. The question arises of what therefinement process does to the backward error. To answer this question we givea general backward error analysis that is applicable to a wide class of linearequation solvers.

Throughout this section, “iterative refinement” means fixedprecision iterative refinement.We assume that the computed solutionto Ax = b satisfies(11.5)where g :and h :have nonnegativeentries. The functions g and h may depend on n and u as well as on the dataA and b. We also assume that the residual r = b – A is computed in such away that(11.6)where t :way, then we can takeis nonnegative. If r is computed in the conventional(11.7)First we give an asymptotic result that does not make any further assumptions on the linear equation solver.Theorem 11.3.

Letbe nonsingular. Suppose the linear systemAx = b is solved in floating point arithmetic using a solver S together with oneI TERATIVE R EFINEMENT236step of iterative refinement. Assume that the computed solutionproduced byS satisfies (11.5) and that the computed residual? satisfies (11.6). Then thecorrected solutionsatisfies(11.8)where q = O(u) ifProof. The residual r = b – Asatisfiesof the original computed solution(11.9)The computed residual iscorrection d satisfies= r + ∆r , where |∆ r| <The computedFinally, for the corrected solution we haveCollecting together the above results we obtainHence(11.12)whereThe claim about the order of q follows sinceorder u.andare all ofTheorem 11.3 shows that, to first order, the componentwise relative backward error w|A|,|b| will be small after one step of iterative refinement as long asandare bounded by a modest scalar multiple ofThis is true for t if the residual is computed in the conventional way (see(11.7)), and in some cases we may take h0, as shown below.

Note that thefunction g of (11.5) does not appear in the first-order term of (11.8). Thisis the essential reason why iterative refinement improves stability: potentialinstability manifested in g is suppressed by the refinement stage.A weakness of Theorem 11.3 is that the bound (11.8) is asymptotic. Sincea strict bound for q is not given, it is difficult to draw firm conclusions about11.2 I TERATIVE R EFINEMENT I MPLIES S TABILITY237the size of w|A|,|b|. The next result overcomes this drawback, at the cost ofsome specialization (and a rather long proof).We introduce a measure of ill scaling of the vector |B||x|,Theorem 11.4. Under the conditions of Theorem 11.3, suppose that g(A, b) =G|A| and h(A, b) = H|b|, where G,have nonnegative entries, andthat the residual is computed in the conventional manner.

Then there is afunctionsuch that ifthenProof. As with the analysis in the previous section, this proof can beskipped without any real loss of understanding. From (11.12) in the proof ofTheorem 11.3, using the formula (11.7) for t, we have(11.13)The inequality (11.9) impliesIf< 1/2 (say) then I – uH isor (I – uH)|b| < (I + uG)nonsingular with a nonnegative inverse satisfying ||(I – uH) –1 || < 2 and wecan solve for |b| to obtain |b| < (I – uH)-1 (I + uG) |A|It follows from thisrelation and consideration of the rest of the proof that the simplifying stepof replacing b by 0 in the analysis has little effect on the bounds—it merelyproduces unimportant perturbations in f in the statement of the theorem.Making this replacement in (11.13) and approximating γ n +1 + uγ n+1, wehave(11.14)Our task is now to boundandulating (11.11) we obtain the inequalityin terms ofBy manip-(11.15)ITERATIVE REFINEMENT238Also, we can boundbyand dropping the |b| terms and using (11.15) givesSubstituting from (11.15) and (11.16) into (11.14) we findwhereNow from (11.10), making use of (11.16),After premultiplying by |A| this may be rearranged as(11.18)whereUsing γn + 1 /u < (n + 1)/(1 – (n + 1)u)n + 1, we have the boundsIf< 1/2 (say) then (I – uM3)-1 > 0 with ||(I – uM3)-1||we can rewrite (11.18) as< 2 and(11.19)11.2 I TERATIVE R EFINEMENT I MPLIES S TABILITY239Substituting this bound into (11.

Характеристики

Список файлов книги