Nash - Compact Numerical Methods for Computers (523163), страница 22
Текст из файла (страница 22)
Therefore algorithms such as that discussed above are widely used to perform stepwiseregression, where the pivots are chosen according to various criteria other thanerror growth. Typically, since the pivot represents a new independent variableentering the regression, we may choose the variable which most reduces theresidual sum of squares at the current stage.
The particular combination of (say) mout of n independent variables chosen by such a forward selection rule is notnecessarily the combination of m variables of the n available which gives thesmallest residual sum of squares. Furthermore, the use of a sum-of-squares andcross-products matrix is subject to all the criticisms of such approaches toleast-squares problems outlined in chapter 5.As an illustration, consider the problem given in example 7.2.
A Data GeneralECLIPSE operating in six hexadecimal digit arithmetic gave a solutionx=-4·64529E-21·02137-0·160467-0·285955206·734when the pivots were chosen according to the residual sum-of-squares reductioncriterion. The average relative difference of these solution elements from those ofsolution (α ) of table 3.2 is 0·79%. Complete (diagonal) pivoting for the largestelement in the remaining submatrix of the Gauss-Jordan working array gave asolution with an average relative difference (root mean square) of 0·41%.
Thereare, of course, 120 pivot permutations, and the differences measured for eachsolution ranged from 0·10% to 0·79%. Thus pivot ordering does not appear to bea serious difficulty in this example.The operations of the Gauss-Jordan algorithm are also of utility in the solutionof linear and quadratic programming problems as well as in methods derived fromsuch techniques (for example, minimum absolute deviation fitting). Unfortunately,these topics, while extremely interesting, will not be discussed further in thismonograph.97The symmetric positive definite matrix again8.2.
THE GAUSS-JORDAN ALGORITHM FOR THE INVERSEOF A SYMMETRIC POSITIVE DEFINITE MATRIXBauer and Reinsch (in Wilkinson and Reinsch 1971, p 45) present a very compactalgorithm for inverting a positive definite symmetric matrix in situ, that is,overwriting itself. The principal advantages of this algorithm are as follows.(i) No pivoting is required. This is a consequence of positive definiteness andsymmetry. Peters and Wilkinson (1975) state that this is ‘well known’, but Ibelieve the full analysis is as yet unpublished.(ii) Only a triangular portion of the matrix need be stored due to symmetry,though a working vector of length n, where n is the order of the matrix, is needed.The algorithm is simply the substitution procedure outlined above. Themodifications which are possible due to symmetry and positive definiteness,however, cause the computational steps to look completely different.Consider an intermediate situation in which the first k of the elements x and bhave been exchanged in solvingAx = b(8.6)by the Gauss-Jordan algorithm.
At this stage the matrix of coefficients will havethe formW XY Z(8.7)with W, k by k; X, k by ( n – k); Y, (n – k) by k; and Z, (n – k) by (n – k). Then Wis the inverse of the corresponding block of A since the equations (8.6) are nowgiven by their equivalents(8.8)Thus by setting xj = 0, for j = (k + 1), (k + 2), . . . , n, in both (8.6) and (8.8) therequired association of W and the leading k by k block of A-1 is established.Likewise, by setting bj = 0, for j = 1, . . . , k, in (8.6) and (8.8), Z is the inverse ofthe corresponding block of A-l. (Note that W and Z are symmetric because A andA-1 are symmetric.)From these results, W and Z are both positive definite for all k since A ispositive definite.
This means that the diagonal elements needed as pivots in theGauss-Jordan steps are always positive. (This follows directly from the definitionof positive definiteness for A, that is, that xT A x > 0 for all x 0. )In order to avoid retaining elements above or below the diagonal in a program,we need one more result. This is thatY = –X T(8.9)in matrix (8.7). This can be proved by induction using the Gauss-Jordan98Compact numerical methods for computerssubstitution formulae for step k of the algorithm (i, jk)(8.10)(8.11)(8.12)(8.13)TFor k = 1, therefore, the condition Y = – X is given as(8.14)from equations (8.11) and (8.12).Now assume(8.15)1 < h < k–1, k < j < n, for any of k = 2, .
. . , n. We will show that the hypothesis isthen true for k. By equation (8.13) we have(8.16)where we use the identity(8.17)since these elements belong to a submatrix Z which is symmetric in accord withthe earlier discussion.It remains to establish thatfor j = (k+1), . . . , n(8.18)but this follows immediately from equations (8.11) and (8.12) and the symmetry ofthe submatrix Z. This completes the induction.There is one more trick needed to make the Bauer-Reinsch algorithm extremely compact. This is a sequential cyclic re-ordering of the rows and columnsof A so that the arithmetic is always performed with k = 1.
This re-numerationrelabels (j + 1) as j for j = 1, 2, . . . , (n - 1) and relabels 1 as n. Letting(8.19)this gives a new Gauss-Jordan step(8.20)(8.21)(8.22)(8.23)for i, j = 2, . . . , n.99The symmetric positive definite matrix againA difficulty in this is that the quantities/p have to be stored or they willbe overwritten by Ai-1,j-1 during the work. This is overcome by using a workingvector to store the needed quantities temporarily.Because of the re-numeration we also have the matrix of coefficients in theform(8.24)This permits, via Z = ZT and Y = – X T, the computation involved in (8.23), whereA j1for j < n-kA1j =(8.25)forj > n-k-Ajlwithout recourse to a full matrix. By using row-wise storage of the lower triangleof A, the algorithm below is obtained. Note that after n re-numerations the arrayelements are in the correct order.
Also by counting backwards (step 1) from n to1 the counter will contain (n - k + 1).Algorithm 9. Bauer-Reinsch inversion of a positive definite symmetric matrixprocedure brspdmi(n : integer; {order of problem}var avector : smatvec; {matrix in vector form}var singmat : boolean); {singular matrix flag}{alg09.pas ==Bauer - Reinsch inversion of a symmetric positive definite matrix insitu, Lower triangle of matrix is stored as a vector using rowordering of the original matrix.Copyright 1988 J.C.Nash}vari,j,k,m,q : integer;s,t : real;X : vector;beginwriteln(‘alg09.pas -- Bauer Reinsch inversion’);singmat := false; {initial setting for singularity flag}for k := n downto 1 do {STEP 1}beginif (not singmat) thenbegins := avector[1]; {STEP 2}if s>0.0 then {STEP 3}beginm := 1; {STEP 4}for i := 2 to n do {STEP 5}begin {STEP 6}q := m; m := m+i; t := avector[q+1]; X[i] := -t/s;{Note the use of the temporary vector X[]}if i>k then X[i] := -X[i];{take account of Eqn.
8.22 -- STEP 7}100Compact numerical methods for computersAlgorithm 9. Bauer-Reinsch inversion of a positive definite symmetric matrix (cont.)for j := (q+2) to m do {STEP 8}beginavector[j-i] := avector[j]+t*X[j-q];end; {loop on j}end; {loop on i} {STEP 9}q := q-1; avector[m] := 1.0/s; {STEP 10}for i := 2 to n do avector[q+i] := X[i]; {STEP 11}end {if s>0.0} {STEP 12}elsesingmat := true; {s<=0.0 and we cannot proceed}end; {if (not singmat)}end, {loop on k}end; {alg09.pas == Bauer Reinsch inversion brspdm}This completes the inversion, the original matrix having been overwritten by its inverse.Example 8.1. The behaviour of the Bauer-Reinsch Gauss-Jordan inversionSince the result of the algorithm is identical in form to the input, re-entering theprocedure will compute the inverse of the inverse, which should be the same asthe original matrix.
The following output was obtained by applying the BauerReinsch algorithm in this fashion to the Frank and Moler matrices (see appendix1) on a Data General NOVA having a 23-bit mantissa.NEWLOAD ENHBRTLOAD ENHMT4RUNENHBRG AUG 19 75BAUER REINSCHORDER? 5FRANK MATRIX11 21 2 31 2 3 41 2 3 4 5NEWLOAD ENHBRTLOAD ENHMT5RUNENHBRG AUG 19 75BAUER REINSCHORDER? 5MOLER MATRIX1-1 2-1 0 3-1 0 1 4-1 0 1 2 5INVERSEINVERSEROW 12ROW 2-1 2ROW 30 -1 2ROW 4-1 2ROW 50 0 0 -1ROW86ROW43ROW22ROW12ROW811222311 646 3 254 2 1 1101The symmetric positive definite matrix againINVERSEROW 11ROW 21 2ROW 3.999999ROW 4.999999ROW 51 2 3OF INVERSE2 32 3 44 5INVERSE OF INVERSEROW 1.999994ROW 2-1 2ROW 3-.999987ROW 4-.999989ROW 5-.9999990 2.999970.9999763.999980.9999781.999984.99998Previous HomeChapter 9THE ALGEBRAIC EIGENVALUE PROBLEM9.1.
INTRODUCTIONThe next three chapters are concerned with the solution of algebraic eigenvalue problemsA x = ex(2.62)andAx = eBx.(2.63)The treatment of these problems will be highly selective, and only methods whichI feel are particularly suitable for small computers will be considered. The readerinterested in other methods is referred to Wilkinson (1965) and Wilkinson andReinsch (1971) for an introduction to the very large literature on methods for thealgebraic eigenproblem. Here I have concentrated on providing methods whichare reliable and easy to implement for matrices which can be stored in thecomputer main memory, except for a short discussion in chapter 19 of two methodssuitable for sparse matrices. The methods have also been chosen to fit in withideas already developed so the reader should not feel on too unfamiliar ground.Thus it is possible, even likely, that better methods exist for most applications andany user should be aware of this.
Section 10.5 discusses some of the possiblemethods for real symmetric methods.9.2. THE POWER METHOD AND INVERSE ITERATIONOne of the simplest methods for finding the eigenvalue of largest magnitude of amatrix A with its associated eigenvector is to iterate using the schemeyi = Ax ixi+1 = yi/|| y i||(9.1a)(9.1b)from some starting vector x1 until successive x i are identical. The eigenvector isthen given by x and the magnitude of the eigenvalue by ||y || where || || representsany vector norm. The simplicity of this algorithm and its requirement for only amatrix-vector product makes it especially suitable for small computers, thoughconvergence may at times be very slow.The method works as follows.
Consider the expansion of x 1 in terms of theeigenvectors φ j, j = 1, 2, . . . n, which span the space. This may not be possible fornon-symmetric matrices. However, the algorithm is generally less useful for suchmatrices and will not be considered further in this context. Wilkinson (1965, chap9) has a good discussion of the difficulties. Therefore, returning to the expansion102NextThe algebraic eigenvalue problem103of x 1, we have(9.2)The first iteration of the power method then gives(9.3)where ej is the eigenvalue of A corresponding to φ j. Denoting the reciprocal ofthe norm of yj by Nj, that isNj = || yj||-1(9.4)the vector xi+1 is given by(9.5)The first factor in this expression is simply a normalisation to prevent numbersbecoming too large to work with in the computer.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
