Nash - Compact Numerical Methods for Computers (523163), страница 21
Текст из файла (страница 21)
The Healy and Price routines mostlikely are slowed down by the manner in which their looping is organised,requiring a backward jump in the program. On a system which interprets, this islikely to require a scan of the program until the appropriate label is found in thecode and it is my opinion that this is the reason for the more rapid execution ofalgorithm 7 in interpreter environments. Such examples as these serve as awarning that programs do show very large changes in their relative, as well asabsolute, performance depending on the hardware and software environment inwhich they are run.Algorithm 7 computes the triangular factor L column by column. Row-by-rowdevelopment is also possible, as are a variety of sequential schemes which performthe central subtraction in STEP 4 in piecemeal fashion. If double-precision arithmetic is possible, the forms of the Choleski decomposition which compute theright-hand side of (7.8) via a single accumulation as in algorithm 7 have theadvantage of incurring a much smaller rounding error.Example 7.1.
The Choleski decomposition of the Moler matrixThe Moler matrix, given byAij = min(i, j) – 2Aii = ifor ijhas the decompositionA = LL Twhere0L ij =1for j > ifor i = j-1for j < i.Note that this is easily verified sincefor i jfor i = j.On a Data General NOVA operating in arithmetic having six hexadecimal digits(that is, a machine precision of 16-5) the correct decomposition was observed forMoler matrices of order 5, 10, 15, 20 and 40. Thus for order 5, the Moler matrix92Compact numerical methods for computershas a lower triangle-1 2-1 0 3-1 0 1 4-1 0 1 2 5which gives the Choleski factor1-1-1-1-11-1-1-11-1 1-1 -11.Example 7.2. Solving least-squares problems via the normal equationsUsing the data in example 3.2, it is straightforward to form the matrix B from thelast four columns of table 3.1 together with a column of ones.
The lower triangleof BT B is then (via a Hewlett-Packard 9830 in 12 decimal digit arithmetic)189268236359705 2 1 6 4 3 7 910985647 3734131 64454373344971 1 1 6 6 5 5 9 2008683 6592261470951478 8 5 9 2926 13Note that, despite the warnings of chapter 5, means have not been subtracted, sincethe program is designed to perform least-squares computations when a constant(column of ones) is not included. This is usually called regression through theorigin in a statistical context.
The Choleski factor L is computed as4350·4968681461·834175165·58938642525·147663258·3731371 48·05831416768·8710282 257·2450797 14·66763457 40·904419643·380993125 1·235276666 0·048499519 0·1948963630·051383414The Choleski decomposition93Using the right-hand side59379382046485TB y = 352641311301775003the forward- and back-substitution algorithm 8 computes a solutionx =-0·0461924351·019386565-0·159822924-0·290376225207·7826146This is to be compared with solution (a) of table 3.2 or the first solution ofexample 4.2 (which is on pp 62 and 63), which shows that the various methods allgive essentially the same solution under the assumption that none of the singularvalues is zero.
This is despite the fact that precautions such as subtracting meanshave been ignored. This is one of the most annoying aspects of numericalcomputation-the foolhardy often get the right answer! To underline, let us usethe above data (that is BT B and BTy ) in the Gauss elimination method, algorithms5 and 6. If a Data General NOVA operating in 23-bit binary arithmetic is used,the largest integer which can be represented exactly is223– 1 = 8388607so that the original matrix of coefficients cannot be represented exactly. However,the solution found by this method, which ignores the symmetry of B T B, isx=-4·62306E-21·01966-0·159942-0·288716207·426While this is not as close as solution (a) of table 3.2 to the solutions computed incomparatively double-length arithmetic on the Hewlett-Packard 9830, it retainsthe character of these solutions and would probably be adequate for manypractitioners. The real advantage of caution in computation is not, in my opinion,that one gets better answers but that the answers obtained are known not to beunnecessarily in error.Previous HomeChapter 8THE SYMMETRIC POSITIVE DEFINITEMATRIX AGAIN8.1.
THE GAUSS-JORDAN REDUCTTONAnother approach to the solution of linear equations (and in fact nonlinear ones)is the method of substitution. Here, one of the unknowns xk is chosen and one ofthe equations in which it has a non-zero coefficient is used to generate anexpression for it in terms of the other xj, j k, and b. For instance, givenAx = band A11(2.2)0, we might use(8.1)The substitution of this into the other equations gives a new set of (n - 1)equations in (n - 1) unknowns which we shall writeA'x' = b'(8.2)in which the indices will run from 2 to n.
In fact x' will consist of the last (n - 1)elements of x. By means of (8.1) it is simple to show that(8.3)and(8.4)for k, j = 2, . . . , n. Notice that if b is included as the (n + 1)th column of A, (8.4) isthe only formula needed, though j must now run to (n + 1).We now have a set of equations of order (n – l), and can continue the processuntil only a set of equations of order 1 remains. Then, using the trivial solution ofthis, a set of substitutions gives the desired solution to the original equations. Thisis entirely equivalent to Gauss elimination (without pivoting) and back-substitution,and all the arithmetic is the same.Consider, however, that the second substitution is made not only of x2 into theremaining (n – 2) equations, but also into the formula (8.1) for x1.
Then the finalorder-1 equations yield all the xj at once. From the viewpoint of elimination itcorresponds to eliminating upper-triangular elements in the matrix R in thesystemRx = f(6.4)94NextThe symmetric positive definite matrix again95then dividing through by diagonal elements. This leaves(8.5)1x = f 'that is, the solution to the set of equations.Yet another way to look at this procedure is as a series of elementary rowoperations (see Wilkinson 1965) designed to replace the pth column of an n by nmatrix A with the p th column of the unit matrix of order n, that is, eP.
Toaccomplish this, the pth row of A must be divided by APP, and Aip times theresulting pth row subtracted from every row i for i p. For this to bepossible, of course, App cannot be zero.A combination of n such steps can be used to solve sets of linear equations. Toavoid build-up of rounding errors, some form of pivoting must be used. This willinvolve one of a variety of schemes for recording pivot positions or else will useexplicit row interchanges.
There are naturally some trade-offs here betweensimplicity of the algorithm and possible efficiency of execution, particularly if theset of equations is presented so that many row exchanges are needed.By using the Gauss-Jordan elimination we avoid the programming needed toperform the back-substitution required in the Gauss elimination method. Theprice we pay for this is that the amount of work rises from roughly n 3 /3operations for a single equation to n 3/2 as n becomes large (see Ralston 1965p 401). For small n the Gauss-Jordan algorithm may be more efficient dependingon factors such as implementation and machine architecture. In particular, it ispossible to arrange to overwrite the ith column of the working matrix with thecorresponding column of the inverse. That is, the substitution operations ofequations (8.1) and (8.4) with 1 replaced by i give elimination formulaeÃij = Aij/AiiÃkj = Akj – Aki (A ij /A ii)(8.1a)(8.4a)for j = 1, 2, .
. . , n, k = 1, 2 , . . . , n, but ki, with the tilde representing thetransformed matrix. However, these operations replace column i with ei, the ithcolumn of the unit matrix 1 n, information which need not be stored. Theright-hand side b is transformed according to(8.1b)(8.3a)for k = 1, 2, . . . , n with k i.
To determine the inverse of a matrix, we could solvethe linear-equation problems for the successive columns of 1n. But now allcolumns ej for j > i will be unaltered by (8.1b) and (8.3a). At the ith stage of thereduction, ei can be substituted on column i of the matrix by storing the pivot Aii,substituting the value of 1 in this diagonal position, then performing the divisionimplied in (8.1a). Row i of the working matrix now contains the multipliersÃij = (Aij/Aii). By performing (8.4a) row-wise, each value Aki can be saved, azero substituted from ei, and the elements of Akj, j = 1, 2, . . .
, n, computed.96Compact numerical methods for computersThis process yields a very compact algorithm for inverting a matrix in theworking storage needed to store only a single matrix. Alas, the lack of pivotingmay be disastrous. The algorithm will not work, for instance, on the matrix0 11 0which is its own inverse. Pivoting causes complications. In this example, interchanging rows to obtain a non-zero pivot implies that the columns of the resultinginverse are also interchanged.The extra work involved in column interchanges which result from partialpivoting is avoided if the matrix is symmetric and positive definite-this specialcase is treated in detail in the next section. In addition, in this case completepivoting becomes diagonal pivoting, which does not disorder the inverse.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
