Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 38
Текст из файла (страница 38)
Fori = 1 the result is immediate, and if it is true for i - 1 then, from (9.16).as required. Note that, similarly, |ui | < |di | + |ci |. Finally,Theorem 9.13. If the nonsingular tridiagonal matrix A is of type (a)-(d) inTheorem 9.11, and if the unit roundoff u is sufficiently small, then GE forsolving Ax = b succeeds and the computed solutionsatisfiesThe same conclusion is true if A is diagonally dominant by rows or columns,with no restriction on u, provided the bound is multiplied by 3.186LU F ACTORIZATIONANDLINEAR E QUATIONSProof. If u is sufficiently small then for types (a)-(c) the diagonal elementsof U will be positive, sinceas u 0.
It is easy to see thatfor all i ensures thatThe argument is similar for type (d). Theresult therefore follows from (9.19) and (9.8). The last part is trivial.A corollary of Theorem 9.13 is that it is not necessary to pivot for thematrices specified in the theorem (and, indeed, pivoting could vitiate the resultof the theorem). Note that large multipliers may occur under the conditionsof the theorem, but they do not affect the stability. For example, consider theM-matrixwhereThe multiplier l32 is unbounded asbut |L||U| = |A|and GE performs very stably, as Theorem 9.13 shows it must.9.6. Historical PerspectiveGE was the first numerical algorithm to be subjected to rounding error analysis, so it is instructive to follow the development of the error analysis fromits beginnings in the 1940s.In the 1940s there were three major papers giving error analyses of GE.Hotelling [583, 1943] presented a short forward error analysis of the LU factorization stage of GE.
Under the assumptions that |aij| < 1 and |bi | < 1 for alli and j and that the pivots are all of modulus unity, Hotelling derives a boundcontaining a fact or 4 n -1 for the error in the elements of the reduced uppertriangular system.
Hotelling’s work was quoted approvingly by Bargmann.Montgomery, and von Newmann [55, 1946], who dismiss elimination met hodsfor the solution of a linear system Ax = b as being numerically unstable. Instead, they recommended computation of A-1 via the Newton Schulz iteration[908, 1933] (which was also discussed by Hotelling). In one paragraph theyout line the alleged shortcomings of elimination methods as follows:In the elimination method a series of n compound operations isperformed each of which depends on the preceding. An error atany stage affects all succeeding results and may become greatlymagnified; this explains roughly why instability should be expected.
It should be noticed that at each step a division is performed by a number whose size cannot be estimated in advanceand which might the so small that any error in it would be greatlymagnified by division. In fact such small divisors must occur if thedeterminant of the matrix is small and may occur even if it is not9.6 H I S T O R I C A L P E R S P E C T I V E187. . . Another reason to expect instability is that once the variablexn is obtained all the other variables are expressed in terms of it.As Wilkinson [1098, 1974, p. 354] notes of this paragraph, “almost every statement in it is either wrong or misleading”.Hotelling’s result led to general pessimism about the practical effectivenessof GE for solving large systems of equations. Three papers later in the samedecade helped to restore confidence in GE.Goldstine [460, 1972, p.
290] says of his discussions with von Neumann:We did not feel it reasonable that so skilled a computer as Gausswould have fallen into the trap that Hotelling thought he had noted. . . Von Neumann remarked one day that even though errors maybuild up during one part of the computation, it was only relevantto ask how effective is the numerically obtained solution, not howclose were some of the auxiliary numbers, calculated on the wayto their correct counterparts.
We sensed that at least for positivedefinite matrices the Gaussian procedure could be shown to bequite stable.von Neumann and Goldstine [1057, 1947] subsequently gave a long and difficult rigorous fixed-point error analysis for the inversion of a symmetric positive definite matrix A via GE. They showed that the computed inversesatisfies< 14.2n 2 uκ2 (A). Parlett [821, 1990] explains that “the joyof this result was getting a polynomial in n, and the pain was obtaining 14.2,a number that reflects little more than the exigencies of the analysis.” Wilkinson [1095, 1971] gives an interesting critique of von Neumann and Goldstine’spaper and points out that the residual bound could hardly be improved usingmodern error analysis techniques. In a later paper [462, 1951], Goldstine andvon Neumann gave a probabilistic analysis, which Goldstine summarizes asshowing that “under reasonable probabilistic assumptions the error estimatesof the previous paper could be reduced from a proportionality of n 2 to n”[460, 1972, p.
291].In his 1970 Turing Award Lecture [1096, 1971], Wilkinson recounts how inthe early 1940s he solved a system of 12 linear equations on a desk calculator,obtaining a small residual. He goes on to describe a later experience:It happened that some time after my arrival [at the National Physical Laboratory in 1946], a system of 18 equations arrived in Mathematics Division and after talking around it for some time we finallydecided to abandon theorizing and to solve it .
. . The operationwas manned by Fox, Goodwin, Turing, and me, and we decidedon Gaussian elimination with complete pivoting . . . Again the system was mildly ill-conditioned, the last equation had a coefficientof order 10-4 (the original coefficients being of order unity) and188LU FACTORIZATIONANDLINEAR E QUATIONSthe residuals were again of order 10-10 , that is of the size corresponding to the exact solution rounded to ten decimals. It isinteresting that in connection with this example we subsequentlyperformed one or two steps of what would now be called “iterativerefinement,” and this convinced us that the first solution had hadalmost six correct figures.(Fox [403, 19 8 7 ] notes that the computation referred to in this quotationtook about two weeks using desk computing equipment!) In a subsequentpaper, Fox, Huskey, and Wilkinson [404, 1948] presented empirical evidencein support of GE, commenting that “in our practical experience on matricesof orders up to the twentieth, some of them very ill-conditioned, the errorswere in fact quite small”.The experiences of Fox, Huskey, and Wilkinson prompted Turing to writea remarkable paper “Rounding-off errors in matrix processes” [1027, 1948].In this paper, Turing made several import ant contributions.
He formulatedthe LU (actually. the LDU) factorization of a matrix, proving the “if” partof Theorem 9.1 and showing that GE computes an LDU factorization. Heintroduced the term “condition number” and defined two matrix conditionnumbers, one of which is n - 1N(A )N( A -1), where N(A) = ||A||F, the “Ncondition number of A”. He used the word “preconditioning” to mean improving the condition of a system of linear equations (a term that did notcome into popular use until the 1970s). He described iterative refinementfor linear systems. He exploited backward error ideas, for example by notingthat *‘the triangular resolution obtained is an exact resolution of a matrixA - S, where M(S) <( M(S) := maxi,j |sij|. Finally, and perhaps mostimportantly, he analysed GEPP for general matrices and obtained a boundforthat contains a term proportional to(By making atrivial change in the analysis, namely replacing A - 1 b by x, Turing’s boundcan be made proportional only toTuring also showed that thefactor 4n -1 in Hotelling’s bound can be improved to 2 n -1 and that still thebound is attained only in exceptional cases.In a review of Turing’s paper, Bodewig [129, 1949] described the errorbounds as “impractical“ and advocated computing the residual of the computed solution and then determining “the exact correct ion by solving a newsystem.” That another researcher could miss the point of Turing’s analysisemphasizes how new the concept of rounding error analysis was in the 1940s.Table 9.1 shows the time for solution of linear systems by GE on someearly computing devices.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
